TABLE OF CONTENTS


ABINIT/m_htetra [ Modules ]

[ Top ] [ 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