TABLE OF CONTENTS
- ABINIT/m_htetra
- m_htetra/get_onetetra_blochl
- m_htetra/get_onetetra_lambinvigneron
- m_htetra/get_onetetra_ppart_lv
- m_htetra/get_ontetratra_lambinvigneron_imag
- m_htetra/htetra_blochl_weights
- m_htetra/htetra_blochl_weights_wvals_zinv
- m_htetra/htetra_free
- m_htetra/htetra_get_delta_mask
- m_htetra/htetra_get_ibz
- m_htetra/htetra_get_onewk_wvals
- m_htetra/htetra_get_onewk_wvals_zinv
- m_htetra/htetra_init
- m_htetra/htetra_init_mapping_ibz
- m_htetra/htetra_print
- m_htetra/htetra_t
- m_htetra/htetra_wvals_weights
- m_htetra/htetra_wvals_weights_delta
- m_htetra/sort_4tetra
- m_htetra/t_htetra_bucket
- m_htetra/tetra_get_onewk
- m_numeric_tools/sort_4tetra_int
ABINIT/m_htetra [ Modules ]
NAME
m_htetra
FUNCTION
Module for tetrahedron integration of DOS and similar quantities Uses some functions from a previous implementation by MJV
COPYRIGHT
Copyright (C) 2010-2022 ABINIT group (HM,MJV) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
TODO
1) Test more carefully the case of degenerate tethraedron 2) Add options to get only delta and/or theta ?
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_htetra 28 29 use defs_basis 30 use m_abicore 31 use m_krank 32 use m_xmpi 33 use m_errors 34 35 use m_fstrings, only : sjoin, itoa, ftoa 36 use m_numeric_tools, only : linspace 37 use m_simtet, only : sim0onei, SIM0TWOI 38 39 implicit none 40 41 private
m_htetra/get_onetetra_blochl [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
get_onetetra_blochl
FUNCTION
Private function to calculate the contributions to the weights due to a single tetrahedron. Extracted from get_tetra_weight
SOURCE
1049 pure subroutine get_onetetra_blochl(eig, energies, nene, bcorr, tweight, dweight) 1050 1051 !Arguments ------------------------------------ 1052 !scalars 1053 integer,intent(in) :: nene,bcorr 1054 real(dp) ,intent(in) :: energies(nene) 1055 !arrays 1056 real(dp), intent(out) :: tweight(4, nene) 1057 real(dp), intent(out) :: dweight(4, nene) 1058 real(dp),intent(in) :: eig(4) 1059 1060 !Local variables------------------------------- 1061 integer :: ieps 1062 real(dp) :: cc,cc1,cc2,cc3 1063 real(dp) :: dcc1de,dcc2de,dcc3de,dccde,eps 1064 real(dp) :: e21,e31,e32,e41,e42,e43 1065 real(dp) :: inv_e32,inv_e41,inv_e42,inv_e43,inv_e21,inv_e31 1066 real(dp) :: deleps1,deleps2,deleps3,deleps4 1067 real(dp) :: e1,e2,e3,e4 1068 real(dp) :: invepsum, cc_pre, dccde_pre 1069 real(dp) :: cc1_pre, cc2_pre, cc3_pre 1070 real(dp) :: dccde_tmp 1071 real(dp) :: bcorr_fact 1072 1073 ! ********************************************************************* 1074 1075 ! This is output 1076 tweight = zero; dweight = zero 1077 1078 ! all notations are from Blochl PRB 49 16223 [[cite:Bloechl1994a]] Appendix B 1079 e1 = eig(1) 1080 e2 = eig(2) 1081 e3 = eig(3) 1082 e4 = eig(4) 1083 e21 = e2-e1 1084 e31 = e3-e1 1085 e41 = e4-e1 1086 e32 = e3-e2 1087 e42 = e4-e2 1088 e43 = e4-e3 1089 inv_e21 = zero; if (e21 > tol14) inv_e21 = 1.d0 / e21 1090 inv_e31 = zero; if (e31 > tol14) inv_e31 = 1.d0 / e31 1091 inv_e41 = zero; if (e41 > tol14) inv_e41 = 1.d0 / e41 1092 inv_e32 = zero; if (e32 > tol14) inv_e32 = 1.d0 / e32 1093 inv_e42 = zero; if (e42 > tol14) inv_e42 = 1.d0 / e42 1094 inv_e43 = zero; if (e43 > tol14) inv_e43 = 1.d0 / e43 1095 1096 do ieps=1,nene 1097 eps = energies(ieps) 1098 1099 ! 1100 ! eps < e1 nothing to do 1101 ! 1102 if (eps < e1) cycle 1103 1104 ! 1105 ! e1 < eps < e2 1106 ! 1107 if (eps < e2) then 1108 deleps1 = eps-e1 1109 invepsum = inv_e21+inv_e31+inv_e41 1110 1111 ! Heaviside 1112 cc = inv_e21*inv_e31*inv_e41*deleps1**3 1113 tweight(1,ieps) = cc*(4.d0-deleps1*invepsum) 1114 tweight(2,ieps) = cc*deleps1*inv_e21 1115 tweight(3,ieps) = cc*deleps1*inv_e31 1116 tweight(4,ieps) = cc*deleps1*inv_e41 1117 1118 ! Delta 1119 dccde_pre = 3.d0*inv_e21*inv_e31*inv_e41 1120 dccde = dccde_pre*deleps1**2 1121 dweight(1,ieps) = dccde*(4.d0-deleps1*invepsum)-cc*invepsum 1122 dweight(2,ieps) = (dccde*deleps1+cc) * inv_e21 1123 dweight(3,ieps) = (dccde*deleps1+cc) * inv_e31 1124 dweight(4,ieps) = (dccde*deleps1+cc) * inv_e41 1125 1126 if (bcorr == 1) then 1127 ! bxu, correction terms based on Bloechl's paper 1128 bcorr_fact = 4.d0/40.d0*dccde_pre*deleps1*deleps1 1129 tweight(1,ieps) = tweight(1,ieps) + bcorr_fact*( e21+e31+e41) 1130 tweight(2,ieps) = tweight(2,ieps) + bcorr_fact*(-e21+e32+e42) 1131 tweight(3,ieps) = tweight(3,ieps) + bcorr_fact*(-e31-e32+e43) 1132 tweight(4,ieps) = tweight(4,ieps) + bcorr_fact*(-e41-e42-e43) 1133 1134 bcorr_fact = 8.d0/40.d0*dccde_pre*deleps1 1135 dweight(1,ieps) = dweight(1,ieps) + bcorr_fact*( e21+e31+e41) 1136 dweight(2,ieps) = dweight(2,ieps) + bcorr_fact*(-e21+e32+e42) 1137 dweight(3,ieps) = dweight(3,ieps) + bcorr_fact*(-e31-e32+e43) 1138 dweight(4,ieps) = dweight(4,ieps) + bcorr_fact*(-e41-e42-e43) 1139 end if 1140 1141 cycle 1142 endif 1143 1144 ! 1145 ! e2 < eps < e3 1146 ! 1147 if (eps < e3) then 1148 deleps1 = eps-e1 1149 deleps2 = eps-e2 1150 deleps3 = e3-eps 1151 deleps4 = e4-eps 1152 cc1_pre = inv_e31*inv_e41 1153 cc2_pre = inv_e41*inv_e32*inv_e31 1154 cc3_pre = inv_e42*inv_e32*inv_e41 1155 1156 ! Heaviside 1157 cc1 = cc1_pre*deleps1*deleps1 1158 cc2 = cc2_pre*deleps1*deleps2*deleps3 1159 cc3 = cc3_pre*deleps2*deleps2*deleps4 1160 1161 tweight(1,ieps) = (cc1)+& 1162 (cc1+cc2)*deleps3*inv_e31+& 1163 (cc1+cc2+cc3)*deleps4*inv_e41 1164 tweight(2,ieps) = (cc1+cc2+cc3)+& 1165 (cc2+cc3)*deleps3*inv_e32+& 1166 (cc3)*deleps4*inv_e42 1167 tweight(3,ieps) = (cc1+cc2)*deleps1*inv_e31+& 1168 (cc2+cc3)*deleps2*inv_e32 1169 tweight(4,ieps) = (cc1+cc2+cc3)*deleps1*inv_e41+& 1170 (cc3)*deleps2*inv_e42 1171 1172 ! Delta 1173 dcc1de = cc1_pre*(2.d0*deleps1) 1174 dcc2de = cc2_pre*( -deleps1*deleps2+deleps1*deleps3+deleps2*deleps3) 1175 dcc3de = cc3_pre*(2.d0*deleps2*deleps4-deleps2*deleps2) 1176 dweight(1,ieps) = dcc1de+& 1177 ((dcc1de+dcc2de)*deleps3-(cc1+cc2))*inv_e31+& 1178 ((dcc1de+dcc2de+dcc3de)*deleps4-(cc1+cc2+cc3))*inv_e41 1179 dweight(2,ieps) = (dcc1de+dcc2de+dcc3de)+& 1180 ((dcc2de+dcc3de)*deleps3-(cc2+cc3))*inv_e32+& 1181 (dcc3de*deleps4-cc3)*inv_e42 1182 dweight(3,ieps) = ((dcc1de+dcc2de)*deleps1+(cc1+cc2))*inv_e31+& 1183 ((dcc2de+dcc3de)*deleps2+(cc2+cc3))*inv_e32 1184 dweight(4,ieps) = ((dcc1de+dcc2de+dcc3de)*deleps1+(cc1+cc2+cc3))*inv_e41+& 1185 (dcc3de*deleps2+cc3)*inv_e42 1186 1187 if (bcorr == 1) then 1188 ! bxu, correction terms based on Bloechl's paper 1189 ! The correction terms may cause the dweight become negative 1190 bcorr_fact = 4.d0/40.d0*cc1_pre*(3.d0*e21+6.d0*deleps2-3.d0*(e31+e42)*deleps2*deleps2*inv_e32*inv_e42) 1191 tweight(1,ieps) = tweight(1,ieps) + bcorr_fact*( e21+e31+e41) 1192 tweight(2,ieps) = tweight(2,ieps) + bcorr_fact*(-e21+e32+e42) 1193 tweight(3,ieps) = tweight(3,ieps) + bcorr_fact*(-e31-e32+e43) 1194 tweight(4,ieps) = tweight(4,ieps) + bcorr_fact*(-e41-e42-e43) 1195 1196 bcorr_fact = 4.d0/40.d0*cc1_pre*(6.d0-6.d0*(e31+e42)*deleps2*inv_e32*inv_e42) 1197 dweight(1,ieps) = dweight(1,ieps) + bcorr_fact*( e21+e31+e41) 1198 dweight(2,ieps) = dweight(2,ieps) + bcorr_fact*(-e21+e32+e42) 1199 dweight(3,ieps) = dweight(3,ieps) + bcorr_fact*(-e31-e32+e43) 1200 dweight(4,ieps) = dweight(4,ieps) + bcorr_fact*(-e41-e42-e43) 1201 end if 1202 1203 cycle 1204 endif 1205 1206 ! 1207 ! e3 < eps < e4 1208 ! 1209 if (eps < e4) then 1210 deleps4 = e4-eps 1211 invepsum = inv_e41+inv_e42+inv_e43 1212 1213 ! Heaviside 1214 cc_pre = inv_e41*inv_e42*inv_e43 1215 cc = cc_pre*deleps4**3 1216 tweight(1,ieps) = one - deleps4*cc*inv_e41 1217 tweight(2,ieps) = one - deleps4*cc*inv_e42 1218 tweight(3,ieps) = one - deleps4*cc*inv_e43 1219 tweight(4,ieps) = one - cc*(4.d0-deleps4*invepsum) 1220 1221 ! Delta 1222 dccde = -3.d0*cc_pre*deleps4**2 1223 dccde_tmp = dccde*deleps4 - cc 1224 dweight(1,ieps) = -dccde_tmp * inv_e41 1225 dweight(2,ieps) = -dccde_tmp * inv_e42 1226 dweight(3,ieps) = -dccde_tmp * inv_e43 1227 dweight(4,ieps) = -4.d0*dccde + dccde_tmp*invepsum 1228 1229 if (bcorr == 1) then 1230 ! bxu, correction terms based on Bloechl's paper 1231 ! The correction terms may cause the dweight become negative 1232 bcorr_fact = 12.d0/40.d0*cc_pre*deleps4*deleps4 1233 tweight(1,ieps) = tweight(1,ieps) + bcorr_fact*( e21+e31+e41) 1234 tweight(2,ieps) = tweight(2,ieps) + bcorr_fact*(-e21+e32+e42) 1235 tweight(3,ieps) = tweight(3,ieps) + bcorr_fact*(-e31-e32+e43) 1236 tweight(4,ieps) = tweight(4,ieps) + bcorr_fact*(-e41-e42-e43) 1237 1238 bcorr_fact = - 24.d0/40.d0*cc_pre*deleps4 1239 dweight(1,ieps) = dweight(1,ieps) + bcorr_fact*( e21+e31+e41) 1240 dweight(2,ieps) = dweight(2,ieps) + bcorr_fact*(-e21+e32+e42) 1241 dweight(3,ieps) = dweight(3,ieps) + bcorr_fact*(-e31-e32+e43) 1242 dweight(4,ieps) = dweight(4,ieps) + bcorr_fact*(-e41-e42-e43) 1243 end if 1244 1245 cycle 1246 endif 1247 1248 ! 1249 ! e4 < eps 1250 ! 1251 if (e4 < eps) then 1252 1253 ! Heaviside 1254 tweight(:,ieps:) = one 1255 1256 ! Delta unchanged by this tetrahedron 1257 exit 1258 end if 1259 1260 ! if we have a fully degenerate tetrahedron, 1261 ! 1) the tweight is a Heaviside (step) function, which is correct above, but 1262 ! 2) the dweight should contain a Dirac function 1263 ! 1264 end do 1265 1266 end subroutine get_onetetra_blochl
m_htetra/get_onetetra_lambinvigneron [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
get_onetetra_lambinvigneron
FUNCTION
Compute the complex weights according to: P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984) This routine is adapted from tdep where it was implemented by Olle Hellman, all credits go to him
INPUTS
OUTPUT
SOURCE
1284 !pure 1285 subroutine get_onetetra_lambinvigneron(eig, z, cw) 1286 1287 ! dispersion values at the corners of the tetrahedron 1288 real(dp), intent(in) :: eig(4) 1289 ! energy to evaluate the weights at 1290 complex(dp), intent(in) :: z 1291 ! complex weights 1292 complex(dp), intent(out) :: cw(4) 1293 complex(dp) :: EZ1,EZ2,EZ3,EZ4 1294 real(dp) :: tol = tol14 1295 !real(dp) :: tol = tol10 1296 !real(dp) :: tol = tol6 1297 real(dp) :: Emin,Emax,Zdist 1298 real(dp) :: E12,E13,E14,E23,E24,E34 1299 real(dp) :: a,b,c,d,e,f 1300 complex(dp) zmE(4), verli(4) !, verm(4) 1301 !integer :: ii, jj 1302 1303 cw = zero 1304 1305 ! Min and max energy 1306 Emin=eig(1) 1307 Emax=eig(4) 1308 1309 ! First the complex energy differences 1310 zmE = z - eig(:) 1311 EZ1=z-eig(1) 1312 EZ2=z-eig(2) 1313 EZ3=z-eig(3) 1314 EZ4=z-eig(4) 1315 ! Smallest distance |z-Ei|, to determine wether I should switch to the 1316 ! asymptotic behavior for numerical stability. 1317 Zdist=huge(Zdist) 1318 Zdist=min(Zdist,abs(EZ2)) 1319 Zdist=min(Zdist,abs(EZ3)) 1320 Zdist=min(Zdist,abs(EZ4)) 1321 !@TODO add asymptotic thing with continued fractions 1322 1323 ! Then the energy differences, for the coefficients. Must always be positive, I hope. 1324 E12=eig(2)-eig(1) 1325 E13=eig(3)-eig(1) 1326 E14=eig(4)-eig(1) 1327 E23=eig(3)-eig(2) 1328 E24=eig(4)-eig(2) 1329 E34=eig(4)-eig(3) 1330 a=zero; if ( E12 .gt. tol ) a=one/E12 1331 b=zero; if ( E13 .gt. tol ) b=one/E13 1332 c=zero; if ( E14 .gt. tol ) c=one/E14 1333 d=zero; if ( E23 .gt. tol ) d=one/E23 1334 e=zero; if ( E24 .gt. tol ) e=one/E24 1335 f=zero; if ( E34 .gt. tol ) f=one/E34 1336 1337 ! Now get the actual weights 1338 ! e1=e2=e3=e4 1339 if ( E12+E23+E34 .lt. tol ) then 1340 #if 0 1341 !print *, "e1=e2=e3=e4" 1342 cw(1)=0.25_dp/EZ1 1343 cw(2)=0.25_dp/EZ2 1344 cw(3)=0.25_dp/EZ3 1345 cw(4)=0.25_dp/EZ4 1346 #else 1347 call SIM0TWOI(cw, VERLI, z-eig) 1348 #endif 1349 ! e2=e3=e4 ! diff wrt simteta 1350 elseif ( E23+E34 .lt. tol ) then 1351 #if 0 1352 !print *, "e2=e3=e4" 1353 cw(1)=-a - (3*a**2*EZ2)*half + 3*a**3*EZ1*EZ2 + 3*a**4*EZ1**2*EZ2*Log(EZ2/EZ1) 1354 cw(2)=-a - (3*a**2*EZ2)*half + 3*a**3*EZ1*EZ2 + 3*a**4*EZ1**2*EZ2*Log(EZ2/EZ1) 1355 cw(3)=-a - (3*a**2*EZ2)*half + 3*a**3*EZ1*EZ2 + 3*a**4*EZ1**2*EZ2*Log(EZ2/EZ1) 1356 cw(4)=-a*third + (a**2*EZ1)*half - a**3*EZ1**2 + a**4*EZ1**3*Log(EZ2/EZ1) 1357 #else 1358 call SIM0TWOI(cw, VERLI, z-eig) 1359 #endif 1360 1361 ! e1=e2=e3 1362 elseif ( E12+E23 .lt. tol ) then 1363 #if 0 1364 !print *, "e1=e2=e3" ! diff wrt simteta 1365 cw(1)=f*third - (EZ4*f**2)*half + EZ4**2*f**3 + EZ4**3*f**4*Log(EZ4/EZ3) 1366 cw(2)=f*third - (EZ4*f**2)*half + EZ4**2*f**3 + EZ4**3*f**4*Log(EZ4/EZ3) 1367 cw(3)=f*third - (EZ4*f**2)*half + EZ4**2*f**3 + EZ4**3*f**4*Log(EZ4/EZ3) 1368 cw(4)=-f + (3*EZ3*f**2)*half - 3*EZ3*EZ4*f**3 + 3*EZ3*EZ4**2*f**4*Log(EZ3/EZ4) 1369 #else 1370 call SIM0TWOI(cw, VERLI, z-eig) 1371 #endif 1372 ! e1=e2 e3=e4 1373 elseif ( E12+E34 .lt. tol ) then 1374 #if 0 1375 !print *, "e1=e2 e3=e4" ! FIXME This is Buggy, does not work even with parabolic dispersion 1376 cw(1)=-d - (3*d**2*EZ2)*half + 3*d**3*EZ2*EZ3 + 3*d**4*EZ2*EZ3**2*Log(EZ2/EZ3) 1377 cw(2)=-d - (3*d**2*EZ2)*half + 3*d**3*EZ2*EZ3 + 3*d**4*EZ2*EZ3**2*Log(EZ2/EZ3) 1378 cw(3)=d - (3*d**2*EZ3)*half - 3*d**3*EZ2*EZ3 + 3*d**4*EZ2**2*EZ3*Log(EZ3/EZ2) 1379 cw(4)=d - (3*d**2*EZ3)*half - 3*d**3*EZ2*EZ3 + 3*d**4*EZ2**2*EZ3*Log(EZ3/EZ2) 1380 #else 1381 !cw(1) = nine * EZ3**2 * EZ2 / E23**4 * log(EZ2/EZ3) * EZ2 * (EZ3 -E23)/E23**3 - one/E23 1382 !cw(2) = cw(1) 1383 !cw(3) = nine * EZ2**2 * EZ3 / E23**4 * log(EZ3/EZ2) * EZ3 * (EZ2 +E23)/E23**3 + one/E23 1384 !cw(4) = cw(3) 1385 call SIM0TWOI(cw, VERLI, z-eig) 1386 !cw = zero 1387 #endif 1388 1389 ! e3=e4 1390 elseif ( E34 .lt. tol ) then 1391 !print *, "e3=e4" 1392 #if 0 1393 cw(1)=-(a*b**2*EZ1**2*(-1 + (a*EZ2 + 2*b*EZ3)*Log(EZ1))) + a**2*d**2*EZ2**3*Log(EZ2) - & 1394 b**2*d*EZ3**2*(1 + (2*b*EZ1 + d*EZ2)*Log(EZ3)) 1395 cw(2)=a**2*b**2*EZ1**3*Log(EZ1) - a*d**2*EZ2**2*(1 + (a*EZ1 - 2*d*EZ3)*Log(EZ2)) - & 1396 b*d**2*EZ3**2*(1 + (b*EZ1 + 2*d*EZ2)*Log(EZ3)) 1397 cw(3)=a*b**3*EZ1**3*Log(EZ1) - a*d**3*EZ2**3*Log(EZ2) + b*d*EZ3*(half + b*EZ1 + d*EZ2 + & 1398 (b**2*EZ1**2 + b*d*EZ1*EZ2 + d**2*EZ2**2)*Log(EZ3)) 1399 cw(4)=a*b**3*EZ1**3*Log(EZ1) - a*d**3*EZ2**3*Log(EZ2) + b*d*EZ3*(half + b*EZ1 + d*EZ2 + & 1400 (b**2*EZ1**2 + b*d*EZ1*EZ2 + d**2*EZ2**2)*Log(EZ3)) 1401 #else 1402 call SIM0TWOI(cw, VERLI, z-eig) 1403 #endif 1404 ! e2=e3 1405 elseif ( E23 .lt. tol ) then 1406 #if 0 1407 !print *, "e2=e3" 1408 cw(1)=-(a**2*c*EZ1**2*(-1 + (2*a*EZ2 + c*EZ4)*Log(EZ1))) + & 1409 a**2*e*EZ2**2*(1 + (2*a*EZ1 - e*EZ4)*Log(EZ2)) + c**2*e**2*EZ4**3*Log(EZ4) 1410 cw(2)=a**3*c*EZ1**3*Log(EZ1) - & 1411 a*e*EZ2*(half + a*EZ1 - e*EZ4 + (a**2*EZ1**2 - a*e*EZ1*EZ4 + e**2*EZ4**2)*Log(EZ2)) + c*e**3*EZ4**3*Log(EZ4) 1412 cw(3)=a**3*c*EZ1**3*Log(EZ1) - & 1413 a*e*EZ2*(half + a*EZ1 - e*EZ4 + (a**2*EZ1**2 - a*e*EZ1*EZ4 + e**2*EZ4**2)*Log(EZ2)) + c*e**3*EZ4**3*Log(EZ4) 1414 cw(4)=a**2*c**2*EZ1**3*Log(EZ1) - & 1415 a*e**2*EZ2**2*(1 + (a*EZ1 - 2*e*EZ4)*Log(EZ2)) - c*e**2*EZ4**2*(1 + (c*EZ1 + 2*e*EZ2)*Log(EZ4)) 1416 #else 1417 call SIM0TWOI(cw, VERLI, z-eig) 1418 #endif 1419 1420 ! e1=e2 1421 elseif ( E12 .lt. tol ) then 1422 !print *, "e1=e2" 1423 #if 0 1424 cw(1)=b*c*EZ1*(half - b*EZ3 - c*EZ4 + (b**2*EZ3**2 + b*c*EZ3*EZ4 + c**2*EZ4**2)*Log(EZ1)) - & 1425 b**3*EZ3**3*f*Log(EZ3) + c**3*EZ4**3*f*Log(EZ4) 1426 cw(2)=b*c*EZ1*(half - b*EZ3 - c*EZ4 + (b**2*EZ3**2 + b*c*EZ3*EZ4 + c**2*EZ4**2)*Log(EZ1)) - & 1427 b**3*EZ3**3*f*Log(EZ3) + c**3*EZ4**3*f*Log(EZ4) 1428 cw(3)=-(b**2*c*EZ1**2*(-1 + (2*b*EZ3 + c*EZ4)*Log(EZ1))) + & 1429 b**2*EZ3**2*f*(1 + (2*b*EZ1 - EZ4*f)*Log(EZ3)) + c**2*EZ4**3*f**2*Log(EZ4) 1430 cw(4)=-(b*c**2*EZ1**2*(-1 + (b*EZ3 + 2*c*EZ4)*Log(EZ1))) + & 1431 b**2*EZ3**3*f**2*Log(EZ3) - c**2*EZ4**2*f*(1 + (2*c*EZ1 + EZ3*f)*Log(EZ4)) 1432 #else 1433 call SIM0TWOI(cw, VERLI, z-eig) 1434 #endif 1435 ! e1<e2<e3<e4 1436 else 1437 !print *, "e1<e2<e3<e4" 1438 #if 0 1439 cw(1)=a**2*d*e*EZ2**3*Log(EZ2/EZ1) - b**2*d*EZ3**3*f*Log(EZ3/EZ1) + c*(a*b*EZ1**2 + c*e*EZ4**3*f*Log(EZ4/EZ1)) 1440 cw(2)=a**2*b*c*EZ1**3*Log(EZ1/EZ2) - b*d**2*EZ3**3*f*Log(EZ3/EZ2) + e*(-(a*d*EZ2**2) + c*e*EZ4**3*f*Log(EZ4/EZ2)) 1441 cw(3)=a*b**2*c*EZ1**3*Log(EZ1/EZ3) - a*d**2*e*EZ2**3*Log(EZ2/EZ3) + f*(b*d*EZ3**2 + c*e*EZ4**3*f*Log(EZ4/EZ3)) 1442 cw(4)=a*b*c**2*EZ1**3*Log(EZ1/EZ4) - a*d*e**2*EZ2**3*Log(EZ2/EZ4) + f*(-(c*e*EZ4**2) + b*d*EZ3**3*f*Log(EZ3/EZ4)) 1443 #else 1444 !do ii=1,4 1445 ! cw(ii) = zme(ii) ** 2 / prod_wo(ii) 1446 ! do jj=1,4 1447 ! if (jj == ii) cycle 1448 ! cw(ii) = cw(ii) + (zme(jj)**3 / prod_wo(jj) * log(zme(jj) / zme(ii)) / (eig(ii) - eig(jj))) 1449 ! end do 1450 !end do 1451 1452 call SIM0TWOI(cw, VERLI, z-eig) 1453 #endif 1454 endif 1455 1456 ! HM:check this 1457 !cw = cw * two 1458 1459 contains 1460 pure real(dp) function prod_wo(ii) 1461 integer,intent(in) :: ii 1462 integer :: kk 1463 1464 prod_wo = one 1465 do kk=1,4 1466 if (kk == ii) cycle 1467 prod_wo = prod_wo * (eig(kk) - eig(ii)) 1468 end do 1469 end function prod_wo 1470 1471 end subroutine get_onetetra_lambinvigneron
m_htetra/get_onetetra_ppart_lv [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
get_onetetra_ppart_lv
FUNCTION
Compute the complex weights according to: P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984)
INPUTS
nw wvals: energy to evaluate the weights at eig: eigenvalues at the corners of the tetrahedron
OUTPUT
rwg(nw, 4)
SOURCE
2440 pure subroutine get_onetetra_ppart_lv(nw, wvals, eig, rwg) 2441 2442 integer,intent(in) :: nw 2443 real(dp), intent(in) :: wvals(nw), eig(4) 2444 real(dp), intent(out) :: rwg(nw, 4) 2445 2446 !Local variables------------------------------- 2447 integer :: ii, iw !jj, 2448 real(dp),parameter :: tol = tol14 2449 !real(dp),parameter :: tol = tol20 2450 !real(dp),parameter :: tol = tol30 2451 real(dp) :: D12,D13,D14,D23,D24,D34 2452 real(dp) :: e10, e20, e30, e31, e32, e21 2453 real(dp) :: inv_e10, inv_e21, inv_e20, inv_e30, inv_e31, inv_e32 2454 real(dp) :: E0(nw), E1(nw), E2(nw), E3(nw) 2455 2456 ! ********************************************************************* 2457 2458 ! Then the energy differences, for the coefficients. Must always be positive, I hope. 2459 D12 = eig(2) - eig(1) 2460 D13 = eig(3) - eig(1) 2461 D14 = eig(4) - eig(1) 2462 D23 = eig(3) - eig(2) 2463 D24 = eig(4) - eig(2) 2464 D34 = eig(4) - eig(3) 2465 2466 ! Now get the actual weights 2467 ! Notations 2468 ! eij = e_i - e_j 2469 ! Ej = E - e_j 2470 2471 e10 = huge(one); e20 = huge(one); e30 = huge(one); 2472 e31 = huge(one); e32 = huge(one); e21 = huge(one) 2473 inv_e10 = huge(one); inv_e21 = huge(one); inv_e20 = huge(one) 2474 inv_e30 = huge(one); inv_e31 = huge(one); inv_e32 = huge(one) 2475 2476 E0 = wvals(:) - eig(1) 2477 E1 = wvals(:) - eig(2) 2478 E2 = wvals(:) - eig(3) 2479 E3 = wvals(:) - eig(4) 2480 2481 #if 0 2482 where (abs(E0) < tol1) 2483 E0 = tol1 2484 end where 2485 where (abs(E1) < tol1) 2486 E1 = tol1 2487 end where 2488 where (abs(E2) < tol1) 2489 E2 = tol1 2490 end where 2491 where (abs(E3) < tol1) 2492 E3 = tol1 2493 end where 2494 #endif 2495 2496 ! e1=e2=e3=e4 2497 if (D12 + D23 + D34 < tol) then 2498 do ii=1,4 2499 !rwg(:, ii) = 0.25_dp / (wvals - eig(ii)) 2500 rwg(:, ii) = 0.25_dp / E0 2501 end do 2502 2503 ! e2=e3=e4 2504 else if (D23 + D34 < tol) then 2505 e10 = eig(2) - eig(1); inv_e10 = one / e10 2506 2507 do iw=1,nw 2508 rwg(iw, 1) = & 2509 three * E0(iw)**2 * E1(iw) * inv_e10**4 * log(abs(E1(iw) / E0(iw))) & 2510 + 1.5_dp * E1(iw) * (two * E0(iw) + e10) * inv_e10**3 + inv_e10 2511 end do 2512 2513 do iw=1,nw 2514 rwg(iw, 2) = & 2515 E0(iw) ** 3 * inv_e10**4 * log(abs(E0(iw) / E1(iw))) & 2516 - (six * E0(iw)**2 + three * E0(iw) * e10 + two * e10**2) * inv_e10**3 / six 2517 end do 2518 rwg(:,3) = rwg(:,2) 2519 rwg(:,4) = rwg(:,2) 2520 2521 !rwg = zero 2522 2523 ! e1=e2=e3 2524 else if (D12 + D23 < tol) then 2525 2526 e30 = eig(4) - eig(1); inv_e30 = one / e30 2527 2528 do iw=1,nw 2529 rwg(iw, 1) = & 2530 E3(iw)**3 * inv_e30**4 * log(abs(E3(iw) / E0(iw))) & 2531 + (six * E3(iw)**2 - three * E3(iw) * e30 + two * e30**2) * inv_e30**3 / six 2532 end do 2533 rwg(:,2) = rwg(:,1) 2534 rwg(:,3) = rwg(:,1) 2535 2536 do iw=1,nw 2537 rwg(iw, 4) = & 2538 three * E0(iw) * E3(iw)**2 * inv_e30**4 * log(abs(E0(iw) / E3(iw))) & 2539 - 1.5_dp * E0(iw) * (two * E3(iw) - e30) * inv_e30**3 - inv_e30 2540 end do 2541 2542 !rwg = zero 2543 2544 ! e1=e2 < e3=e4 2545 else if (D12 + D34 < tol) then 2546 2547 e20 = eig(3) - eig(1); inv_e20 = one / e20 2548 2549 do iw=1,nw 2550 rwg(iw, 1) = & 2551 three * E0(iw) * E2(iw)**2 * inv_e20**4 * log(abs(E0(iw) / E2(iw))) & 2552 - 1.5_dp * E0(iw) * (two * E2(iw) - e20) * inv_e20**3 - inv_e20 2553 end do 2554 rwg(:,2) = rwg(:,1) 2555 2556 do iw=1,nw 2557 rwg(iw, 3) = & 2558 three * E0(iw)** 2 * E2(iw) * inv_e20**4 * log(abs(E2(iw) / E0(iw))) & 2559 + 1.5_dp * E2(iw) * (two * E0(iw) + e20) * inv_e20**3 + inv_e20 2560 end do 2561 rwg(:,4) = rwg(:,3) 2562 2563 !rwg = zero 2564 2565 ! e3=e4 2566 else if (D34 < tol) then 2567 2568 e10 = eig(2) - eig(1); inv_e10 = one / e10 2569 e20 = eig(3) - eig(1); inv_e20 = one / e20 2570 e21 = eig(3) - eig(2); inv_e21 = one / e21 2571 2572 do iw=1,nw 2573 rwg(iw, 1) = & 2574 E0(iw)**2 * inv_e20**2 * inv_e10 & 2575 * (one + (-two * E2(iw) * inv_e20 - E1(iw) * inv_e10) * log(abs(E0(iw)))) & 2576 - E2(iw)**2 * inv_e20**2 * inv_e21 & 2577 * (one + (two * E0(iw) * inv_e20 + E1(iw) * inv_e21) * log(abs(E2(iw)))) & 2578 + E1(iw)**3 * inv_e10**2 * inv_e21**2 * log(abs(E1(iw))) 2579 end do 2580 2581 do iw=1,nw 2582 rwg(iw, 2) = & 2583 -E1(iw)**2 * inv_e21**2 * inv_e10 & 2584 * (one + (-two * E2(iw) * inv_e21 + E0(iw) * inv_e10) * log(abs(E1(iw)))) & 2585 - E2(iw)**2 * inv_e21**2 * inv_e20 & 2586 * (one + (two * E1(iw) * inv_e21 + E0(iw) * inv_e20) * log(abs(E2(iw)))) & 2587 + E0(iw)**3 * inv_e10**2 * inv_e20**2 * log(abs(E0(iw))) 2588 end do 2589 2590 do iw=1,nw 2591 rwg(iw, 3) = & 2592 +E0(iw)**3 * inv_e10 * inv_e20**3 * log(abs(E0(iw))) & 2593 -E1(iw)**3 * inv_e10 * inv_e21**3 * log(abs(E1(iw))) & 2594 +E2(iw) * inv_e20 * inv_e21 & 2595 * (half + E0(iw) * inv_e20 + E1(iw) * inv_e21 & 2596 +(E0(iw)**2 * inv_e20**2 + E1(iw)**2 * inv_e21**2 + & 2597 E0(iw) * E1(iw) * inv_e20 * inv_e21) * log(abs(E2(iw)))) 2598 end do 2599 2600 rwg(:,4) = rwg(:,3) 2601 2602 !rwg = zero 2603 2604 ! e2=e3 2605 else if (D23 < tol) then 2606 2607 e10 = eig(2) - eig(1); inv_e10 = one / e10 2608 e30 = eig(4) - eig(1); inv_e30 = one / e30 2609 e31 = eig(4) - eig(2); inv_e31 = one / e31 2610 2611 do iw=1,nw 2612 rwg(iw, 1) = & 2613 E0(iw)**2 * inv_e10**2 * inv_e30 & 2614 * (one - (two * E1(iw) * inv_e10 + E3(iw) * inv_e30) * log(abs(E0(iw)))) & 2615 + E1(iw)**2 * inv_e10**2 * inv_e31 & 2616 * (one + (+two * E0(iw) * inv_e10 - E3(iw) * inv_e31) * log(abs(E1(iw)))) & 2617 + E3(iw)**3 * inv_e30**2 * inv_e31**2 * log(abs(E3(iw))) 2618 end do 2619 2620 do iw=1,nw 2621 rwg(iw, 2) = & 2622 E0(iw)**3 * inv_e30 * inv_e10**3 * log(abs(E0(iw))) & 2623 + E3(iw)**3 * inv_e30 * inv_e31**3 * log(abs(E3(iw))) & 2624 - E1(iw) * inv_e10 * inv_e31 & 2625 * (half + E0(iw) * inv_e10 - E3(iw) * inv_e31 + & 2626 (E0(iw)**2 * inv_e10**2 + E3(iw)**2 * inv_e31**2 - E0(iw) * E3(iw) * inv_e10 * inv_e31) * log(abs(E1(iw)))) 2627 end do 2628 rwg(:,3) = rwg(:,2) 2629 2630 do iw=1,nw 2631 rwg(iw, 4) = & 2632 -E3(iw)**2 * inv_e31**2 * inv_e30 & 2633 * (one + (two * E1(iw) * inv_e31 + E0(iw) * inv_e30) * log(abs(E3(iw)))) & 2634 - E1(iw)**2 * inv_e31**2 * inv_e10 & 2635 * (one + (-two * E3(iw) * inv_e31 + E0(iw) * inv_e10) * log(abs(E1(iw)))) & 2636 + E0(iw)**3 * inv_e30**2 * inv_e10**2 * log(abs(E0(iw))) 2637 end do 2638 2639 !rwg = zero 2640 2641 ! e1=e2 2642 else if (D12 < tol) then 2643 2644 e20 = eig(3) - eig(1); inv_e20 = one / e20 2645 e30 = eig(4) - eig(1); inv_e30 = one / e30 2646 e32 = eig(4) - eig(3); inv_e32 = one / e32 2647 2648 do iw=1,nw 2649 rwg(iw, 1) = & 2650 -E2(iw)**3 * inv_e32 * inv_e20**3 * log(abs(E2(iw))) & 2651 + E3(iw)**3 * inv_e32 * inv_e30**3 * log(abs(E3(iw))) & 2652 + E0(iw) * inv_e20 * inv_e30 & 2653 * (half - E2(iw) * inv_e20 - E3(iw) * inv_e30 + & 2654 (E2(iw)**2 * inv_e20**2 + E3(iw)**2 * inv_e30**2 + E2(iw) * E3(iw) * inv_e20 * inv_e30) * log(abs(E0(iw)))) 2655 end do 2656 rwg(:,2) = rwg(:,1) 2657 2658 do iw=1,nw 2659 rwg(iw, 3) = & 2660 E2(iw)**2 * inv_e20**2 * inv_e32 & 2661 * (one + (two * E0(iw) * inv_e20 - E3(iw) * inv_e32) * log(abs(E2(iw)))) & 2662 + E0(iw)**2 * inv_e20**2 * inv_e30 & 2663 * (one - (two * E2(iw) * inv_e20 + E3(iw) * inv_e30) * log(abs(E0(iw)))) & 2664 + (E3(iw)**3 * inv_e32**2 * inv_e30**2 * log(abs(E3(iw)))) 2665 end do 2666 ! This was wrong due to a misplaced parantes 2667 !rwg(:,3) = zero 2668 2669 do iw=1,nw 2670 rwg(iw, 4) = & 2671 -E3(iw)**2 * inv_e30**2 * inv_e32 & 2672 * (one + (two * E0(iw) * inv_e30 + E2(iw) * inv_e32) * log(abs(E3(iw)))) & 2673 + E0(iw)**2 * inv_e30**2 * inv_e20 & 2674 * (one - (two * E3(iw) * inv_e30 + E2(iw) * inv_e20) * log(abs(E0(iw)))) & 2675 + (E2(iw)**3 * inv_e32**2 * inv_e20**2 * log(abs(E2(iw)))) 2676 end do 2677 2678 !rwg = zero 2679 2680 ! e1<e2<e3<e4 2681 else 2682 2683 e10 = eig(2) - eig(1); inv_e10 = one / e10 2684 e20 = eig(3) - eig(1); inv_e20 = one / e20 2685 e21 = eig(3) - eig(2); inv_e21 = one / e21 2686 e30 = eig(4) - eig(1); inv_e30 = one / e30 2687 e31 = eig(4) - eig(2); inv_e31 = one / e31 2688 e32 = eig(4) - eig(3); inv_e32 = one / e32 2689 2690 do iw=1,nw 2691 rwg(iw, 1) = & 2692 E0(iw)**2 * inv_e10 * inv_e20 * inv_e30 & 2693 * (one - (E1(iw) * inv_e10 + E2(iw) * inv_e20 + E3(iw) * inv_e30) * log(abs(E0(iw)))) & 2694 + E1(iw)**3 * inv_e10**2 * inv_e21 * inv_e31 * log(abs(E1(iw))) & 2695 - E2(iw)**3 * inv_e20**2 * inv_e21 * inv_e32 * log(abs(E2(iw))) & 2696 + E3(iw)**3 * inv_e30**2 * inv_e31 * inv_e32 * log(abs(E3(iw))) 2697 end do 2698 2699 do iw=1,nw 2700 rwg(iw, 2) = & 2701 -E1(iw)**2 * inv_e10 * inv_e21 * inv_e31 & 2702 * (one + (E0(iw) * inv_e10 - E2(iw) * inv_e21 - E3(iw) * inv_e31) * log(abs(E1(iw)))) & 2703 + E0(iw)**3 * inv_e10**2 * inv_e20 * inv_e30 * log(abs(E0(iw))) & 2704 - E2(iw)**3 * inv_e20 * inv_e21**2 * inv_e32 * log(abs(E2(iw))) & 2705 + E3(iw)**3 * inv_e30 * inv_e31**2 * inv_e32 * log(abs(E3(iw))) 2706 end do 2707 2708 do iw=1,nw 2709 rwg(iw, 3) = & 2710 E2(iw)**2 * inv_e20 * inv_e21 * inv_e32 & 2711 * (one + (E0(iw) * inv_e20 + E1(iw) * inv_e21 - E3(iw) * inv_e32) * log(abs(E2(iw)))) & 2712 + E0(iw)**3 * inv_e10 * inv_e20**2 * inv_e30 * log(abs(E0(iw))) & 2713 - E1(iw)**3 * inv_e10 * inv_e21**2 * inv_e31 * log(abs(E1(iw))) & 2714 + E3(iw)**3 * inv_e30 * inv_e31 * inv_e32**2 * log(abs(E3(iw))) 2715 end do 2716 2717 do iw=1,nw 2718 rwg(iw, 4) = & 2719 -E3(iw)**2 * inv_e30 * inv_e31 * inv_e32 & 2720 * (one + (E0(iw) * inv_e30 + E1(iw) * inv_e31 + E2(iw) * inv_e32) * log(abs(E3(iw)))) & 2721 + E0(iw)**3 * inv_e10 * inv_e20 * inv_e30**2 * log(abs(E0(iw))) & 2722 - E1(iw)**3 * inv_e10 * inv_e21 * inv_e31**2 * log(abs(E1(iw))) & 2723 + E2(iw)**3 * inv_e20 * inv_e21 * inv_e32**2 * log(abs(E2(iw))) 2724 end do 2725 !rwg = zero 2726 2727 end if 2728 2729 end subroutine get_onetetra_ppart_lv
m_htetra/get_ontetratra_lambinvigneron_imag [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
get_ontetratra_lambinvigneron_imag
FUNCTION
Compute the complex weights according to: P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984) This routine is adapted from tdep where it was implemented by Olle Hellman, all credits go to him
INPUTS
OUTPUT
SOURCE
1490 pure subroutine get_onetetetra_lambinvigneron_imag(eig, energies, nene, wt) 1491 1492 ! dispersion values at the corners of the tetrahedron 1493 real(dp), intent(in), dimension(4) :: eig 1494 ! number of energies 1495 integer, intent(in) :: nene 1496 ! energy to evaulate the weights at 1497 real(dp), intent(in) :: energies(nene) 1498 ! integration weights 1499 real(dp), intent(out) :: wt(4,nene) 1500 1501 integer :: ie 1502 real(dp) :: z 1503 real(dp) :: EZ1,EZ2,EZ3,EZ4 1504 real(dp) :: Emin,Emax 1505 real(dp) :: E12,E13,E14,E23,E24,E34 1506 real(dp) :: a,b,c,d,e,f,ff0,ff1,ff2,ff3,gg0,gg1,gg2,gg3,hh0,hh1,hh2,hh3,ii0,ii1,ii2,ii3 1507 1508 wt = zero 1509 Emin = eig(1) 1510 Emax = eig(4) 1511 1512 do ie=1,nene 1513 z = energies(ie) 1514 if (z<eig(1)) then ! e<e1<e2<e3<e4 1515 cycle 1516 else if (z .lt. eig(2)) then ! e1<e<e2<e3<e4 1517 EZ1=z-eig(1) 1518 EZ2=z-eig(2) 1519 EZ3=z-eig(3) 1520 EZ4=z-eig(4) 1521 E12=eig(2)-eig(1) 1522 E13=eig(3)-eig(1) 1523 E14=eig(4)-eig(1) 1524 a=one/E12 1525 b=one/E13 1526 c=one/E14 1527 wt(1,ie)=a*b*c*EZ1**2*(-a*EZ2 - b*EZ3 - c*EZ4) 1528 wt(2,ie)=a**2*b*c*EZ1**3 1529 wt(3,ie)=a*b**2*c*EZ1**3 1530 wt(4,ie)=a*b*c**2*EZ1**3 1531 cycle 1532 else if (z .lt. eig(3)) then ! e1<e2<e<e3<e4 1533 EZ1=z-eig(1) 1534 EZ2=z-eig(2) 1535 EZ3=z-eig(3) 1536 EZ4=z-eig(4) 1537 E13=eig(3)-eig(1) 1538 E14=eig(4)-eig(1) 1539 E23=eig(3)-eig(2) 1540 E24=eig(4)-eig(2) 1541 b=one/E13 1542 c=one/E14 1543 d=one/E23 1544 e=one/E24 1545 ff0=-b**2*EZ3 1546 ff2=-c**2*EZ4 1547 gg0=-d**2*EZ3 1548 gg2=-e**2*EZ4 1549 hh0=d**2*EZ2 1550 hh2=b**2*EZ1 1551 ii0=e**2*EZ2 1552 ii2=c**2*EZ1 1553 ff1=-c*d*EZ1*EZ3-d*e*EZ2*EZ3-c*e*EZ1*EZ4 1554 ff3=-b*d*EZ1*EZ3-b*e*EZ1*EZ4-d*e*EZ2*EZ4 1555 gg1=-b*c*EZ1*EZ3-b*e*EZ2*EZ3-c*e*EZ2*EZ4 1556 gg3=-b*d*EZ2*EZ3-b*c*EZ1*EZ4-c*d*EZ2*EZ4 1557 hh1=-b*c*EZ1*EZ3-b*e*EZ2*EZ3-c*e*EZ2*EZ4 1558 hh3=-c*d*EZ1*EZ3-d*e*EZ2*EZ3-c*e*EZ1*EZ4 1559 ii1=-b*d*EZ2*EZ3-b*c*EZ1*EZ4-c*d*EZ2*EZ4 1560 ii3=-b*d*EZ1*EZ3-b*e*EZ1*EZ4-d*e*EZ2*EZ4 1561 wt(1,ie)=half*(ff0*ff1+ff2*ff3) 1562 wt(2,ie)=half*(gg0*gg1+gg2*gg3) 1563 wt(3,ie)=half*(hh0*hh1+hh2*hh3) 1564 wt(4,ie)=half*(ii0*ii1+ii2*ii3) 1565 cycle 1566 else if (z .lt. eig(4)) then ! e1<e2<e3<e<e4 1567 EZ1=z-eig(1) 1568 EZ2=z-eig(2) 1569 EZ3=z-eig(3) 1570 EZ4=z-eig(4) 1571 E14=eig(4)-eig(1) 1572 E24=eig(4)-eig(2) 1573 E34=eig(4)-eig(3) 1574 c=one/E14 1575 e=one/E24 1576 f=one/E34 1577 wt(1,ie)=-(c**2*e*EZ4**3*f) 1578 wt(2,ie)=-(c*e**2*EZ4**3*f) 1579 wt(3,ie)=-(c*e*EZ4**3*f**2) 1580 wt(4,ie)=c*e*EZ4**2*f*(c*EZ1 + e*EZ2 + EZ3*f) 1581 cycle 1582 else 1583 exit 1584 end if 1585 end do 1586 wt = wt*4.0_dp 1587 1588 end subroutine get_onetetetra_lambinvigneron_imag
m_htetra/htetra_blochl_weights [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_blochl_weights
FUNCTION
Emulates the behaviour of the previous tetrahedron implementation. IBZ weights are included.
INPUTS
OUTPUT
SOURCE
2100 subroutine htetra_blochl_weights(tetra, eig_ibz, enemin, enemax, max_occ, nw, nkpt, bcorr, tweight, dweight, comm) 2101 2102 !Arguments ------------------------------------ 2103 !scalars 2104 class(htetra_t), intent(in) :: tetra 2105 integer,intent(in) :: nw,nkpt, bcorr, comm 2106 real(dp) ,intent(in) :: enemax, enemin, max_occ 2107 !arrays 2108 real(dp) ,intent(in) :: eig_ibz(nkpt) 2109 real(dp) ,intent(out) :: dweight(nw,nkpt), tweight(nw,nkpt) 2110 2111 !Local variables------------------------------- 2112 real(dp) :: wvals(nw) 2113 2114 ! ********************************************************************* 2115 2116 wvals = linspace(enemin, enemax, nw) 2117 call htetra_wvals_weights(tetra,eig_ibz,nw,wvals,max_occ,nkpt,bcorr,tweight,dweight,comm) 2118 2119 end subroutine htetra_blochl_weights
m_htetra/htetra_blochl_weights_wvals_zinv [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_blochl_weights_wvals_zinv
FUNCTION
The same as htetra_get_onewk_wvals_zinv but looping over tetrahedra which is more efficient
INPUTS
nz: Number of frequencies zvals(nw): z-values max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2) nkpt=number of irreducible kpoints zinv_opt: 1 for S. Kaprzyk routines, 2 for Lambin-Vigneron. [erange(2)]: if present, weights are computed with a standard quadrature method if real(z) is outside of this interval and with tetra if inside. comm=MPI communicator
OUTPUT
SOURCE
2148 subroutine htetra_weights_wvals_zinv(tetra, eig_ibz, nz, zvals, max_occ, nkpt, zinv_opt, cweight, comm, erange) 2149 2150 !Arguments ------------------------------------ 2151 !scalars 2152 integer,intent(in) :: nz, nkpt, zinv_opt, comm 2153 class(htetra_t), intent(in) :: tetra 2154 real(dp) ,intent(in) :: max_occ 2155 !arrays 2156 real(dp),intent(in) :: eig_ibz(nkpt) 2157 real(dp),optional,intent(in) :: erange(2) 2158 complex(dp),intent(in) :: zvals(nz) 2159 complex(dp),intent(out) :: cweight(nz, nkpt) 2160 2161 !Local variables------------------------------- 2162 !scalars 2163 integer :: ik_ibz, iz, multiplicity, nprocs, my_rank, ierr, ii, jj, kk, esumk 2164 integer :: tetra_count, itetra, isummit, ihash 2165 !arrays 2166 integer :: ind_ibz(4) 2167 real(dp) :: eig(4), my_erange(2) 2168 complex(dp) :: cw(4), verli(4), verm(4), aw(4), bw(4) !, cw_lw(4) 2169 real(dp) :: rwg(nz, 4) 2170 ! ********************************************************************* 2171 2172 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2173 cweight = zero 2174 2175 my_erange = [-huge(one), huge(one)]; if (present(erange)) my_erange = erange 2176 2177 ! For each bucket of tetrahedra 2178 do ihash=1,tetra%nbuckets 2179 if (mod(ihash, nprocs) /= my_rank) cycle 2180 2181 ! For each tetrahedron that belongs to this k-point 2182 tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2) 2183 do itetra=1,tetra_count 2184 2185 ! Get mapping of each summit to eig_ibz 2186 do isummit=1,4 2187 ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit, itetra) 2188 eig(isummit) = eig_ibz(ind_ibz(isummit)) 2189 end do 2190 2191 ! Get multiplicity 2192 multiplicity = tetra%unique_tetra(ihash)%indexes(0, itetra) 2193 2194 ! Sort energies before calling get_onetetra_lambinvigneron 2195 ! SIM0TWOI does not require sorted energies but since we call the routine 2196 ! to fix get_onetetra_lambinvigneron we need to sort here. 2197 call sort_4tetra(eig, ind_ibz) 2198 2199 if (zinv_opt == 1) then 2200 ! Loop over frequencies 2201 do iz=1,nz 2202 2203 ! Get tetrahedron weights 2204 if (real(zvals(iz)) >= my_erange(1) .and. real(zvals(iz)) <= my_erange(2)) then 2205 2206 select case(zinv_opt) 2207 case(1) 2208 verm = zvals(iz) - eig 2209 call SIM0TWOI(cw, VERLI, VERM) 2210 2211 !call get_onetetra_lambinvigneron(eig, zvals(iz), cw_lw) 2212 !if (any(abs(real(cw_lw(:)) / (real(cw(:)))) > 1.1)) then 2213 ! do ierr=1,4 2214 ! !write(std_out, *) "simte vs lw:", cw(ierr), cw_lw(ierr) 2215 ! write(std_out, *) "cw_lw / (simtet)", cw_lw(ierr) / (cw(ierr)) 2216 ! end do 2217 !end if 2218 2219 case(2) 2220 call get_onetetra_lambinvigneron(eig, zvals(iz), cw) 2221 end select 2222 2223 else 2224 ! Use asymptotic expansion of integral for large z. 2225 aw = (eig + sum(eig)) / five 2226 do ii=1,4 2227 bw(ii) = zero 2228 do jj=1,4 2229 if (jj == ii) cycle 2230 esumk = zero 2231 do kk=1,4 2232 if (kk == jj) cycle 2233 esumk = esumk + (eig(kk) - eig(jj)) ** 2 2234 end do 2235 bw(ii) = bw(ii) + three * (eig(jj) - eig(ii)) ** 2 + esumk 2236 end do 2237 end do 2238 bw = bw / 300_dp 2239 cw = one / (zvals(iz) - aw - bw / zvals(iz)) / four 2240 !This for naive integration 2241 !cw = (one / (zvals(iz) - eig)) / four 2242 end if 2243 2244 ! Accumulate contributions 2245 do isummit=1,4 2246 ik_ibz = ind_ibz(isummit) 2247 cweight(iz, ik_ibz) = cweight(iz, ik_ibz) + cw(isummit) * multiplicity * max_occ 2248 end do 2249 2250 end do ! iz 2251 2252 else 2253 call get_onetetra_ppart_lv(nz, real(zvals), eig, rwg) 2254 2255 do iz=1,nz 2256 ! Accumulate contributions 2257 do isummit=1,4 2258 ik_ibz = ind_ibz(isummit) 2259 cweight(iz, ik_ibz) = cweight(iz, ik_ibz) + rwg(iz, isummit) * multiplicity * max_occ 2260 end do 2261 end do 2262 endif 2263 2264 end do ! itetra 2265 end do 2266 2267 ! Rescale weights 2268 select case(tetra%opt) 2269 case(1) 2270 do ik_ibz=1,tetra%nkibz 2271 cweight(:,ik_ibz) = cweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%nkbz / tetra%tetra_total(ik_ibz) 2272 end do 2273 case(2) 2274 cweight = cweight * tetra%vv 2275 end select 2276 2277 call xmpi_sum(cweight, comm, ierr) 2278 2279 end subroutine htetra_weights_wvals_zinv
m_htetra/htetra_free [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_free
FUNCTION
deallocate tetrahedra pointers if needed
SOURCE
1011 subroutine htetra_free(tetra) 1012 1013 class(htetra_t), intent(inout) :: tetra 1014 integer :: ikibz,ihash 1015 1016 ABI_SFREE(tetra%tetra_count) 1017 ABI_SFREE(tetra%tetra_total) 1018 ABI_SFREE(tetra%ibz_multiplicity) 1019 1020 if (allocated(tetra%unique_tetra)) then 1021 do ihash=1,tetra%nbuckets 1022 ABI_SFREE(tetra%unique_tetra(ihash)%indexes) 1023 end do 1024 ABI_FREE(tetra%unique_tetra) 1025 end if 1026 1027 if (allocated(tetra%ibz)) then 1028 do ikibz=1,tetra%nkibz 1029 ABI_SFREE(tetra%ibz(ikibz)%indexes) 1030 end do 1031 ABI_FREE(tetra%ibz) 1032 end if 1033 1034 end subroutine htetra_free
m_htetra/htetra_get_delta_mask [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_get_delta_mask
FUNCTION
Get a mask for the kpoints where the delta is finite
m_htetra/htetra_get_ibz [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_get_ibz
FUNCTION
Get the itetra tetrahedron contributing to the ikibz k-point
SOURCE
949 pure subroutine htetra_get_ibz(tetra, ikibz, itetra, tetra_mibz) 950 951 class(htetra_t), intent(in) :: tetra 952 integer,intent(in) :: ikibz, itetra 953 integer,intent(out) :: tetra_mibz(0:4) 954 integer :: ihash, jtetra 955 956 ihash = tetra%ibz(ikibz)%indexes(1,itetra) 957 jtetra = tetra%ibz(ikibz)%indexes(2,itetra) 958 tetra_mibz = tetra%unique_tetra(ihash)%indexes(:,jtetra) 959 960 end subroutine htetra_get_ibz
m_htetra/htetra_get_onewk_wvals [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_get_onewk_wvals
FUNCTION
Calculate integration weights and their derivatives for a single k-point in the IBZ.
INPUTS
tetra<htetra_t>=Object with tables for tetrahedron method. ik_ibz=Index of the k-point in the IBZ array bcorr=1 to include Blochl correction else 0. nw=number of energies in wvals nibz=number of irreducible kpoints wvals(nw)=Frequency points. eigen_ibz(nkibz)=eigenenergies for each k point
OUTPUT
weights(nw,2) = integration weights for Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function) for a given (band, k-point, spin).
SOURCE
1616 subroutine htetra_get_onewk_wvals(tetra, ik_ibz, opt, nw, wvals, max_occ, nkibz, eig_ibz, weights) 1617 1618 !Arguments ------------------------------------ 1619 !scalars 1620 integer,intent(in) :: ik_ibz,nw,nkibz,opt 1621 real(dp) ,intent(in) :: max_occ 1622 class(htetra_t), intent(inout) :: tetra 1623 !arrays 1624 real(dp),intent(in) :: wvals(nw) 1625 real(dp),intent(in) :: eig_ibz(nkibz) 1626 real(dp),intent(out) :: weights(nw, 2) 1627 1628 !Local variables------------------------------- 1629 !scalars 1630 integer :: itetra,isummit,tetra_count,tetra_total 1631 real(dp) :: tweight 1632 !arrays 1633 integer :: ind_ibz(4),tetra_mibz(0:4) 1634 real(dp) :: eig(4) 1635 real(dp) :: tweight_tmp(4,nw),dweight_tmp(4,nw) 1636 1637 ! ********************************************************************* 1638 1639 weights = zero 1640 ! lazy evaluation of the mapping from k-points to tetrahedra 1641 if (.not.allocated(tetra%ibz)) call htetra_init_mapping_ibz(tetra) 1642 1643 ! For each tetrahedron that belongs to this k-point 1644 tetra_count = tetra%tetra_count(ik_ibz) 1645 tetra_total = tetra%tetra_total(ik_ibz) 1646 do itetra=1,tetra_count 1647 1648 call htetra_get_ibz(tetra, ik_ibz, itetra, tetra_mibz) 1649 tweight = one*tetra_mibz(0) / tetra_total 1650 do isummit=1,4 1651 ! Get mapping of each summit to eig_ibz 1652 ind_ibz(isummit) = tetra_mibz(isummit) 1653 eig(isummit) = eig_ibz(ind_ibz(isummit)) 1654 end do 1655 1656 ! Sort energies before calling get_onetetra_blochl 1657 call sort_4tetra(eig, ind_ibz) 1658 1659 ! HM: Here we should only compute what we will use! 1660 select case (opt) 1661 case(0:1) 1662 call get_onetetra_blochl(eig, wvals, nw, opt, tweight_tmp, dweight_tmp) 1663 case(2) 1664 call get_onetetetra_lambinvigneron_imag(eig, wvals, nw, dweight_tmp) 1665 tweight_tmp = zero 1666 end select 1667 1668 ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz) 1669 do isummit=1,4 1670 if (ind_ibz(isummit) /= ik_ibz) cycle 1671 weights(:,1) = weights(:,1) + dweight_tmp(isummit,:)*tweight*max_occ 1672 weights(:,2) = weights(:,2) + tweight_tmp(isummit,:)*tweight*max_occ 1673 ! HM: This exit is important, avoids summing the same contribution more than once 1674 exit 1675 end do 1676 end do ! itetra 1677 1678 end subroutine htetra_get_onewk_wvals
m_htetra/htetra_get_onewk_wvals_zinv [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_get_onewk_wvals_zinv
FUNCTION
Calculate integration weights for 1/(z-E(k)) for a single k-point in the IBZ. Using either the implementation from: S. Kaprzyk, Computer Physics Communications 183, 347 (2012). or (TODO) P. Lambin and J.P. Vigneron, Phys. Rev. B 29, 3430 (1984).
INPUTS
tetra<htetra_t>=Object with tables for tetrahedron method. ik_ibz=Index of the k-point in the IBZ array bcorr=1 to include Blochl correction else 0. nw=number of energies in wvals nibz=number of irreducible kpoints wvals(nw)=Frequency points. eigen_ibz(nkibz)=eigenenergies for each k point opt: 1 for S. Kaprzyk routines, 2 for Lambin.
OUTPUT
weights(nw,2) = integration weights for Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function) for a given (band, k-point, spin). [erange(2)]: if present, weights are computed with an approximated asyntotic expression if real(z) is outside of this interval and with tetra if inside.
SOURCE
1755 subroutine htetra_get_onewk_wvals_zinv(tetra, ik_ibz, nz, zvals, max_occ, nkibz, eig_ibz, opt, cweights, erange) 1756 1757 !Arguments ------------------------------------ 1758 !scalars 1759 integer,intent(in) :: ik_ibz,nz,nkibz,opt 1760 real(dp) ,intent(in) :: max_occ 1761 class(htetra_t), intent(inout) :: tetra 1762 !arrays 1763 complex(dp),intent(in) :: zvals(nz) 1764 real(dp),optional,intent(in) :: erange(2) 1765 real(dp),intent(in) :: eig_ibz(nkibz) 1766 complex(dp),intent(out) :: cweights(nz) 1767 1768 !Local variables------------------------------- 1769 !scalars 1770 integer :: itetra,isummit,tetra_total,tetra_count,iz 1771 real(dp) :: tweight 1772 !arrays 1773 integer :: ind_ibz(4),tetra_mibz(0:4) 1774 real(dp) :: eig(4), my_erange(2) 1775 complex(dp) :: verm(4), cw(4), verli(4) 1776 ! ********************************************************************* 1777 1778 cweights = zero 1779 ! lazy evaluation of the mapping from k-points to tetrahedra 1780 if (.not.allocated(tetra%ibz)) call htetra_init_mapping_ibz(tetra) 1781 1782 if (all(opt /= [1, 2])) then 1783 ABI_ERROR(sjoin("Invalid opt:", itoa(opt))) 1784 end if 1785 1786 my_erange = [-huge(one), huge(one)]; if (present(erange)) my_erange = erange 1787 1788 ! For each tetrahedron that belongs to this k-point 1789 tetra_count = tetra%tetra_count(ik_ibz) 1790 tetra_total = tetra%tetra_total(ik_ibz) 1791 do itetra=1,tetra_count 1792 1793 call htetra_get_ibz(tetra, ik_ibz, itetra, tetra_mibz) 1794 tweight = one * tetra_mibz(0) / tetra_total 1795 do isummit=1,4 1796 ! Get mapping of each summit to eig_ibz 1797 ind_ibz(isummit) = tetra_mibz(isummit) 1798 eig(isummit) = eig_ibz(ind_ibz(isummit)) 1799 end do 1800 1801 ! Loop over frequencies 1802 do iz=1,nz 1803 1804 if (real(zvals(iz)) >= my_erange(1) .and. real(zvals(iz)) <= my_erange(2)) then 1805 select case(opt) 1806 case(1) 1807 verm = zvals(iz) - eig 1808 call SIM0TWOI(cw, VERLI, VERM) 1809 case(2) 1810 call get_onetetra_lambinvigneron(eig, zvals(iz), cw) 1811 end select 1812 else 1813 cw = (one / (zvals(iz) - eig)) / four !* tetra%vv 1814 end if 1815 1816 do isummit=1,4 1817 if (ind_ibz(isummit) /= ik_ibz) cycle 1818 cweights(iz) = cweights(iz) + cw(isummit) * tweight * max_occ 1819 ! HM: This exit is important, avoids summing the same contribution more than once 1820 exit 1821 end do 1822 end do 1823 end do ! itetra 1824 1825 end subroutine htetra_get_onewk_wvals_zinv
m_htetra/htetra_init [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_init
FUNCTION
get tetrahedra characterized by apexes
INPUTS
bz2ibz(nkpt_fullbz)=indexes of irred kpoints equivalent to kpt_fullbz gprimd(3,3) = reciprocal space vectors klatt(3,3)=reciprocal of lattice vectors for full kpoint grid kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone nkpt_fullbz=number of kpoints in full brillouin zone options=1.generate 24 tetrahedra per k-point faster but gives different results depending on the IBZ, small error for large grids 2.generate tetrahedra on the FBZ and map to IBZ slower but same results for IBZ and FBZ. comm= MPI communicator [opt]= 1 for Togo's version, 2 for Blochl's version (default)
OUTPUT
tetra%ibz(4,24,nkibz)=for each k-point, the indexes in the IBZ tetra%vv = tetrahedron volume divided by full BZ volume
SOURCE
182 subroutine htetra_init(tetra, bz2ibz, gprimd, klatt, kpt_fullbz, nkpt_fullbz, kpt_ibz, nkpt_ibz, & 183 ierr, errorstring, comm, & 184 opt) ! optional 185 186 !Arguments ------------------------------------ 187 !scalars 188 integer,intent(in) :: nkpt_fullbz, nkpt_ibz, comm 189 integer,optional,intent(in) :: opt 190 integer,intent(out) :: ierr 191 character(len=80),intent(out) :: errorstring 192 class(htetra_t),intent(out),target :: tetra 193 !arrays 194 integer,intent(in) :: bz2ibz(nkpt_fullbz) 195 real(dp),intent(in) :: gprimd(3,3),klatt(3,3),kpt_fullbz(3,nkpt_fullbz),kpt_ibz(3,nkpt_ibz) 196 197 !Local variables------------------------------- 198 !scalars 199 !type(octree_t) :: oct 200 type(krank_t) :: krank 201 integer :: ikpt2,isummit,itetra,jtetra 202 integer :: ikibz,ikbz,idiag,ihash,min_idiag,my_rank,nprocs 203 integer :: max_ntetra, ntetra 204 real(dp) :: rcvol,length,min_length 205 !arrays 206 integer,allocatable,target :: indexes(:,:), tetra_hash_count(:) 207 integer :: tetra_ibz(4) 208 integer :: tetra_shifts(3,4,24,4) ! 3 dimensions, 4 summits, 24 tetrahedra, 4 main diagonals 209 integer :: tetra_shifts_6(3,4,6,1) ! 3 dimensions, 4 summits, 6 tetrahedra, 4 main diagonals 210 integer :: main_diagonals(3,4) 211 integer :: tetra_mibz(0:4) 212 real(dp) :: k1(3),k2(3),k3(3),diag(3) 213 214 ! ********************************************************************* 215 216 ! Use the shifts from kpclib developed by Atsushi Togo 217 ! This part is produced by a python script 218 ! This implementation is based on spglib and kpclib by Atsushi Togo 219 ! after a discussion with hin on the APS 2019 where he provided 220 ! details of his implementation. 221 ! Note that we don't use it in production as we found that is approach, although faster than the original 222 ! one proposed by Blochl (and implemented by MJV) does not preserve symmetries that is calculations done on the full BZ 223 ! and the IBZ do not produce the same result. The diff, however, decreases if the sampling is densified. 224 225 tetra_shifts(:, 1, 1,1) = [ 0, 0, 0] 226 tetra_shifts(:, 2, 1,1) = [ 1, 0, 0] 227 tetra_shifts(:, 3, 1,1) = [ 1, 1, 0] 228 tetra_shifts(:, 4, 1,1) = [ 1, 1, 1] 229 tetra_shifts(:, 1, 2,1) = [ 0, 0, 0] 230 tetra_shifts(:, 2, 2,1) = [ 1, 0, 0] 231 tetra_shifts(:, 3, 2,1) = [ 1, 0, 1] 232 tetra_shifts(:, 4, 2,1) = [ 1, 1, 1] 233 tetra_shifts(:, 1, 3,1) = [ 0, 0, 0] 234 tetra_shifts(:, 2, 3,1) = [ 0, 1, 0] 235 tetra_shifts(:, 3, 3,1) = [ 1, 1, 0] 236 tetra_shifts(:, 4, 3,1) = [ 1, 1, 1] 237 tetra_shifts(:, 1, 4,1) = [ 0, 0, 0] 238 tetra_shifts(:, 2, 4,1) = [ 0, 1, 0] 239 tetra_shifts(:, 3, 4,1) = [ 0, 1, 1] 240 tetra_shifts(:, 4, 4,1) = [ 1, 1, 1] 241 tetra_shifts(:, 1, 5,1) = [ 0, 0, 0] 242 tetra_shifts(:, 2, 5,1) = [ 0, 0, 1] 243 tetra_shifts(:, 3, 5,1) = [ 1, 0, 1] 244 tetra_shifts(:, 4, 5,1) = [ 1, 1, 1] 245 tetra_shifts(:, 1, 6,1) = [ 0, 0, 0] 246 tetra_shifts(:, 2, 6,1) = [ 0, 0, 1] 247 tetra_shifts(:, 3, 6,1) = [ 0, 1, 1] 248 tetra_shifts(:, 4, 6,1) = [ 1, 1, 1] 249 tetra_shifts(:, 1, 7,1) = [ 0, 0, 0] 250 tetra_shifts(:, 2, 7,1) = [ 0, 1, 0] 251 tetra_shifts(:, 3, 7,1) = [ 0, 1, 1] 252 tetra_shifts(:, 4, 7,1) = [ -1, 0, 0] 253 tetra_shifts(:, 1, 8,1) = [ 0, 0, 0] 254 tetra_shifts(:, 2, 8,1) = [ 0, 0, 1] 255 tetra_shifts(:, 3, 8,1) = [ 0, 1, 1] 256 tetra_shifts(:, 4, 8,1) = [ -1, 0, 0] 257 tetra_shifts(:, 1, 9,1) = [ 0, 0, 0] 258 tetra_shifts(:, 2, 9,1) = [ 1, 0, 0] 259 tetra_shifts(:, 3, 9,1) = [ 1, 0, 1] 260 tetra_shifts(:, 4, 9,1) = [ 0, -1, 0] 261 tetra_shifts(:, 1,10,1) = [ 0, 0, 0] 262 tetra_shifts(:, 2,10,1) = [ 0, 0, 1] 263 tetra_shifts(:, 3,10,1) = [ 1, 0, 1] 264 tetra_shifts(:, 4,10,1) = [ 0, -1, 0] 265 tetra_shifts(:, 1,11,1) = [ 0, 0, 0] 266 tetra_shifts(:, 2,11,1) = [ 0, 0, 1] 267 tetra_shifts(:, 3,11,1) = [ -1, -1, 0] 268 tetra_shifts(:, 4,11,1) = [ 0, -1, 0] 269 tetra_shifts(:, 1,12,1) = [ 0, 0, 0] 270 tetra_shifts(:, 2,12,1) = [ 0, 0, 1] 271 tetra_shifts(:, 3,12,1) = [ -1, -1, 0] 272 tetra_shifts(:, 4,12,1) = [ -1, 0, 0] 273 tetra_shifts(:, 1,13,1) = [ 0, 0, 0] 274 tetra_shifts(:, 2,13,1) = [ 1, 0, 0] 275 tetra_shifts(:, 3,13,1) = [ 1, 1, 0] 276 tetra_shifts(:, 4,13,1) = [ 0, 0, -1] 277 tetra_shifts(:, 1,14,1) = [ 0, 0, 0] 278 tetra_shifts(:, 2,14,1) = [ 0, 1, 0] 279 tetra_shifts(:, 3,14,1) = [ 1, 1, 0] 280 tetra_shifts(:, 4,14,1) = [ 0, 0, -1] 281 tetra_shifts(:, 1,15,1) = [ 0, 0, 0] 282 tetra_shifts(:, 2,15,1) = [ 0, 1, 0] 283 tetra_shifts(:, 3,15,1) = [ -1, 0, -1] 284 tetra_shifts(:, 4,15,1) = [ 0, 0, -1] 285 tetra_shifts(:, 1,16,1) = [ 0, 0, 0] 286 tetra_shifts(:, 2,16,1) = [ 0, 1, 0] 287 tetra_shifts(:, 3,16,1) = [ -1, 0, -1] 288 tetra_shifts(:, 4,16,1) = [ -1, 0, 0] 289 tetra_shifts(:, 1,17,1) = [ 0, 0, 0] 290 tetra_shifts(:, 2,17,1) = [ 1, 0, 0] 291 tetra_shifts(:, 3,17,1) = [ 0, -1, -1] 292 tetra_shifts(:, 4,17,1) = [ 0, 0, -1] 293 tetra_shifts(:, 1,18,1) = [ 0, 0, 0] 294 tetra_shifts(:, 2,18,1) = [ 1, 0, 0] 295 tetra_shifts(:, 3,18,1) = [ 0, -1, -1] 296 tetra_shifts(:, 4,18,1) = [ 0, -1, 0] 297 tetra_shifts(:, 1,19,1) = [ 0, 0, 0] 298 tetra_shifts(:, 2,19,1) = [ -1, -1, -1] 299 tetra_shifts(:, 3,19,1) = [ 0, -1, -1] 300 tetra_shifts(:, 4,19,1) = [ 0, 0, -1] 301 tetra_shifts(:, 1,20,1) = [ 0, 0, 0] 302 tetra_shifts(:, 2,20,1) = [ -1, -1, -1] 303 tetra_shifts(:, 3,20,1) = [ 0, -1, -1] 304 tetra_shifts(:, 4,20,1) = [ 0, -1, 0] 305 tetra_shifts(:, 1,21,1) = [ 0, 0, 0] 306 tetra_shifts(:, 2,21,1) = [ -1, -1, -1] 307 tetra_shifts(:, 3,21,1) = [ -1, 0, -1] 308 tetra_shifts(:, 4,21,1) = [ 0, 0, -1] 309 tetra_shifts(:, 1,22,1) = [ 0, 0, 0] 310 tetra_shifts(:, 2,22,1) = [ -1, -1, -1] 311 tetra_shifts(:, 3,22,1) = [ -1, 0, -1] 312 tetra_shifts(:, 4,22,1) = [ -1, 0, 0] 313 tetra_shifts(:, 1,23,1) = [ 0, 0, 0] 314 tetra_shifts(:, 2,23,1) = [ -1, -1, -1] 315 tetra_shifts(:, 3,23,1) = [ -1, -1, 0] 316 tetra_shifts(:, 4,23,1) = [ 0, -1, 0] 317 tetra_shifts(:, 1,24,1) = [ 0, 0, 0] 318 tetra_shifts(:, 2,24,1) = [ -1, -1, -1] 319 tetra_shifts(:, 3,24,1) = [ -1, -1, 0] 320 tetra_shifts(:, 4,24,1) = [ -1, 0, 0] 321 tetra_shifts(:, 1, 1,2) = [ 0, 0, 0] 322 tetra_shifts(:, 2, 1,2) = [ 1, 0, 0] 323 tetra_shifts(:, 3, 1,2) = [ 0, 1, 0] 324 tetra_shifts(:, 4, 1,2) = [ 0, 1, 1] 325 tetra_shifts(:, 1, 2,2) = [ 0, 0, 0] 326 tetra_shifts(:, 2, 2,2) = [ 1, 0, 0] 327 tetra_shifts(:, 3, 2,2) = [ 0, 0, 1] 328 tetra_shifts(:, 4, 2,2) = [ 0, 1, 1] 329 tetra_shifts(:, 1, 3,2) = [ 0, 0, 0] 330 tetra_shifts(:, 2, 3,2) = [ -1, 1, 0] 331 tetra_shifts(:, 3, 3,2) = [ -1, 1, 1] 332 tetra_shifts(:, 4, 3,2) = [ -1, 0, 0] 333 tetra_shifts(:, 1, 4,2) = [ 0, 0, 0] 334 tetra_shifts(:, 2, 4,2) = [ -1, 0, 1] 335 tetra_shifts(:, 3, 4,2) = [ -1, 1, 1] 336 tetra_shifts(:, 4, 4,2) = [ -1, 0, 0] 337 tetra_shifts(:, 1, 5,2) = [ 0, 0, 0] 338 tetra_shifts(:, 2, 5,2) = [ -1, 1, 0] 339 tetra_shifts(:, 3, 5,2) = [ 0, 1, 0] 340 tetra_shifts(:, 4, 5,2) = [ -1, 1, 1] 341 tetra_shifts(:, 1, 6,2) = [ 0, 0, 0] 342 tetra_shifts(:, 2, 6,2) = [ 0, 1, 0] 343 tetra_shifts(:, 3, 6,2) = [ -1, 1, 1] 344 tetra_shifts(:, 4, 6,2) = [ 0, 1, 1] 345 tetra_shifts(:, 1, 7,2) = [ 0, 0, 0] 346 tetra_shifts(:, 2, 7,2) = [ -1, 0, 1] 347 tetra_shifts(:, 3, 7,2) = [ 0, 0, 1] 348 tetra_shifts(:, 4, 7,2) = [ -1, 1, 1] 349 tetra_shifts(:, 1, 8,2) = [ 0, 0, 0] 350 tetra_shifts(:, 2, 8,2) = [ 0, 0, 1] 351 tetra_shifts(:, 3, 8,2) = [ -1, 1, 1] 352 tetra_shifts(:, 4, 8,2) = [ 0, 1, 1] 353 tetra_shifts(:, 1, 9,2) = [ 0, 0, 0] 354 tetra_shifts(:, 2, 9,2) = [ 0, 0, 1] 355 tetra_shifts(:, 3, 9,2) = [ 0, -1, 0] 356 tetra_shifts(:, 4, 9,2) = [ 1, -1, 0] 357 tetra_shifts(:, 1,10,2) = [ 0, 0, 0] 358 tetra_shifts(:, 2,10,2) = [ 1, 0, 0] 359 tetra_shifts(:, 3,10,2) = [ 0, 0, 1] 360 tetra_shifts(:, 4,10,2) = [ 1, -1, 0] 361 tetra_shifts(:, 1,11,2) = [ 0, 0, 0] 362 tetra_shifts(:, 2,11,2) = [ -1, 0, 1] 363 tetra_shifts(:, 3,11,2) = [ 0, -1, 0] 364 tetra_shifts(:, 4,11,2) = [ -1, 0, 0] 365 tetra_shifts(:, 1,12,2) = [ 0, 0, 0] 366 tetra_shifts(:, 2,12,2) = [ -1, 0, 1] 367 tetra_shifts(:, 3,12,2) = [ 0, 0, 1] 368 tetra_shifts(:, 4,12,2) = [ 0, -1, 0] 369 tetra_shifts(:, 1,13,2) = [ 0, 0, 0] 370 tetra_shifts(:, 2,13,2) = [ 0, 1, 0] 371 tetra_shifts(:, 3,13,2) = [ 0, 0, -1] 372 tetra_shifts(:, 4,13,2) = [ 1, 0, -1] 373 tetra_shifts(:, 1,14,2) = [ 0, 0, 0] 374 tetra_shifts(:, 2,14,2) = [ 1, 0, 0] 375 tetra_shifts(:, 3,14,2) = [ 0, 1, 0] 376 tetra_shifts(:, 4,14,2) = [ 1, 0, -1] 377 tetra_shifts(:, 1,15,2) = [ 0, 0, 0] 378 tetra_shifts(:, 2,15,2) = [ -1, 1, 0] 379 tetra_shifts(:, 3,15,2) = [ 0, 0, -1] 380 tetra_shifts(:, 4,15,2) = [ -1, 0, 0] 381 tetra_shifts(:, 1,16,2) = [ 0, 0, 0] 382 tetra_shifts(:, 2,16,2) = [ -1, 1, 0] 383 tetra_shifts(:, 3,16,2) = [ 0, 1, 0] 384 tetra_shifts(:, 4,16,2) = [ 0, 0, -1] 385 tetra_shifts(:, 1,17,2) = [ 0, 0, 0] 386 tetra_shifts(:, 2,17,2) = [ 0, -1, -1] 387 tetra_shifts(:, 3,17,2) = [ 1, -1, -1] 388 tetra_shifts(:, 4,17,2) = [ 0, 0, -1] 389 tetra_shifts(:, 1,18,2) = [ 0, 0, 0] 390 tetra_shifts(:, 2,18,2) = [ 0, -1, -1] 391 tetra_shifts(:, 3,18,2) = [ 1, -1, -1] 392 tetra_shifts(:, 4,18,2) = [ 0, -1, 0] 393 tetra_shifts(:, 1,19,2) = [ 0, 0, 0] 394 tetra_shifts(:, 2,19,2) = [ 1, -1, -1] 395 tetra_shifts(:, 3,19,2) = [ 0, 0, -1] 396 tetra_shifts(:, 4,19,2) = [ 1, 0, -1] 397 tetra_shifts(:, 1,20,2) = [ 0, 0, 0] 398 tetra_shifts(:, 2,20,2) = [ 1, 0, 0] 399 tetra_shifts(:, 3,20,2) = [ 1, -1, -1] 400 tetra_shifts(:, 4,20,2) = [ 1, 0, -1] 401 tetra_shifts(:, 1,21,2) = [ 0, 0, 0] 402 tetra_shifts(:, 2,21,2) = [ 1, -1, -1] 403 tetra_shifts(:, 3,21,2) = [ 0, -1, 0] 404 tetra_shifts(:, 4,21,2) = [ 1, -1, 0] 405 tetra_shifts(:, 1,22,2) = [ 0, 0, 0] 406 tetra_shifts(:, 2,22,2) = [ 1, 0, 0] 407 tetra_shifts(:, 3,22,2) = [ 1, -1, -1] 408 tetra_shifts(:, 4,22,2) = [ 1, -1, 0] 409 tetra_shifts(:, 1,23,2) = [ 0, 0, 0] 410 tetra_shifts(:, 2,23,2) = [ 0, -1, -1] 411 tetra_shifts(:, 3,23,2) = [ 0, 0, -1] 412 tetra_shifts(:, 4,23,2) = [ -1, 0, 0] 413 tetra_shifts(:, 1,24,2) = [ 0, 0, 0] 414 tetra_shifts(:, 2,24,2) = [ 0, -1, -1] 415 tetra_shifts(:, 3,24,2) = [ 0, -1, 0] 416 tetra_shifts(:, 4,24,2) = [ -1, 0, 0] 417 tetra_shifts(:, 1, 1,3) = [ 0, 0, 0] 418 tetra_shifts(:, 2, 1,3) = [ 1, 0, 0] 419 tetra_shifts(:, 3, 1,3) = [ 0, 1, 0] 420 tetra_shifts(:, 4, 1,3) = [ 1, 0, 1] 421 tetra_shifts(:, 1, 2,3) = [ 0, 0, 0] 422 tetra_shifts(:, 2, 2,3) = [ 0, 1, 0] 423 tetra_shifts(:, 3, 2,3) = [ 0, 0, 1] 424 tetra_shifts(:, 4, 2,3) = [ 1, 0, 1] 425 tetra_shifts(:, 1, 3,3) = [ 0, 0, 0] 426 tetra_shifts(:, 2, 3,3) = [ -1, 1, 0] 427 tetra_shifts(:, 3, 3,3) = [ 0, 0, 1] 428 tetra_shifts(:, 4, 3,3) = [ -1, 0, 0] 429 tetra_shifts(:, 1, 4,3) = [ 0, 0, 0] 430 tetra_shifts(:, 2, 4,3) = [ -1, 1, 0] 431 tetra_shifts(:, 3, 4,3) = [ 0, 1, 0] 432 tetra_shifts(:, 4, 4,3) = [ 0, 0, 1] 433 tetra_shifts(:, 1, 5,3) = [ 0, 0, 0] 434 tetra_shifts(:, 2, 5,3) = [ 1, -1, 1] 435 tetra_shifts(:, 3, 5,3) = [ 0, -1, 0] 436 tetra_shifts(:, 4, 5,3) = [ 1, -1, 0] 437 tetra_shifts(:, 1, 6,3) = [ 0, 0, 0] 438 tetra_shifts(:, 2, 6,3) = [ 0, -1, 1] 439 tetra_shifts(:, 3, 6,3) = [ 1, -1, 1] 440 tetra_shifts(:, 4, 6,3) = [ 0, -1, 0] 441 tetra_shifts(:, 1, 7,3) = [ 0, 0, 0] 442 tetra_shifts(:, 2, 7,3) = [ 1, 0, 0] 443 tetra_shifts(:, 3, 7,3) = [ 1, -1, 1] 444 tetra_shifts(:, 4, 7,3) = [ 1, -1, 0] 445 tetra_shifts(:, 1, 8,3) = [ 0, 0, 0] 446 tetra_shifts(:, 2, 8,3) = [ 1, 0, 0] 447 tetra_shifts(:, 3, 8,3) = [ 1, -1, 1] 448 tetra_shifts(:, 4, 8,3) = [ 1, 0, 1] 449 tetra_shifts(:, 1, 9,3) = [ 0, 0, 0] 450 tetra_shifts(:, 2, 9,3) = [ 0, -1, 1] 451 tetra_shifts(:, 3, 9,3) = [ 1, -1, 1] 452 tetra_shifts(:, 4, 9,3) = [ 0, 0, 1] 453 tetra_shifts(:, 1,10,3) = [ 0, 0, 0] 454 tetra_shifts(:, 2,10,3) = [ 1, -1, 1] 455 tetra_shifts(:, 3,10,3) = [ 0, 0, 1] 456 tetra_shifts(:, 4,10,3) = [ 1, 0, 1] 457 tetra_shifts(:, 1,11,3) = [ 0, 0, 0] 458 tetra_shifts(:, 2,11,3) = [ 0, -1, 1] 459 tetra_shifts(:, 3,11,3) = [ 0, -1, 0] 460 tetra_shifts(:, 4,11,3) = [ -1, 0, 0] 461 tetra_shifts(:, 1,12,3) = [ 0, 0, 0] 462 tetra_shifts(:, 2,12,3) = [ 0, -1, 1] 463 tetra_shifts(:, 3,12,3) = [ 0, 0, 1] 464 tetra_shifts(:, 4,12,3) = [ -1, 0, 0] 465 tetra_shifts(:, 1,13,3) = [ 0, 0, 0] 466 tetra_shifts(:, 2,13,3) = [ 1, 0, 0] 467 tetra_shifts(:, 3,13,3) = [ 0, 0, -1] 468 tetra_shifts(:, 4,13,3) = [ 0, 1, -1] 469 tetra_shifts(:, 1,14,3) = [ 0, 0, 0] 470 tetra_shifts(:, 2,14,3) = [ 1, 0, 0] 471 tetra_shifts(:, 3,14,3) = [ 0, 1, 0] 472 tetra_shifts(:, 4,14,3) = [ 0, 1, -1] 473 tetra_shifts(:, 1,15,3) = [ 0, 0, 0] 474 tetra_shifts(:, 2,15,3) = [ -1, 0, -1] 475 tetra_shifts(:, 3,15,3) = [ 0, 0, -1] 476 tetra_shifts(:, 4,15,3) = [ -1, 1, -1] 477 tetra_shifts(:, 1,16,3) = [ 0, 0, 0] 478 tetra_shifts(:, 2,16,3) = [ -1, 0, -1] 479 tetra_shifts(:, 3,16,3) = [ -1, 1, -1] 480 tetra_shifts(:, 4,16,3) = [ -1, 0, 0] 481 tetra_shifts(:, 1,17,3) = [ 0, 0, 0] 482 tetra_shifts(:, 2,17,3) = [ 0, 0, -1] 483 tetra_shifts(:, 3,17,3) = [ -1, 1, -1] 484 tetra_shifts(:, 4,17,3) = [ 0, 1, -1] 485 tetra_shifts(:, 1,18,3) = [ 0, 0, 0] 486 tetra_shifts(:, 2,18,3) = [ 0, 1, 0] 487 tetra_shifts(:, 3,18,3) = [ -1, 1, -1] 488 tetra_shifts(:, 4,18,3) = [ 0, 1, -1] 489 tetra_shifts(:, 1,19,3) = [ 0, 0, 0] 490 tetra_shifts(:, 2,19,3) = [ -1, 1, 0] 491 tetra_shifts(:, 3,19,3) = [ -1, 1, -1] 492 tetra_shifts(:, 4,19,3) = [ -1, 0, 0] 493 tetra_shifts(:, 1,20,3) = [ 0, 0, 0] 494 tetra_shifts(:, 2,20,3) = [ -1, 1, 0] 495 tetra_shifts(:, 3,20,3) = [ 0, 1, 0] 496 tetra_shifts(:, 4,20,3) = [ -1, 1, -1] 497 tetra_shifts(:, 1,21,3) = [ 0, 0, 0] 498 tetra_shifts(:, 2,21,3) = [ 0, 0, -1] 499 tetra_shifts(:, 3,21,3) = [ 0, -1, 0] 500 tetra_shifts(:, 4,21,3) = [ 1, -1, 0] 501 tetra_shifts(:, 1,22,3) = [ 0, 0, 0] 502 tetra_shifts(:, 2,22,3) = [ 1, 0, 0] 503 tetra_shifts(:, 3,22,3) = [ 0, 0, -1] 504 tetra_shifts(:, 4,22,3) = [ 1, -1, 0] 505 tetra_shifts(:, 1,23,3) = [ 0, 0, 0] 506 tetra_shifts(:, 2,23,3) = [ -1, 0, -1] 507 tetra_shifts(:, 3,23,3) = [ 0, 0, -1] 508 tetra_shifts(:, 4,23,3) = [ 0, -1, 0] 509 tetra_shifts(:, 1,24,3) = [ 0, 0, 0] 510 tetra_shifts(:, 2,24,3) = [ -1, 0, -1] 511 tetra_shifts(:, 3,24,3) = [ 0, -1, 0] 512 tetra_shifts(:, 4,24,3) = [ -1, 0, 0] 513 tetra_shifts(:, 1, 1,4) = [ 0, 0, 0] 514 tetra_shifts(:, 2, 1,4) = [ 1, 0, 0] 515 tetra_shifts(:, 3, 1,4) = [ 1, 1, 0] 516 tetra_shifts(:, 4, 1,4) = [ 0, 0, 1] 517 tetra_shifts(:, 1, 2,4) = [ 0, 0, 0] 518 tetra_shifts(:, 2, 2,4) = [ 0, 1, 0] 519 tetra_shifts(:, 3, 2,4) = [ 1, 1, 0] 520 tetra_shifts(:, 4, 2,4) = [ 0, 0, 1] 521 tetra_shifts(:, 1, 3,4) = [ 0, 0, 0] 522 tetra_shifts(:, 2, 3,4) = [ 0, 1, 0] 523 tetra_shifts(:, 3, 3,4) = [ -1, 0, 1] 524 tetra_shifts(:, 4, 3,4) = [ -1, 0, 0] 525 tetra_shifts(:, 1, 4,4) = [ 0, 0, 0] 526 tetra_shifts(:, 2, 4,4) = [ 0, 1, 0] 527 tetra_shifts(:, 3, 4,4) = [ -1, 0, 1] 528 tetra_shifts(:, 4, 4,4) = [ 0, 0, 1] 529 tetra_shifts(:, 1, 5,4) = [ 0, 0, 0] 530 tetra_shifts(:, 2, 5,4) = [ 1, 0, 0] 531 tetra_shifts(:, 3, 5,4) = [ 0, -1, 1] 532 tetra_shifts(:, 4, 5,4) = [ 0, -1, 0] 533 tetra_shifts(:, 1, 6,4) = [ 0, 0, 0] 534 tetra_shifts(:, 2, 6,4) = [ 1, 0, 0] 535 tetra_shifts(:, 3, 6,4) = [ 0, -1, 1] 536 tetra_shifts(:, 4, 6,4) = [ 0, 0, 1] 537 tetra_shifts(:, 1, 7,4) = [ 0, 0, 0] 538 tetra_shifts(:, 2, 7,4) = [ -1, -1, 1] 539 tetra_shifts(:, 3, 7,4) = [ -1, -1, 0] 540 tetra_shifts(:, 4, 7,4) = [ 0, -1, 0] 541 tetra_shifts(:, 1, 8,4) = [ 0, 0, 0] 542 tetra_shifts(:, 2, 8,4) = [ -1, -1, 1] 543 tetra_shifts(:, 3, 8,4) = [ -1, -1, 0] 544 tetra_shifts(:, 4, 8,4) = [ -1, 0, 0] 545 tetra_shifts(:, 1, 9,4) = [ 0, 0, 0] 546 tetra_shifts(:, 2, 9,4) = [ -1, -1, 1] 547 tetra_shifts(:, 3, 9,4) = [ 0, -1, 1] 548 tetra_shifts(:, 4, 9,4) = [ 0, -1, 0] 549 tetra_shifts(:, 1,10,4) = [ 0, 0, 0] 550 tetra_shifts(:, 2,10,4) = [ -1, -1, 1] 551 tetra_shifts(:, 3,10,4) = [ -1, 0, 1] 552 tetra_shifts(:, 4,10,4) = [ -1, 0, 0] 553 tetra_shifts(:, 1,11,4) = [ 0, 0, 0] 554 tetra_shifts(:, 2,11,4) = [ -1, -1, 1] 555 tetra_shifts(:, 3,11,4) = [ 0, -1, 1] 556 tetra_shifts(:, 4,11,4) = [ 0, 0, 1] 557 tetra_shifts(:, 1,12,4) = [ 0, 0, 0] 558 tetra_shifts(:, 2,12,4) = [ -1, -1, 1] 559 tetra_shifts(:, 3,12,4) = [ -1, 0, 1] 560 tetra_shifts(:, 4,12,4) = [ 0, 0, 1] 561 tetra_shifts(:, 1,13,4) = [ 0, 0, 0] 562 tetra_shifts(:, 2,13,4) = [ 0, 0, -1] 563 tetra_shifts(:, 3,13,4) = [ 1, 0, -1] 564 tetra_shifts(:, 4,13,4) = [ 1, 1, -1] 565 tetra_shifts(:, 1,14,4) = [ 0, 0, 0] 566 tetra_shifts(:, 2,14,4) = [ 0, 0, -1] 567 tetra_shifts(:, 3,14,4) = [ 0, 1, -1] 568 tetra_shifts(:, 4,14,4) = [ 1, 1, -1] 569 tetra_shifts(:, 1,15,4) = [ 0, 0, 0] 570 tetra_shifts(:, 2,15,4) = [ 1, 0, 0] 571 tetra_shifts(:, 3,15,4) = [ 1, 0, -1] 572 tetra_shifts(:, 4,15,4) = [ 1, 1, -1] 573 tetra_shifts(:, 1,16,4) = [ 0, 0, 0] 574 tetra_shifts(:, 2,16,4) = [ 0, 1, 0] 575 tetra_shifts(:, 3,16,4) = [ 0, 1, -1] 576 tetra_shifts(:, 4,16,4) = [ 1, 1, -1] 577 tetra_shifts(:, 1,17,4) = [ 0, 0, 0] 578 tetra_shifts(:, 2,17,4) = [ 1, 0, 0] 579 tetra_shifts(:, 3,17,4) = [ 1, 1, 0] 580 tetra_shifts(:, 4,17,4) = [ 1, 1, -1] 581 tetra_shifts(:, 1,18,4) = [ 0, 0, 0] 582 tetra_shifts(:, 2,18,4) = [ 0, 1, 0] 583 tetra_shifts(:, 3,18,4) = [ 1, 1, 0] 584 tetra_shifts(:, 4,18,4) = [ 1, 1, -1] 585 tetra_shifts(:, 1,19,4) = [ 0, 0, 0] 586 tetra_shifts(:, 2,19,4) = [ 0, 0, -1] 587 tetra_shifts(:, 3,19,4) = [ 0, 1, -1] 588 tetra_shifts(:, 4,19,4) = [ -1, 0, 0] 589 tetra_shifts(:, 1,20,4) = [ 0, 0, 0] 590 tetra_shifts(:, 2,20,4) = [ 0, 1, 0] 591 tetra_shifts(:, 3,20,4) = [ 0, 1, -1] 592 tetra_shifts(:, 4,20,4) = [ -1, 0, 0] 593 tetra_shifts(:, 1,21,4) = [ 0, 0, 0] 594 tetra_shifts(:, 2,21,4) = [ 0, 0, -1] 595 tetra_shifts(:, 3,21,4) = [ 1, 0, -1] 596 tetra_shifts(:, 4,21,4) = [ 0, -1, 0] 597 tetra_shifts(:, 1,22,4) = [ 0, 0, 0] 598 tetra_shifts(:, 2,22,4) = [ 1, 0, 0] 599 tetra_shifts(:, 3,22,4) = [ 1, 0, -1] 600 tetra_shifts(:, 4,22,4) = [ 0, -1, 0] 601 tetra_shifts(:, 1,23,4) = [ 0, 0, 0] 602 tetra_shifts(:, 2,23,4) = [ 0, 0, -1] 603 tetra_shifts(:, 3,23,4) = [ -1, -1, 0] 604 tetra_shifts(:, 4,23,4) = [ 0, -1, 0] 605 tetra_shifts(:, 1,24,4) = [ 0, 0, 0] 606 tetra_shifts(:, 2,24,4) = [ 0, 0, -1] 607 tetra_shifts(:, 3,24,4) = [ -1, -1, 0] 608 tetra_shifts(:, 4,24,4) = [ -1, 0, 0] 609 610 ! These shifts are taken from previous tetrahedron implementation by MJV and BXU 611 ! TODO: implement shifts for the other diagonals 612 tetra_shifts_6(:,1,1,1) = [0,0,0] 613 tetra_shifts_6(:,2,1,1) = [1,0,0] 614 tetra_shifts_6(:,3,1,1) = [0,1,0] 615 tetra_shifts_6(:,4,1,1) = [1,0,1] 616 tetra_shifts_6(:,1,2,1) = [1,0,0] 617 tetra_shifts_6(:,2,2,1) = [1,1,0] 618 tetra_shifts_6(:,3,2,1) = [0,1,0] 619 tetra_shifts_6(:,4,2,1) = [1,0,1] 620 tetra_shifts_6(:,1,3,1) = [0,1,0] 621 tetra_shifts_6(:,2,3,1) = [1,1,0] 622 tetra_shifts_6(:,3,3,1) = [1,0,1] 623 tetra_shifts_6(:,4,3,1) = [1,1,1] 624 tetra_shifts_6(:,1,4,1) = [0,0,0] 625 tetra_shifts_6(:,2,4,1) = [0,1,0] 626 tetra_shifts_6(:,3,4,1) = [0,0,1] 627 tetra_shifts_6(:,4,4,1) = [1,0,1] 628 tetra_shifts_6(:,1,5,1) = [0,0,1] 629 tetra_shifts_6(:,2,5,1) = [1,0,1] 630 tetra_shifts_6(:,3,5,1) = [0,1,0] 631 tetra_shifts_6(:,4,5,1) = [0,1,1] 632 tetra_shifts_6(:,1,6,1) = [0,1,0] 633 tetra_shifts_6(:,2,6,1) = [1,0,1] 634 tetra_shifts_6(:,3,6,1) = [0,1,1] 635 tetra_shifts_6(:,4,6,1) = [1,1,1] 636 637 main_diagonals(:,1) = [ 1, 1, 1] ! 0-7 638 main_diagonals(:,2) = [-1, 1, 1] ! 1-6 639 main_diagonals(:,3) = [ 1,-1, 1] ! 2-5 640 main_diagonals(:,4) = [ 1, 1,-1] ! 3-4 641 642 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 643 tetra%nkibz = nkpt_ibz 644 tetra%nkbz = nkpt_fullbz 645 ! HM: this value should be important for performance 646 ! more buckets means faster queries for unique tetrahedra 647 ! but more memory due to the initial size TETRA_SIZE of the buckets 648 ! when changing the number of buckets one should also change the hash function 649 ! to distribute the tetrahedra in the buckets as uniformly as possible 650 ! the simplest hash (not the best!) is: 651 ! ihash = mod(sum(tetra_ibz),nbuckts) 652 ! the value of sum(tetra_ibz) is between 1 and 4*nkibz so I use nkibz nbuckets 653 ! a larger number of buckets should speed up finding the irreducible tetrahedra 654 ! but with more memory allocated 655 tetra%nbuckets = nkpt_ibz 656 ierr = 0 657 ABI_CALLOC(tetra_hash_count,(tetra%nbuckets)) 658 659 ! Determine the smallest diagonal in k-space 660 min_length = huge(min_length) 661 do idiag = 1,4 662 diag(:) = gprimd(:,1)*main_diagonals(1,idiag)+& 663 gprimd(:,2)*main_diagonals(2,idiag)+& 664 gprimd(:,3)*main_diagonals(3,idiag) 665 length = sqrt(diag(1)*diag(1) + diag(2)*diag(2) + diag(3)*diag(3)) 666 if (length < min_length) then 667 min_length = length 668 min_idiag = idiag 669 end if 670 end do 671 672 ! HM TODO: Avoid krank and map the k-point grid to indexes 673 ! Make full k-point rank arrays 674 !oct = octree_init(kpt_fullbz,2**4,[-one,-one,-one],[two,two,two]) 675 krank = krank_new(nkpt_fullbz, kpt_fullbz) 676 677 ! 678 ! HM (13/04/2019): I implement two different versions: 679 ! 1. I only use 24 tetrahedra around the IBZ k-point 680 ! following the approach of A. Togo (phonopy, spglib, kspclib). 681 ! 2. I generate tetrahedra on the full Brillouin zone 682 ! and keep track of how many are contributing to the IBZ point and the multiplicities, 683 ! these can be more than 24 tetrahedra. 684 ! This is equivalent to what Matthieu implemented but avoids large memory allocations. 685 ! 686 ! The two implementations differ specially when using low k-point sampling 687 ! (the second yields the same results using IBZ or FBZ, the first one not). 688 ! For large sampling the two approaches yield similar results, with the first 689 ! one using less memory, faster to generate and compute 690 ! 691 ABI_MALLOC(tetra%unique_tetra,(tetra%nbuckets)) 692 do ihash=1,tetra%nbuckets 693 ABI_CALLOC(tetra%unique_tetra(ihash)%indexes,(0:4,TETRA_SIZE)) 694 end do 695 tetra%opt = 2; if (present(opt)) tetra%opt = opt 696 697 select case(tetra%opt) 698 case(1) 699 ! For each k-point in the IBZ store 24 tetrahedra each refering to 4 k-points 700 do ikibz=1,tetra%nkibz 701 !if (mod(ikibz,nprocs) /= my_rank) cycle 702 k1 = kpt_ibz(:,ikibz) 703 tetra_loop1: do itetra=1,24 704 do isummit=1,4 705 ! Find the index of the neighbouring k-points in the BZ 706 k2 = k1 + tetra_shifts(1,isummit,itetra,min_idiag)*klatt(:,1) + & 707 tetra_shifts(2,isummit,itetra,min_idiag)*klatt(:,2) + & 708 tetra_shifts(3,isummit,itetra,min_idiag)*klatt(:,3) 709 ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full ! 710 !ikpt2 = octree_find(oct,k2,dist) 711 !ikpt2 = octree_find_nearest_pbc(oct,k2,dist,shift) 712 !if (dist>tol12) call exit(1) 713 ikpt2 = krank%get_index(k2) 714 ! Find the index of those points in the BZ and IBZ 715 tetra_ibz(isummit) = bz2ibz(ikpt2) 716 end do 717 ! Sort index of irr k-point edges (need this so the comparison works) 718 call sort_4tetra_int(tetra_ibz) 719 720 ! Store only unique tetrahedra 721 ! Compute a very simple hash for each tetrahedron 722 ihash = compute_hash(tetra,tetra_ibz) !mod(sum(tetra_ibz),tetra%nbuckets)+1 723 ! Loop over all tetrahedrons that contain this ikibz as first element 724 do jtetra=1,tetra_hash_count(ihash) 725 ! if tetrahedron already exists add multiplicity 726 if (tetra%unique_tetra(ihash)%indexes(1,jtetra)/=tetra_ibz(1)) cycle 727 if (tetra%unique_tetra(ihash)%indexes(2,jtetra)/=tetra_ibz(2)) cycle 728 if (tetra%unique_tetra(ihash)%indexes(3,jtetra)/=tetra_ibz(3)) cycle 729 if (tetra%unique_tetra(ihash)%indexes(4,jtetra)/=tetra_ibz(4)) cycle 730 tetra%unique_tetra(ihash)%indexes(0,jtetra) = tetra%unique_tetra(ihash)%indexes(0,jtetra)+1 731 cycle tetra_loop1 732 end do 733 ! Otherwise store new tetrahedron 734 tetra_hash_count(ihash) = tetra_hash_count(ihash)+1 735 max_ntetra = size(tetra%unique_tetra(ihash)%indexes,2) 736 ! The contents don't fit the array so I have to resize it 737 if (tetra_hash_count(ihash)>max_ntetra) then 738 ABI_MALLOC(indexes,(0:4,max_ntetra+TETRA_STEP)) 739 indexes(0:4,:max_ntetra) = tetra%unique_tetra(ihash)%indexes 740 indexes(:,max_ntetra+1:) = 0 741 ABI_MOVE_ALLOC(indexes,tetra%unique_tetra(ihash)%indexes) 742 end if 743 tetra%unique_tetra(ihash)%indexes(1:,tetra_hash_count(ihash)) = tetra_ibz(:) 744 tetra%unique_tetra(ihash)%indexes(0, tetra_hash_count(ihash)) = 1 745 end do tetra_loop1 746 end do 747 ! HM: The multiplicity of the tetrahedrons computed so far is wrong because we are using IBZ 748 ! I compute the k-point multiplicity so I can fix this later on 749 ! Only needed in blochl_weights* interface for looping over tetrahedra. 750 ! in the onewk routines the weight is known outside 751 ABI_CALLOC(tetra%ibz_multiplicity,(tetra%nkibz)) 752 do ikbz=1,nkpt_fullbz 753 ikibz = bz2ibz(ikbz) 754 tetra%ibz_multiplicity(ikibz) = tetra%ibz_multiplicity(ikibz) + 1 755 end do 756 757 case(2) 758 min_idiag = 1 759 ! For each k-point in the BZ generate the 6 tetrahedra that tesselate a microzone 760 do ikbz=1,tetra%nkbz 761 k1 = kpt_fullbz(:,ikbz) 762 tetra_loop2: do itetra=1,6 763 ! Determine tetrahedron 764 do isummit=1,4 765 ! Find the index of the neighbouring k-points in the BZ 766 k2 = k1 + tetra_shifts_6(1,isummit,itetra,min_idiag)*klatt(:,1) + & 767 tetra_shifts_6(2,isummit,itetra,min_idiag)*klatt(:,2) + & 768 tetra_shifts_6(3,isummit,itetra,min_idiag)*klatt(:,3) 769 ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full ! 770 !ikpt2 = octree_find(oct,k2,dist) 771 !ikpt2 = octree_find_nearest_pbc(oct,k2,dist,shift) 772 !if (dist>tol12) call exit(1) 773 ikpt2 = krank%get_index(k2) 774 ! Find the index of those points in the BZ and IBZ 775 tetra_ibz(isummit) = bz2ibz(ikpt2) 776 end do 777 ! Sort index of irr k-point edges (need this so the comparison works) 778 call sort_4tetra_int(tetra_ibz) 779 780 ! Store only unique tetrahedra 781 ! Compute a very simple hash for each tetrahedron 782 ihash = compute_hash(tetra,tetra_ibz) !mod(sum(tetra_ibz),tetra%nbuckets)+1 783 ! Loop over all tetrahedrons that contain this ikibz as first element 784 do jtetra=1,tetra_hash_count(ihash) 785 ! if tetrahedron already exists add multiplicity 786 if (tetra%unique_tetra(ihash)%indexes(1,jtetra)/=tetra_ibz(1)) cycle 787 if (tetra%unique_tetra(ihash)%indexes(2,jtetra)/=tetra_ibz(2)) cycle 788 if (tetra%unique_tetra(ihash)%indexes(3,jtetra)/=tetra_ibz(3)) cycle 789 if (tetra%unique_tetra(ihash)%indexes(4,jtetra)/=tetra_ibz(4)) cycle 790 tetra%unique_tetra(ihash)%indexes(0,jtetra) = tetra%unique_tetra(ihash)%indexes(0,jtetra)+1 791 cycle tetra_loop2 792 end do 793 ! Otherwise store new tetrahedron 794 tetra_hash_count(ihash) = tetra_hash_count(ihash)+1 795 max_ntetra = size(tetra%unique_tetra(ihash)%indexes,2) 796 ! The contents don't fit the array so I have to resize it 797 if (tetra_hash_count(ihash)>max_ntetra) then 798 ABI_MALLOC(indexes,(0:4,max_ntetra+TETRA_STEP)) 799 indexes(0:4,:max_ntetra) = tetra%unique_tetra(ihash)%indexes 800 indexes(:,max_ntetra+1:) = 0 801 ABI_MOVE_ALLOC(indexes,tetra%unique_tetra(ihash)%indexes) 802 end if 803 tetra%unique_tetra(ihash)%indexes(1:,tetra_hash_count(ihash)) = tetra_ibz(:) 804 tetra%unique_tetra(ihash)%indexes(0, tetra_hash_count(ihash)) = 1 805 end do tetra_loop2 806 end do 807 case default 808 ierr = 1 809 write(errorstring,*) 'Invalid option for the generation of tetrahedra,',ch10,& 810 'possible options are:',ch10,& 811 '1. Generate 24 tetrahedra per k-point',ch10,& 812 '2. Generate tetrahedra in the FBZ a map to IBZ (default)' 813 return 814 end select 815 816 !ierr = octree_free(oct) 817 ABI_FREE(tetra_hash_count) 818 call krank%free() 819 820 ! Do some maintenance: free unused memory and count unique tetrahedra per IBZ point 821 tetra%nunique_tetra = 0 822 do ihash=1,tetra%nbuckets 823 ! Count tetrahedra in this bucket 824 ntetra = count(tetra%unique_tetra(ihash)%indexes(0,:)>0) 825 tetra%nunique_tetra = tetra%nunique_tetra + ntetra 826 ! Allocate array with right size 827 ABI_MALLOC(indexes,(0:4,ntetra)) 828 indexes = tetra%unique_tetra(ihash)%indexes(:,:ntetra) 829 ABI_MOVE_ALLOC(indexes, tetra%unique_tetra(ihash)%indexes) 830 end do 831 832 ! Sum the multiplicity 833 ABI_MALLOC(tetra%tetra_count,(tetra%nkibz)) 834 ABI_MALLOC(tetra%tetra_total,(tetra%nkibz)) 835 tetra%tetra_count = 0 836 tetra%tetra_total = 0 837 do ihash=1,tetra%nbuckets 838 ntetra = size(tetra%unique_tetra(ihash)%indexes,2) 839 do itetra=1,ntetra 840 tetra_mibz = tetra%unique_tetra(ihash)%indexes(:,itetra) 841 do isummit=1,4 842 ikibz = tetra_mibz(isummit) 843 tetra%tetra_total(ikibz) = tetra%tetra_total(ikibz) + tetra_mibz(0) 844 tetra%tetra_count(ikibz) = tetra%tetra_count(ikibz) + 1 845 end do 846 end do 847 end do 848 tetra%nibz_tetra = sum(tetra%tetra_count) 849 850 ! HM: This was being allocated here, however this is only used when we loop over kpoints 851 ! I will only allocate this memory if the htetra_get_onewk_* routines are called (lazy evaluation) 852 !call htetra_init_mapping_ibz(tetra) 853 854 ! Calculate the volume of the tetrahedra 855 rcvol = abs(gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3))- & 856 gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3))+ & 857 gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3))) 858 859 ! Volume of all tetrahedra should be the same as that of tetra 1 860 ! this is the volume of 1 tetrahedron, should be coherent with notation in Lehmann & Taut 861 k1(:) = gprimd(:,1)*klatt(1,1) + gprimd(:,2)*klatt(2,1) + gprimd(:,3)*klatt(3,1) 862 k2(:) = gprimd(:,1)*klatt(1,2) + gprimd(:,2)*klatt(2,2) + gprimd(:,3)*klatt(3,2) 863 k3(:) = gprimd(:,1)*klatt(1,3) + gprimd(:,2)*klatt(2,3) + gprimd(:,3)*klatt(3,3) 864 tetra%vv = abs(k1(1)*(k2(2)*k3(3)-k2(3)*k3(2))- & 865 k1(2)*(k2(1)*k3(3)-k2(3)*k3(1))+ & 866 k1(3)*(k2(1)*k3(2)-k2(2)*k3(1))) / 6.d0 / rcvol 867 868 contains 869 integer function compute_hash(tetra,t) result(ihash) 870 class(htetra_t),intent(in) :: tetra 871 integer,intent(in) :: t(4) 872 ihash = mod(sum(t), tetra%nbuckets) + 1 873 ! TODO: should use a more general hash function that supports more buckets 874 ! Something like: 875 ! id = t(1)*nk3+t(2)*nk2+t(3)*nk1+t(4) 876 ! where nk is the number of points in the IBZ 877 ! Computing this leads to overflow so should use 878 ! mod comutation operations 879 ! (A + B) mod C = (A mod C + B mod C) mod C 880 ! (A * B) mod C = (A mod C * B mod C) mod C 881 ! A^B mod C = ( (A mod C)^B ) mod C 882 end function compute_hash 883 884 end subroutine htetra_init
m_htetra/htetra_init_mapping_ibz [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_init_mapping_ibz
FUNCTION
The mapping to the IBZ has its own allocation routine. I will only allocate this memory if the htetra_get_onewk_* routines are called (lazy evaluation)
SOURCE
899 subroutine htetra_init_mapping_ibz(tetra) 900 901 class(htetra_t),intent(inout) :: tetra 902 integer :: ikibz, itetra, isummit, ihash, ntetra 903 integer :: tetra_count(tetra%nkibz),tetra_mibz(0:4) 904 905 real(dp) :: mem_mb 906 907 ! Only execute the following if not yet alocated 908 if (allocated(tetra%ibz)) return 909 910 ! Allocate IBZ to tetrahedron mapping 911 ABI_MALLOC(tetra%ibz, (tetra%nkibz)) 912 mem_mb = ABI_MEM_MB(tetra%ibz) 913 do ikibz=1,tetra%nkibz 914 ABI_MALLOC(tetra%ibz(ikibz)%indexes, (2, tetra%tetra_count(ikibz))) 915 mem_mb = mem_mb + 2 * tetra%tetra_count(ikibz) * 4 * b2Mb 916 end do 917 918 call wrtout(std_out, sjoin(" Allocating tetra%ibz%indexes with memory:", ftoa(mem_mb, fmt="f8.1"), " [Mb] <<< MEM")) 919 920 ! Create mapping from IBZ to unique tetrahedra 921 tetra_count = 0 922 do ihash=1,tetra%nbuckets 923 ntetra = size(tetra%unique_tetra(ihash)%indexes, dim=2) 924 do itetra=1,ntetra 925 tetra_mibz = tetra%unique_tetra(ihash)%indexes(:,itetra) 926 do isummit=1,4 927 ikibz = tetra_mibz(isummit) 928 tetra_count(ikibz) = tetra_count(ikibz) + 1 929 tetra%ibz(ikibz)%indexes(1, tetra_count(ikibz)) = ihash 930 tetra%ibz(ikibz)%indexes(2, tetra_count(ikibz)) = itetra 931 end do 932 end do 933 end do 934 935 end subroutine htetra_init_mapping_ibz
m_htetra/htetra_print [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_print
FUNCTION
write information about the tetrahedra object
SOURCE
974 subroutine htetra_print(self, unit) 975 976 class(htetra_t), intent(in) :: self 977 integer,intent(in) :: unit 978 979 real(dp) :: total_size, unique_tetra_size, ibz_pointer_size 980 981 if (unit == dev_null) return 982 983 unique_tetra_size = self%nunique_tetra * 5* four / 1024 ** 2 984 total_size = unique_tetra_size 985 !write(unit,'(a,i0)') ' htetra unique_tetra:', self%nunique_tetra 986 !write(unit,'(a,f12.1,a)') ' htetra unique_tetra_size ', unique_tetra_size, ' [Mb] <<< MEM' 987 if (allocated(self%ibz)) then 988 ibz_pointer_size = self%nibz_tetra*2*four / 1024 ** 2 989 write(unit,'(a,i0)') ' ibz_tetra: ', self%nibz_tetra 990 write(unit,'(a,f12.1,a)') ' ibz_tetra_size: ', ibz_pointer_size, ' [Mb] <<< MEM' 991 total_size = total_size + ibz_pointer_size 992 end if 993 ! integer arrays 994 total_size = total_size + 3 * (self%nkibz * 4) / 1024**2 995 write(unit,'(a,f12.1,a)') ' htetra total size:', total_size, ' [Mb] <<< MEM' 996 997 end subroutine htetra_print
m_htetra/htetra_t [ Types ]
[ Top ] [ m_htetra ] [ Types ]
NAME
htetra_t
FUNCTION
tetrahedron geometry object
SOURCE
72 type, public :: htetra_t 73 74 integer :: opt 75 ! Option for the generation of tetrahedra 76 77 integer :: nkibz 78 ! Number of points in the irreducible Brillouin zone 79 80 integer :: nkbz 81 ! Number of points in the full Brillouin zone 82 83 integer :: nbuckets 84 ! Number of buckets for the hash table 85 86 integer :: nunique_tetra 87 ! Number of unique tetrahedron 88 89 integer :: nibz_tetra 90 ! Number of ibz tetrahedron 91 92 integer,allocatable :: tetra_total(:) 93 ! (%nbkibz) 94 ! Equivalent tetrahedra per kpoint (number of irred tetra times multiplicity) 95 96 integer,allocatable :: tetra_count(:) 97 ! (%nbkibz) 98 ! Inequivalent tetrahedra per kpoint (number of irred tetra) 99 100 integer,allocatable :: ibz_multiplicity(:) 101 ! (%nbkibz) 102 ! Multiplicity of each k-point 103 104 real(dp) :: vv 105 ! volume of the tetrahedra 106 107 real(dp) :: klatt(3, 3) 108 ! reciprocal of lattice vectors for full kpoint grid 109 110 type(htetra_bucket),allocatable :: ibz(:) 111 ! indexes of the tetrahedra for each k-point 112 113 type(htetra_bucket),allocatable :: unique_tetra(:) 114 ! indexes of the unique tetrahedra 115 116 contains 117 118 procedure :: free => htetra_free 119 ! Free memory 120 121 procedure :: print => htetra_print 122 ! Print information about tetrahedron object 123 124 procedure :: get_onewk => htetra_get_onewk 125 ! Calculate integration weights and their derivatives for a single k-point in the IBZ. 126 127 procedure :: get_onewk_wvals => htetra_get_onewk_wvals 128 ! Similar to tetra_get_onewk_wvals but receives arbitrary list of frequency points. 129 130 procedure :: get_onewk_wvals_zinv => htetra_get_onewk_wvals_zinv 131 ! Calculate integration weights for 1/(z-E(k)) for a single k-point in the IBZ. 132 133 procedure :: weights_wvals_zinv => htetra_weights_wvals_zinv 134 ! Same as above but return the weight on all the kpoints by looping over tetrahedra 135 136 procedure :: wvals_weights => htetra_wvals_weights 137 ! Compute delta and theta on a list of energies for all kpoints 138 139 procedure :: wvals_weights_delta => htetra_wvals_weights_delta 140 ! Compute delta on a list of energies for all kpoints 141 142 procedure :: blochl_weights => htetra_blochl_weights 143 ! And interface to help to facilitate the transition to the new tetrahedron implementation 144 145 end type htetra_t
m_htetra/htetra_wvals_weights [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_wvals_weights
FUNCTION
Emulates the behaviour of the previous tetrahedron implementation but taking a list of energies as input. HM: I find that in many routines its better to change the implementation and accumulate the tetrahedron weights in the same way as the gaussian smearing weights using htetra_get_onewk_wvals. However this requires some refactoring of the code. I provide this routine to make it easier to transition to the new tetrahedron implementation without refactoring. Looping over tetrahedra (i.e. using tetra_blochl_weights) is currently faster than looping over k-points. MG: Note, however, that tetra_blochl_weights requires more memory as one has to allocate dweight(nw,nkpt),tweight(nw,nkpt) and the size of the arrays increases quickly with the k-mesh and the number of frequencies (propto mband)
INPUTS
OUTPUT
SOURCE
1919 subroutine htetra_wvals_weights(tetra, eig_ibz, nw, wvals, max_occ, nkpt, opt, tweight, dweight, comm) 1920 1921 !Arguments ------------------------------------ 1922 !scalars 1923 integer,intent(in) :: nw,nkpt,opt,comm 1924 class(htetra_t), intent(in) :: tetra 1925 real(dp) ,intent(in) :: max_occ 1926 !arrays 1927 real(dp),intent(in) :: eig_ibz(nkpt) 1928 real(dp),intent(out) :: dweight(nw,nkpt),tweight(nw,nkpt) 1929 1930 !Local variables------------------------------- 1931 !scalars 1932 integer :: ik_ibz,multiplicity,nprocs,my_rank,ierr 1933 integer :: tetra_count, itetra, isummit, ihash 1934 !arrays 1935 integer :: ind_ibz(4) 1936 real(dp) :: eig(4) 1937 real(dp) :: wvals(nw) 1938 real(dp) :: dweight_tmp(4,nw),tweight_tmp(4,nw) 1939 1940 ! ********************************************************************* 1941 1942 tweight = zero; dweight = zero 1943 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1944 1945 ! For each bucket of tetrahedra 1946 do ihash=1,tetra%nbuckets 1947 if (mod(ihash, nprocs) /= my_rank) cycle 1948 1949 ! For each tetrahedron 1950 tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2) 1951 do itetra=1,tetra_count 1952 1953 ! Get mapping of each summit to eig_ibz 1954 do isummit=1,4 1955 ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit, itetra) 1956 eig(isummit) = eig_ibz(ind_ibz(isummit)) 1957 end do 1958 1959 ! Sort energies before calling get_onetetra_blochl 1960 call sort_4tetra(eig, ind_ibz) 1961 1962 ! Get tetrahedron weights 1963 select case (opt) 1964 case(0:1) 1965 call get_onetetra_blochl(eig, wvals, nw, opt, tweight_tmp, dweight_tmp) 1966 case(2) 1967 call get_onetetetra_lambinvigneron_imag(eig, wvals, nw, dweight_tmp) 1968 tweight_tmp = zero 1969 end select 1970 1971 ! Acumulate the contributions 1972 multiplicity = tetra%unique_tetra(ihash)%indexes(0, itetra) 1973 do isummit=1,4 1974 ik_ibz = ind_ibz(isummit) 1975 dweight(:,ik_ibz) = dweight(:,ik_ibz) + dweight_tmp(isummit,:) * multiplicity * max_occ 1976 tweight(:,ik_ibz) = tweight(:,ik_ibz) + tweight_tmp(isummit,:) * multiplicity * max_occ 1977 end do 1978 end do ! itetra 1979 end do 1980 1981 ! Rescale weights 1982 select case(tetra%opt) 1983 case(1) 1984 do ik_ibz=1,tetra%nkibz 1985 dweight(:,ik_ibz) = dweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%tetra_total(ik_ibz) / tetra%nkbz 1986 tweight(:,ik_ibz) = tweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%tetra_total(ik_ibz) / tetra%nkbz 1987 end do 1988 case(2) 1989 dweight = dweight*tetra%vv / 4.0_dp 1990 tweight = tweight*tetra%vv / 4.0_dp 1991 end select 1992 1993 call xmpi_sum(dweight, comm, ierr) 1994 call xmpi_sum(tweight, comm, ierr) 1995 1996 end subroutine htetra_wvals_weights
m_htetra/htetra_wvals_weights_delta [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
htetra_wvals_weights_delta
FUNCTION
Same as above but computing only delta for performance and memory HM: Should find a clean way to avoid copy paste routine
SOURCE
2011 subroutine htetra_wvals_weights_delta(tetra, eig_ibz, nw, wvals, max_occ, nkpt, opt, dweight, comm) 2012 2013 !Arguments ------------------------------------ 2014 !scalars 2015 integer,intent(in) :: nw, nkpt, opt, comm 2016 class(htetra_t), intent(in) :: tetra 2017 real(dp),intent(in) :: max_occ 2018 !arrays 2019 real(dp),intent(in) :: eig_ibz(nkpt), wvals(nw) 2020 real(dp),intent(out) :: dweight(nw, nkpt) 2021 2022 !Local variables------------------------------- 2023 !scalars 2024 integer :: ik_ibz,multiplicity,nprocs,my_rank,ierr 2025 integer :: tetra_count, itetra, isummit, ihash 2026 !arrays 2027 integer :: ind_ibz(4) 2028 real(dp) :: eig(4), dweight_tmp(4,nw),tweight_tmp(4,nw) 2029 2030 ! ********************************************************************* 2031 2032 dweight = zero 2033 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2034 2035 ! For each bucket of tetrahedra 2036 do ihash=1,tetra%nbuckets 2037 if (mod(ihash, nprocs) /= my_rank) cycle 2038 2039 ! For each tetrahedron 2040 tetra_count = size(tetra%unique_tetra(ihash)%indexes, dim=2) 2041 do itetra=1,tetra_count 2042 2043 ! Get mapping of each summit to eig_ibz 2044 do isummit=1,4 2045 ind_ibz(isummit) = tetra%unique_tetra(ihash)%indexes(isummit,itetra) 2046 eig(isummit) = eig_ibz(ind_ibz(isummit)) 2047 end do 2048 2049 ! Sort energies before calling get_onetetra_blochl 2050 call sort_4tetra(eig, ind_ibz) 2051 2052 ! Get tetrahedron weights 2053 select case (opt) 2054 case(0:1) 2055 call get_onetetra_blochl(eig, wvals, nw, opt, tweight_tmp, dweight_tmp) 2056 case(2) 2057 call get_onetetetra_lambinvigneron_imag(eig, wvals, nw, dweight_tmp) 2058 end select 2059 2060 ! Acumulate the contributions 2061 multiplicity = tetra%unique_tetra(ihash)%indexes(0,itetra) 2062 do isummit=1,4 2063 ik_ibz = ind_ibz(isummit) 2064 dweight(:,ik_ibz) = dweight(:,ik_ibz) + dweight_tmp(isummit,:)*multiplicity*max_occ 2065 end do 2066 end do ! itetra 2067 end do 2068 2069 ! Rescale weights 2070 select case(tetra%opt) 2071 case(1) 2072 do ik_ibz=1,tetra%nkibz 2073 dweight(:,ik_ibz) = dweight(:,ik_ibz) * tetra%ibz_multiplicity(ik_ibz) / tetra%tetra_total(ik_ibz) / tetra%nkbz 2074 end do 2075 case(2) 2076 dweight = dweight * tetra%vv / 4.0_dp 2077 end select 2078 2079 call xmpi_sum(dweight, comm, ierr) 2080 2081 end subroutine htetra_wvals_weights_delta
m_htetra/sort_4tetra [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
sort_4tetra
FUNCTION
Sort double precision array list(4) into ascending numerical order while making corresponding rearrangement of the integer array iperm. Taken from: https://stackoverflow.com/questions/6145364/sort-4-number-with-few-comparisons
INPUTS
list(4) intent(inout) list of double precision numbers to be sorted perm(4) intent(inout) iperm(i)=i (very important)
OUTPUT
list(4) sorted list perm(4) index of permutation given the right ascending order
SOURCE
2304 pure subroutine sort_4tetra(list, perm) 2305 2306 integer, intent(inout) :: perm(4) 2307 real(dp), intent(inout) :: list(4) 2308 2309 !Local variables------------------------------- 2310 integer :: ia,ib,ic,id,ilow1,ilow2,ihigh1,ihigh2 2311 integer :: ilowest,ihighest,imiddle1,imiddle2 2312 real(dp) :: va,vb,vc,vd,vlow1,vlow2,vhigh1,vhigh2 2313 real(dp) :: vlowest,vhighest,vmiddle1,vmiddle2 2314 2315 va = list(1); ia = perm(1) 2316 vb = list(2); ib = perm(2) 2317 vc = list(3); ic = perm(3) 2318 vd = list(4); id = perm(4) 2319 2320 if (va < vb) then 2321 vlow1 = va; vhigh1 = vb 2322 ilow1 = ia; ihigh1 = ib 2323 else 2324 vlow1 = vb; vhigh1 = va 2325 ilow1 = ib; ihigh1 = ia 2326 endif 2327 2328 if (vc < vd) then 2329 vlow2 = vc; vhigh2 = vd 2330 ilow2 = ic; ihigh2 = id 2331 else 2332 vlow2 = vd; vhigh2 = vc 2333 ilow2 = id; ihigh2 = ic 2334 endif 2335 2336 if (vlow1 < vlow2) then 2337 vlowest = vlow1; vmiddle1 = vlow2 2338 ilowest = ilow1; imiddle1 = ilow2 2339 else 2340 vlowest = vlow2; vmiddle1 = vlow1 2341 ilowest = ilow2; imiddle1 = ilow1 2342 endif 2343 2344 if (vhigh1 > vhigh2) then 2345 vhighest = vhigh1; vmiddle2 = vhigh2 2346 ihighest = ihigh1; imiddle2 = ihigh2 2347 else 2348 vhighest = vhigh2; vmiddle2 = vhigh1 2349 ihighest = ihigh2; imiddle2 = ihigh1 2350 endif 2351 2352 if (vmiddle1 < vmiddle2) then 2353 list = [vlowest, vmiddle1, vmiddle2, vhighest] 2354 perm = [ilowest, imiddle1, imiddle2, ihighest] 2355 else 2356 list = [vlowest, vmiddle2, vmiddle1, vhighest] 2357 perm = [ilowest, imiddle2, imiddle1, ihighest] 2358 endif 2359 2360 end subroutine sort_4tetra
m_htetra/t_htetra_bucket [ Types ]
[ Top ] [ m_htetra ] [ Types ]
NAME
t_htetra_bucket
FUNCTION
Store a bunch of tetrahedra
SOURCE
56 type :: htetra_bucket 57 58 integer,allocatable :: indexes(:,:) 59 60 end type htetra_bucket
m_htetra/tetra_get_onewk [ Functions ]
[ Top ] [ m_htetra ] [ Functions ]
NAME
tetra_get_onewk
FUNCTION
Calculate integration weights and their derivatives for a single k-point in the IBZ. Same as above but different calling arguments. IBZ Weights are not included HM: The above is prefered but I keep this one to ease the transition
INPUTS
OUTPUT
SOURCE
1699 subroutine htetra_get_onewk(tetra, ik_ibz, bcorr, nw, nkibz, eig_ibz, enemin, enemax, max_occ, weights) 1700 1701 !Arguments ------------------------------------ 1702 !scalars 1703 integer,intent(in) :: ik_ibz,nw,nkibz,bcorr 1704 class(htetra_t), intent(inout) :: tetra 1705 real(dp) ,intent(in) :: enemin,enemax,max_occ 1706 !arrays 1707 real(dp),intent(in) :: eig_ibz(nkibz) 1708 real(dp),intent(out) :: weights(nw,2) 1709 1710 !Local variables------------------------------- 1711 !scalars 1712 real(dp) :: wvals(nw) 1713 1714 ! ********************************************************************* 1715 1716 weights = zero 1717 wvals = linspace(enemin, enemax, nw) 1718 call htetra_get_onewk_wvals(tetra, ik_ibz, bcorr, nw, wvals, max_occ, nkibz, eig_ibz, weights) 1719 1720 end subroutine htetra_get_onewk
m_numeric_tools/sort_4tetra_int [ Functions ]
[ Top ] [ m_numeric_tools ] [ Functions ]
NAME
sort_4tetra_int
FUNCTION
INPUTS
OUTPUT
SOURCE
2377 pure subroutine sort_4tetra_int(list) 2378 2379 integer, intent(inout) :: list(4) 2380 2381 !Local variables------------------------------- 2382 integer :: va,vb,vc,vd, vlow1,vlow2,vhigh1,vhigh2 2383 integer :: vlowest,vhighest, vmiddle1,vmiddle2 2384 2385 va = list(1) 2386 vb = list(2) 2387 vc = list(3) 2388 vd = list(4) 2389 2390 if (va < vb) then 2391 vlow1 = va; vhigh1 = vb 2392 else 2393 vlow1 = vb; vhigh1 = va 2394 endif 2395 2396 if (vc < vd) then 2397 vlow2 = vc; vhigh2 = vd 2398 else 2399 vlow2 = vd; vhigh2 = vc 2400 endif 2401 2402 if (vlow1 < vlow2) then 2403 vlowest = vlow1; vmiddle1 = vlow2 2404 else 2405 vlowest = vlow2; vmiddle1 = vlow1 2406 endif 2407 2408 if (vhigh1 > vhigh2) then 2409 vhighest = vhigh1; vmiddle2 = vhigh2 2410 else 2411 vhighest = vhigh2; vmiddle2 = vhigh1 2412 endif 2413 2414 if (vmiddle1 < vmiddle2) then 2415 list = [vlowest, vmiddle1, vmiddle2, vhighest] 2416 else 2417 list = [vlowest, vmiddle2, vmiddle1, vhighest] 2418 endif 2419 2420 end subroutine sort_4tetra_int