TABLE OF CONTENTS


ABINIT/m_tetrahedron [ Modules ]

[ Top ] [ Modules ]

NAME

 m_tetrahedron

FUNCTION

  module for tetrahedron interpolation of DOS and similar quantities
  depends on sort_tetra and on m_kpt_rank

COPYRIGHT

  Copyright (C) 2010-2018 ABINIT group (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 carefull the case of degenerate tethraedron (add
  2) Change API so that we can pass the energy mesh instead of omega_min and omega_max

PARENTS

CHILDREN

SOURCE

26 #if defined HAVE_CONFIG_H
27 #include "config.h"
28 #endif
29 
30 #include "libtetra.h"
31 
32 module m_tetrahedron
33 
34  USE_MEMORY_PROFILING
35  USE_MSG_HANDLING
36  use m_kptrank
37 #ifdef HAVE_LIBTETRA_ABINIT
38  use m_io_tools, only : open_file
39 #endif
40 #if defined HAVE_MPI2
41  use mpi
42 #endif
43 
44 implicit none
45 
46 #if defined HAVE_MPI1
47  include 'mpif.h'
48 #endif
49 
50 private

m_tetrahedron/destroy_tetra [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 destroy_tetra

FUNCTION

 deallocate tetrahedra pointers if needed

PARENTS

      ep_el_weights,ep_fs_weights,ep_ph_weights,gstate,m_ebands,m_epjdos
      m_fstab,m_gruneisen,m_phgamma,m_phonons,thmeig,wfk_analyze

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

133 subroutine destroy_tetra (tetra)
134 
135 
136 !This section has been created automatically by the script Abilint (TD).
137 !Do not modify the following lines by hand.
138 #undef ABI_FUNC
139 #define ABI_FUNC 'destroy_tetra'
140 !End of the abilint section
141 
142  implicit none
143 
144  type(t_tetrahedron), intent(inout) :: tetra
145 
146  if (allocated(tetra%tetra_full))  then
147    TETRA_DEALLOCATE(tetra%tetra_full)
148  end if
149  if (allocated(tetra%tetra_mult))  then
150    TETRA_DEALLOCATE(tetra%tetra_mult)
151  end if
152  if (allocated(tetra%tetra_wrap))  then
153    TETRA_DEALLOCATE(tetra%tetra_wrap)
154  end if
155 
156 end subroutine destroy_tetra

m_tetrahedron/get_dbl_tetra_weight [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 get_dbl_tetra_weight

FUNCTION

 calculate integration weights and their derivatives
 for double tetrahedron method from Allen Phys Stat Sol B 120 529 (1983) [[cite:Allen1983b]]
 the k-points and tetrahedra must be the same for both grids, of course,
 but the range of energies is arbitrary

 Omega is called eigen1 here
 E is called eigen2 here
 indexing goes from 1 to 4 for the tetrahedron corners, in order of increasing eigen1
  in Allen, from 0 to 3...

INPUTS

 eigen1_in(nkpt)=eigenenergies for each k point
 eigen2_in(nkpt)=eigenenergies for each k point
 enemin1=minimal energy for DOS in energy 1
 enemax1=maximal energy for DOS
 enemin2=minimal energy for DOS in energy 2
 enemax2=maximal energy for DOS
 max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)
 nene1=number of energies for DOS in energy 1
 nene2=number of energies for DOS in energy 2
 nkpt=number of irreducible kpoints
 tetra%ntetra=number of tetra
 tetra%tetra_full(4,2,ntetra)=for each irred tetrahedron, the list of k point vertices
   1 -> irred kpoint   2 -> fullkpt
 tetra%tetra_mult(ntetra)=for each irred tetrahedron, its multiplicity
 tetra%vv = ratio of volume of one tetrahedron in reciprocal space to full BZ volume
 ierr = error code on exit

OUTPUT

  tweight(nkpt,nene1,nene2) = integration weights for each irred kpoint from all adjacent tetrahedra
  dtweightde(nkpt,nene1,nene2) = derivative of tweight wrt energy

PARENTS

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

 833 subroutine get_dbl_tetra_weight(eigen1_in,eigen2_in,enemin1,enemax1,enemin2,enemax2,&
 834 &    max_occ,nene1,nene2,nkpt,tetra,tweight,dtweightde, ierr)
 835 
 836 
 837 !This section has been created automatically by the script Abilint (TD).
 838 !Do not modify the following lines by hand.
 839 #undef ABI_FUNC
 840 #define ABI_FUNC 'get_dbl_tetra_weight'
 841 !End of the abilint section
 842 
 843  implicit none
 844 
 845 !Arguments ------------------------------------
 846 !scalars
 847  integer,intent(in) :: nene1,nene2,nkpt
 848  integer,intent(out) :: ierr
 849  type(t_tetrahedron), intent(in) :: tetra
 850  real(dp_),intent(in) :: enemax1,enemin1
 851  real(dp_),intent(in) :: enemax2,enemin2
 852  real(dp_),intent(in) :: max_occ
 853 !arrays
 854  real(dp_),intent(in) :: eigen1_in(nkpt)
 855  real(dp_),intent(in) :: eigen2_in(nkpt)
 856  real(dp_),intent(out) :: dtweightde(nkpt,nene1,nene2),tweight(nkpt,nene1,nene2)
 857 
 858 !Local variables-------------------------------
 859 !  needed for gaussian replacement of Dirac functions
 860 !  the three coefficients of the DOS as quadratic form,
 861 !    in the interval [eig(ikpt-1), eig(ikpt)]
 862 !    for ikpt = 1 we add a point below eigen(1) which doesnt
 863 !    contribute to the DOS in any tetrahedron
 864 !scalars
 865  integer :: ieps1,ieps2,itetra
 866  integer :: nn1_1,nn1_2,nn1_3,nn1_4
 867  integer :: nn2_1,nn2_2,nn2_3
 868  integer :: ind_a(3), ind_b(3), ind_c(3)
 869  real(dp_)  :: deltaene1,eps1
 870  real(dp_)  :: deltaene2,eps2
 871 ! real(dp_)  :: gau_prefactor,gau_width,gau_width2
 872  real(dp_)  :: epsilon1(4,4)
 873  real(dp_)  :: epsilon2(4,4)
 874  real(dp_)  :: inv_epsilon1(4,4)
 875  real(dp_)  :: aa(3),bb(3),cc(3)
 876  real(dp_)  :: delaa(3),delbb(3),delcc(3)
 877  real(dp_)  :: delaa0,delbb0,delcc0
 878  real(dp_)  :: inv_delaa(3),inv_delbb(3),inv_delcc(3)
 879  real(dp_)  :: deleps1, deleps2
 880  real(dp_)  :: inv_deleps1
 881  real(dp_)  :: dccde1, dccde1_pre
 882  real(dp_)  :: volconst,volconst_mult
 883  real(dp_)  :: ii0, ii1, ii3
 884 !arrays
 885  integer :: ind_k(4)
 886  real(dp_), allocatable :: tweight_tmp(:,:,:)
 887  real(dp_), allocatable :: dtweightde_tmp(:,:,:)
 888  real(dp_)  :: eigen1_1tetra(4)
 889  real(dp_)  :: eigen2_1tetra(4)
 890 
 891 ! *********************************************************************
 892 
 893  ierr = 0
 894  if (nene1 <= 1 .or. nene2 <= 1)  then
 895    !'get_dbl_tetra_weight: nene must be at least 2'
 896    ierr = 1
 897    return
 898  else
 899    deltaene1 = (enemax1-enemin1) / (nene1-1)
 900    deltaene2 = (enemax2-enemin2) / (nene2-1)
 901  end if
 902 
 903  TETRA_ALLOCATE(tweight_tmp, (4, nene2, nene1))
 904  TETRA_ALLOCATE(dtweightde_tmp, (4, nene2, nene1))
 905 
 906 !print *, "warning: for the moment, heaviside weights are 0. The delta function / DOS weights are the only ones calculated "
 907 
 908  volconst = tetra%vv/4.d0
 909 
 910  ! for each tetrahedron
 911  do itetra=1,tetra%ntetra
 912    ! these are for 1 tetrahedron only.
 913    tweight_tmp = zero
 914    dtweightde_tmp = zero
 915 
 916    volconst_mult = max_occ*volconst*float(tetra%tetra_mult(itetra))
 917 
 918    ! Here we need the original ordering to reference the correct irred kpoints
 919    ! ind_k refers to the index in the full k list of the summits of the present tetrahedra
 920    ! we can forget the order of the summits within the tetrahedron, because eigen1 fixes that
 921    ! order with its increasing value
 922    ind_k(1) = tetra%tetra_full(1,1,itetra)
 923    ind_k(2) = tetra%tetra_full(2,1,itetra)
 924    ind_k(3) = tetra%tetra_full(3,1,itetra)
 925    ind_k(4) = tetra%tetra_full(4,1,itetra)
 926    eigen1_1tetra(1) = eigen1_in(ind_k(1))
 927    eigen1_1tetra(2) = eigen1_in(ind_k(2))
 928    eigen1_1tetra(3) = eigen1_in(ind_k(3))
 929    eigen1_1tetra(4) = eigen1_in(ind_k(4))
 930    call sort_tetra(4,eigen1_1tetra,ind_k,tol14)
 931 
 932    ! re-sort eigen2 values according to order chosen for eigen1. Eigen2 are _not_ in increasing order!
 933    eigen2_1tetra(1) = eigen2_in(ind_k(1))
 934    eigen2_1tetra(2) = eigen2_in(ind_k(2))
 935    eigen2_1tetra(3) = eigen2_in(ind_k(3))
 936    eigen2_1tetra(4) = eigen2_in(ind_k(4))
 937 
 938    ! the epsilons are energy differences for the two eigenvalue sets
 939    epsilon1 = zero
 940    epsilon2 = zero
 941    do ieps1 = 1, 4
 942      do ieps2 = ieps1+1, 4
 943        epsilon1(ieps1,ieps2) = eigen1_1tetra(ieps1)-eigen1_1tetra(ieps2)
 944        epsilon1(ieps2,ieps1) = -epsilon1(ieps1,ieps2)
 945        epsilon2(ieps1,ieps2) = eigen2_1tetra(ieps1)-eigen2_1tetra(ieps2)
 946        epsilon2(ieps2,ieps1) = -epsilon2(ieps1,ieps2)
 947      end do
 948    end do
 949 
 950    ! we precalculate the inverses to avoid doing tons of divisions in the energy loops below
 951    ! Allen formulae only require the inverses of the differences of eigen1 + the a b c below
 952    inv_epsilon1 = zero
 953    do ieps1 = 1, 4
 954      do ieps2 = ieps1+1, 4
 955        if (abs(epsilon1(ieps1,ieps2)) > tol6) then
 956          inv_epsilon1(ieps1,ieps2) = 1.d0 / epsilon1(ieps1,ieps2)
 957          inv_epsilon1(ieps2,ieps1) = -inv_epsilon1(ieps1,ieps2)
 958        end if
 959      end do
 960    end do
 961 
 962    ! these bounds determine the intervals for Omega in Allen paper, and cases A, B, C
 963    nn1_1 = int((eigen1_1tetra(1)-enemin1)/deltaene1)+1
 964    nn1_2 = int((eigen1_1tetra(2)-enemin1)/deltaene1)+1
 965    nn1_3 = int((eigen1_1tetra(3)-enemin1)/deltaene1)+1
 966    nn1_4 = int((eigen1_1tetra(4)-enemin1)/deltaene1)+1
 967 
 968    nn1_1 = max(1,nn1_1)
 969    nn1_1 = min(nn1_1,nene1)
 970    nn1_2 = max(1,nn1_2)
 971    nn1_2 = min(nn1_2,nene1)
 972    nn1_3 = max(1,nn1_3)
 973    nn1_3 = min(nn1_3,nene1)
 974    nn1_4 = max(1,nn1_4)
 975    nn1_4 = min(nn1_4,nene1)
 976 
 977    ! calculate Allen a_i b_i and c_i parameters
 978    ! sort the a_i b_i c_i
 979    !
 980    ! NOTE: indices here go from 1 to 4 instead of 0 to 3 as in Allen...
 981    aa(1) = epsilon2(2,1) * inv_epsilon1(2,1)
 982    aa(2) = epsilon2(3,1) * inv_epsilon1(3,1)
 983    aa(3) = epsilon2(4,1) * inv_epsilon1(4,1)
 984    ind_a = (/2,3,4/)
 985    call sort_tetra(3,aa,ind_a,tol14)
 986    ! aa are now in order a_s a_m a_l !!! Preserve the hash function ind_a to order the positions of k below
 987    delaa(1) = aa(2)-aa(1)
 988    delaa(2) = aa(3)-aa(1)
 989    delaa(3) = aa(3)-aa(2)
 990    inv_delaa = zero
 991    if(delaa(1)> tol6) inv_delaa(1)= 1.0d0 / delaa(1)
 992    if(delaa(2)> tol6) inv_delaa(2)= 1.0d0 / delaa(2)
 993    if(delaa(3)> tol6) inv_delaa(3)= 1.0d0 / delaa(3)
 994 
 995    bb(1) = epsilon2(1,2) * inv_epsilon1(1,2)
 996    bb(2) = epsilon2(3,2) * inv_epsilon1(3,2)
 997    bb(3) = epsilon2(4,2) * inv_epsilon1(4,2)
 998    ind_b = (/1,3,4/)
 999    call sort_tetra(3,bb,ind_b,tol14)
1000    delbb(1) = bb(2)-bb(1)
1001    delbb(2) = bb(3)-bb(1)
1002    delbb(3) = bb(3)-bb(2)
1003    inv_delbb = zero
1004    if(delbb(1)> tol6) inv_delbb(1)= 1.0d0 / delbb(1)
1005    if(delbb(2)> tol6) inv_delbb(2)= 1.0d0 / delbb(2)
1006    if(delbb(3)> tol6) inv_delbb(3)= 1.0d0 / delbb(3)
1007 
1008    cc(1) = epsilon2(1,4) * inv_epsilon1(1,4)
1009    cc(2) = epsilon2(2,4) * inv_epsilon1(2,4)
1010    cc(3) = epsilon2(3,4) * inv_epsilon1(3,4)
1011    ind_c = (/1,2,3/)
1012    call sort_tetra(3,cc,ind_c,tol14)
1013    delcc(1) = cc(2)-cc(1)
1014    delcc(2) = cc(3)-cc(1)
1015    delcc(3) = cc(3)-cc(2)
1016    inv_delcc = zero
1017    if(delcc(1)> tol6) inv_delcc(1)= 1.0d0 / delcc(1)
1018    if(delcc(2)> tol6) inv_delcc(2)= 1.0d0 / delcc(2)
1019    if(delcc(3)> tol6) inv_delcc(3)= 1.0d0 / delcc(3)
1020 
1021    !----------------------------------------------------------------------
1022    ! start main loop A B C over eps1
1023    !----------------------------------------------------------------------
1024 
1025    !
1026    !  interval enemin1 < eps1 < e1 nothing to do
1027    !
1028    !
1029    !  interval e1 < eps1 < e3   CASE A in Allen + first term in B
1030    !
1031    ! NB: eps1 is not updated inside the loop, only between the loops
1032    eps1 = enemin1+nn1_1*deltaene1
1033    deleps1 = eps1-eigen1_1tetra(1) ! this is Omega - omega_0
1034    dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(2,1)*inv_epsilon1(3,1)*inv_epsilon1(4,1)
1035 
1036    ! note we go to nn1_3
1037    do ieps1=nn1_1+1,nn1_3
1038 
1039      dccde1 = dccde1_pre * deleps1  ! this is f_0(Omega)*6*v
1040 
1041      ! at fixed ieps1 we can find the pivot indices for the ieps2 loop
1042      nn2_1 = int((eigen2_1tetra(1)+deleps1*aa(1) -enemin2)/deltaene2)+1
1043      nn2_2 = int((eigen2_1tetra(1)+deleps1*aa(2) -enemin2)/deltaene2)+1
1044      nn2_3 = int((eigen2_1tetra(1)+deleps1*aa(3) -enemin2)/deltaene2)+1
1045 
1046      nn2_1 = max(1,nn2_1)
1047      nn2_1 = min(nn2_1,nene2)
1048      nn2_2 = max(1,nn2_2)
1049      nn2_2 = min(nn2_2,nene2)
1050      nn2_3 = max(1,nn2_3)
1051      nn2_3 = min(nn2_3,nene2)
1052 
1053      inv_deleps1 = 1.0d0 / deleps1
1054 
1055      eps2 = enemin2+nn2_1*deltaene2 ! this is E
1056      deleps2 = eps2 - eigen2_1tetra(1) ! this is E-epsilon_0
1057 
1058      !-----------------------------------------------------------------------
1059      ! This is case AI
1060      !-----------------------------------------------------------------------
1061      do ieps2 = nn2_1+1, nn2_2
1062        ! calculate running value of del "a"  = a-a_s: first term should really mix eps1 and eps2
1063        delaa0 = deleps2*inv_deleps1 - aa(1) ! a - a_s
1064 
1065        ii0 = dccde1*delaa0*inv_delaa(1)*inv_delaa(2) ! this is I_0(Omega E)
1066 
1067        dtweightde_tmp(1,ieps2,ieps1) = dtweightde_tmp(1,ieps2,ieps1) + &
1068 &         ii0*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_a(1),1)* &
1069 &                 (-2.0d0 + delaa0*inv_delaa(1)*epsilon1(ind_a(2),ind_a(1))*inv_epsilon1(ind_a(2),1) &
1070 &                         + delaa0*inv_delaa(2)*epsilon1(ind_a(3),ind_a(1))*inv_epsilon1(ind_a(3),1)))
1071        dtweightde_tmp(ind_a(1),ieps2,ieps1) = dtweightde_tmp(ind_a(1),ieps2,ieps1) + &
1072 &         ii0*0.5d0*deleps1*inv_epsilon1(ind_a(1),1)*(2.0d0 - delaa0*inv_delaa(1) - delaa0*inv_delaa(2))
1073        dtweightde_tmp(ind_a(2),ieps2,ieps1) = dtweightde_tmp(ind_a(2),ieps2,ieps1) + &
1074 &         ii0*0.5d0*delaa0*inv_delaa(1)*deleps1*inv_epsilon1(ind_a(2),1)
1075        dtweightde_tmp(ind_a(3),ieps2,ieps1) = dtweightde_tmp(ind_a(3),ieps2,ieps1) + &
1076 &         ii0*0.5d0*delaa0*inv_delaa(2)*deleps1*inv_epsilon1(ind_a(3),1)
1077        deleps2 = deleps2 + deltaene2
1078      end do
1079 
1080 
1081      eps2 = enemin2+nn2_2*deltaene2 ! this is E
1082      deleps2 = eps2 - eigen2_1tetra(1)  ! E-E_0
1083 
1084      !-----------------------------------------------------------------------
1085      ! This is case AII
1086      !-----------------------------------------------------------------------
1087      do ieps2 = nn2_2+1, nn2_3
1088        ! calculate running value of del "a"  = a_l-a: first term should really mix eps1 and eps2
1089        delaa0 = aa(3) - deleps2*inv_deleps1 ! a_l - a
1090 
1091        ii0 = dccde1*delaa0*inv_delaa(3)*inv_delaa(2) ! this is I_0(Omega E)
1092 
1093        dtweightde_tmp(1,ieps2,ieps1) = dtweightde_tmp(1,ieps2,ieps1) + &
1094 &         ii0*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_a(3),1)* &
1095 &                 (-2.0d0 + delaa0*inv_delaa(3)*epsilon1(ind_a(2),ind_a(3))*inv_epsilon1(ind_a(2),1) &
1096 &                         + delaa0*inv_delaa(2)*epsilon1(ind_a(1),ind_a(3))*inv_epsilon1(ind_a(1),1)))
1097        dtweightde_tmp(ind_a(3),ieps2,ieps1) = dtweightde_tmp(ind_a(3),ieps2,ieps1) + &
1098 &         ii0*0.5d0*deleps1*inv_epsilon1(ind_a(3),1)*(2.0d0 - delaa0*inv_delaa(3) - delaa0*inv_delaa(2))
1099        dtweightde_tmp(ind_a(2),ieps2,ieps1) = dtweightde_tmp(ind_a(2),ieps2,ieps1) + &
1100 &         ii0*0.5d0*delaa0*inv_delaa(3)*deleps1*inv_epsilon1(ind_a(2),1)
1101        dtweightde_tmp(ind_a(1),ieps2,ieps1) = dtweightde_tmp(ind_a(1),ieps2,ieps1) + &
1102 &         ii0*0.5d0*delaa0*inv_delaa(2)*deleps1*inv_epsilon1(ind_a(1),1)
1103 
1104        deleps2 = deleps2 + deltaene2
1105      end do
1106      deleps1 = deleps1 + deltaene1
1107    end do
1108    !
1109    !  interval e2 < eps < e3
1110    !
1111    eps1 = eps1 + (nn1_2-nn1_1)*deltaene1
1112 
1113    deleps1 = eps1-eigen1_1tetra(2) ! Omega - omega_1
1114 
1115    dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(2,1)*inv_epsilon1(3,2)*inv_epsilon1(4,2) ! f1 function
1116    do ieps1=nn1_2+1,nn1_3
1117 
1118      dccde1 = dccde1_pre * deleps1 ! f2(Omega) * 6 * v
1119 
1120      ! at fixed ieps1 we can find the pivot indices for the ieps2 loop
1121      nn2_1 = int((eigen2_1tetra(2)+deleps1*bb(1) -enemin2)/deltaene2)+1
1122      nn2_2 = int((eigen2_1tetra(2)+deleps1*bb(2) -enemin2)/deltaene2)+1
1123      nn2_3 = int((eigen2_1tetra(2)+deleps1*bb(3) -enemin2)/deltaene2)+1
1124 
1125      nn2_1 = max(1,nn2_1)
1126      nn2_1 = min(nn2_1,nene2)
1127      nn2_2 = max(1,nn2_2)
1128      nn2_2 = min(nn2_2,nene2)
1129      nn2_3 = max(1,nn2_3)
1130      nn2_3 = min(nn2_3,nene2)
1131 
1132      inv_deleps1 = 1.0d0 / deleps1
1133 
1134      eps2 = enemin2+nn2_1*deltaene2 ! starting value for E
1135      deleps2 = eps2 - eigen2_1tetra(2) ! E - epsilon_1
1136 
1137      !-----------------------------------------------------------------------
1138      ! This is case BI
1139      !-----------------------------------------------------------------------
1140      do ieps2 = nn2_1+1, nn2_2
1141        ! calculate running value of del "b"  = b-b_s: first term should really mix eps1 and eps2
1142        delbb0 = deleps2*inv_deleps1 - bb(1)
1143 
1144        ii1 = dccde1*delbb0*inv_delbb(1)*inv_delbb(2) ! this is I_1(Omega E)
1145 
1146        ! note negative sign here - we are correcting the I0 a0 term already calculated above
1147        dtweightde_tmp(2,ieps2,ieps1) = dtweightde_tmp(2,ieps2,ieps1) - &
1148 &         ii1*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_b(1),2)* &
1149 &                 (-2.0d0 + delbb0*inv_delbb(1)*epsilon1(ind_b(2),ind_b(1))*inv_epsilon1(ind_b(2),2) &
1150 &                         + delbb0*inv_delbb(2)*epsilon1(ind_b(3),ind_b(1))*inv_epsilon1(ind_b(3),2)))
1151        dtweightde_tmp(ind_b(1),ieps2,ieps1) = dtweightde_tmp(ind_b(1),ieps2,ieps1) - &
1152 &         ii1*0.5d0*deleps1*inv_epsilon1(ind_b(1),2)*(2.0d0 - delbb0*inv_delbb(1) - delbb0*inv_delbb(2))
1153        dtweightde_tmp(ind_b(2),ieps2,ieps1) = dtweightde_tmp(ind_b(2),ieps2,ieps1) - &
1154 &         ii1*0.5d0*delbb0*inv_delbb(1)*deleps1*inv_epsilon1(ind_b(2),2)
1155        dtweightde_tmp(ind_b(3),ieps2,ieps1) = dtweightde_tmp(ind_b(3),ieps2,ieps1) - &
1156 &         ii1*0.5d0*delbb0*inv_delbb(2)*deleps1*inv_epsilon1(ind_b(3),2)
1157        deleps2 = deleps2 + deltaene2
1158      end do
1159 
1160      eps2 = enemin2+nn2_2*deltaene2
1161      deleps2 = eps2 - eigen2_1tetra(2)
1162 
1163      !-----------------------------------------------------------------------
1164      ! This is case BII
1165      !-----------------------------------------------------------------------
1166      do ieps2 = nn2_2+1, nn2_3
1167        ! calculate running value of del "b"  = b_l-b: first term should really mix eps1 and eps2
1168        delbb0 = bb(3) - deleps2*inv_deleps1
1169 
1170        ii1 = dccde1*delbb0*inv_delbb(3)*inv_delbb(2) ! this is I_1(Omega E)
1171 
1172        ! note negative sign here - we are correcting the I0 a0 term already calculated above
1173        dtweightde_tmp(2,ieps2,ieps1) = dtweightde_tmp(2,ieps2,ieps1) - &
1174 &         ii1*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_b(3),2)* &
1175 &                 (-2.0d0 + delbb0*inv_delbb(3)*epsilon1(ind_b(2),ind_b(3))*inv_epsilon1(ind_b(2),2) &
1176 &                         + delbb0*inv_delbb(2)*epsilon1(ind_b(1),ind_b(3))*inv_epsilon1(ind_b(1),2)))
1177        dtweightde_tmp(ind_b(3),ieps2,ieps1) = dtweightde_tmp(ind_b(3),ieps2,ieps1) - &
1178 &         ii1*0.5d0*deleps1*inv_epsilon1(ind_b(3),2)*(2.0d0 - delbb0*inv_delbb(3) - delbb0*inv_delbb(2))
1179        dtweightde_tmp(ind_b(2),ieps2,ieps1) = dtweightde_tmp(ind_b(2),ieps2,ieps1) - &
1180 &         ii1*0.5d0*delbb0*inv_delbb(3)*deleps1*inv_epsilon1(ind_b(2),2)
1181        dtweightde_tmp(ind_b(1),ieps2,ieps1) = dtweightde_tmp(ind_b(1),ieps2,ieps1) - &
1182 &         ii1*0.5d0*delbb0*inv_delbb(2)*deleps1*inv_epsilon1(ind_b(1),2)
1183 
1184        deleps2 = deleps2 + deltaene2
1185      end do
1186 
1187      deleps1 = deleps1 + deltaene1
1188    end do
1189 
1190    !
1191    !  interval e3 < eps < e4
1192    !
1193    eps1 = eps1 + (nn1_3-nn1_2)*deltaene1
1194    deleps1 = eps1-eigen1_1tetra(4)
1195    dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(4,1)*inv_epsilon1(4,2)*inv_epsilon1(4,3)
1196    do ieps1=nn1_3+1,nn1_4
1197      ! note - sign from definition of f3
1198      dccde1 = -dccde1_pre *  deleps1 ! f3(Omega) * 6 * v
1199 
1200      ! at fixed ieps1 we can find the pivot indices for the ieps2 loop
1201      ! NB: order is inverted for cc because deleps1 is defined negative (Omega is always less than omega_3)
1202      nn2_1 = int((eigen2_1tetra(4)+deleps1*cc(3) -enemin2)/deltaene2)+1
1203      nn2_2 = int((eigen2_1tetra(4)+deleps1*cc(2) -enemin2)/deltaene2)+1
1204      nn2_3 = int((eigen2_1tetra(4)+deleps1*cc(1) -enemin2)/deltaene2)+1
1205 
1206      nn2_1 = max(1,nn2_1)
1207      nn2_1 = min(nn2_1,nene2)
1208      nn2_2 = max(1,nn2_2)
1209      nn2_2 = min(nn2_2,nene2)
1210      nn2_3 = max(1,nn2_3)
1211      nn2_3 = min(nn2_3,nene2)
1212      inv_deleps1 = 1.0d0 / deleps1
1213 
1214      eps2 = enemin2+nn2_1*deltaene2 ! starting value for E
1215      deleps2 = eps2 - eigen2_1tetra(4) ! E - epsilon_3
1216 
1217      !-----------------------------------------------------------------------
1218      ! This is case CII
1219      !-----------------------------------------------------------------------
1220      do ieps2 = nn2_1+1, nn2_2
1221        ! calculate running value of del "c"  = c_l-c: first term should really mix eps1 and eps2
1222        delcc0 = cc(3) - deleps2*inv_deleps1
1223 
1224        ii3 = dccde1*delcc0*inv_delcc(3)*inv_delcc(2) ! this is I_3(Omega E)
1225 
1226        dtweightde_tmp(4,ieps2,ieps1) = dtweightde_tmp(4,ieps2,ieps1) + &
1227 &         ii3*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_c(3),4)* &
1228 &                 (-2.0d0 + delcc0*inv_delcc(3)*epsilon1(ind_c(2),ind_c(3))*inv_epsilon1(ind_c(2),4) &
1229 &                         + delcc0*inv_delcc(2)*epsilon1(ind_c(1),ind_c(3))*inv_epsilon1(ind_c(1),4)))
1230        dtweightde_tmp(ind_c(3),ieps2,ieps1) = dtweightde_tmp(ind_c(3),ieps2,ieps1) + &
1231 &         ii3*0.5d0*deleps1*inv_epsilon1(ind_c(3),4)*(2.0d0 - delcc0*inv_delcc(3) - delcc0*inv_delcc(2))
1232        dtweightde_tmp(ind_c(2),ieps2,ieps1) = dtweightde_tmp(ind_c(2),ieps2,ieps1) + &
1233 &         ii3*0.5d0*delcc0*inv_delcc(3)*deleps1*inv_epsilon1(ind_c(2),4)
1234        dtweightde_tmp(ind_c(1),ieps2,ieps1) = dtweightde_tmp(ind_c(1),ieps2,ieps1) + &
1235 &         ii3*0.5d0*delcc0*inv_delcc(2)*deleps1*inv_epsilon1(ind_c(1),4)
1236 
1237        deleps2 = deleps2 + deltaene2
1238      end do
1239 
1240 
1241      eps2 = enemin2+nn2_2*deltaene2
1242      deleps2 = eps2 - eigen2_1tetra(4)
1243 
1244      !-----------------------------------------------------------------------
1245      ! This is case CI
1246      !-----------------------------------------------------------------------
1247      do ieps2 = nn2_2+1, nn2_3
1248        ! calculate running value of del "c"  = c-c_s: first term should really mix eps1 and eps2
1249        delcc0 = deleps2*inv_deleps1 - cc(1) ! c - c_s
1250 
1251        ii3 = dccde1*delcc0*inv_delcc(1)*inv_delcc(2) ! this is I_3(Omega E)
1252 
1253        dtweightde_tmp(4,ieps2,ieps1) = dtweightde_tmp(4,ieps2,ieps1) + &
1254 &         ii3*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_c(1),4)* &
1255 &                 (-2.0d0 + delcc0*inv_delcc(1)*epsilon1(ind_c(2),ind_c(1))*inv_epsilon1(ind_c(2),4) &
1256 &                         + delcc0*inv_delcc(2)*epsilon1(ind_c(3),ind_c(1))*inv_epsilon1(ind_c(3),4)))
1257        dtweightde_tmp(ind_c(1),ieps2,ieps1) = dtweightde_tmp(ind_c(1),ieps2,ieps1) + &
1258 &         ii3*0.5d0*deleps1*inv_epsilon1(ind_c(1),4)*(2.0d0 - delcc0*inv_delcc(1) - delcc0*inv_delcc(2))
1259        dtweightde_tmp(ind_c(2),ieps2,ieps1) = dtweightde_tmp(ind_c(2),ieps2,ieps1) + &
1260 &         ii3*0.5d0*delcc0*inv_delcc(1)*deleps1*inv_epsilon1(ind_c(2),4)
1261        dtweightde_tmp(ind_c(3),ieps2,ieps1) = dtweightde_tmp(ind_c(3),ieps2,ieps1) + &
1262 &         ii3*0.5d0*delcc0*inv_delcc(2)*deleps1*inv_epsilon1(ind_c(3),4)
1263        deleps2 = deleps2 + deltaene2
1264      end do
1265 
1266      deleps1 = deleps1 + deltaene1
1267    end do
1268 
1269    eps1 = eps1 + (nn1_4-nn1_3)*deltaene1
1270    !
1271    !
1272    !  interval e4 < eps < enemax
1273    !
1274    do ieps1=nn1_4+1,nene1
1275      ! dtweightde unchanged by this tetrahedron
1276    end do
1277 
1278    !
1279    !  if we have a fully degenerate tetrahedron,
1280    !  1) the tweight is a Heaviside (step) function, which is correct above, but
1281    !  2) the dtweightde should contain a Dirac function: add a Gaussian here
1282 
1283    ! TODO : add treatment in double tetra case
1284    !  end degenerate tetrahedron if
1285 
1286    ! NOTE: the following blas calls are not working systematically, or do not give speed ups, strange...
1287    !   call daxpy (nene, 1.d0, dtweightde_tmp(:,1), 1, dtweightde_t(:,ind_ibz(1)), 1)
1288    !   call daxpy (nene, 1.d0, dtweightde_tmp(:,2), 1, dtweightde_t(:,ind_ibz(2)), 1)
1289    !   call daxpy (nene, 1.d0, dtweightde_tmp(:,3), 1, dtweightde_t(:,ind_ibz(3)), 1)
1290    !   call daxpy (nene, 1.d0, dtweightde_tmp(:,4), 1, dtweightde_t(:,ind_ibz(4)), 1)
1291 
1292    do ieps2 = 1, nene2
1293      dtweightde(ind_k(1),:,ieps2) = dtweightde(ind_k(1),:,ieps2) + dtweightde_tmp(1,ieps2,:)
1294      dtweightde(ind_k(2),:,ieps2) = dtweightde(ind_k(2),:,ieps2) + dtweightde_tmp(2,ieps2,:)
1295      dtweightde(ind_k(3),:,ieps2) = dtweightde(ind_k(3),:,ieps2) + dtweightde_tmp(3,ieps2,:)
1296      dtweightde(ind_k(4),:,ieps2) = dtweightde(ind_k(4),:,ieps2) + dtweightde_tmp(4,ieps2,:)
1297      !tweight(nkpt,nene1,nene2)
1298    end do
1299 
1300  end do ! itetra
1301 
1302  ! transpose: otherwise the data access is crap and the code slows by an order of magnitude
1303  TETRA_DEALLOCATE(tweight_tmp)
1304  TETRA_DEALLOCATE(dtweightde_tmp)
1305 
1306 end subroutine get_dbl_tetra_weight

m_tetrahedron/get_onetetra_ [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 get_onetetra_

FUNCTION

 Private function to calculate the contributions to the weights due to a single tetrahedron.
 Extracted from get_tetra_weight

PARENTS

CHILDREN

SOURCE

1551 pure subroutine get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,nene,bcorr, &
1552 &  tweight_tmp,dtweightde_tmp)
1553 
1554 
1555 !This section has been created automatically by the script Abilint (TD).
1556 !Do not modify the following lines by hand.
1557 #undef ABI_FUNC
1558 #define ABI_FUNC 'get_onetetra_'
1559 !End of the abilint section
1560 
1561  implicit none
1562 
1563 !Arguments ------------------------------------
1564 !scalars
1565  integer,intent(in) :: nene,bcorr,itetra
1566  type(t_tetrahedron), intent(in) :: tetra
1567  real(dp_) ,intent(in) :: enemax,enemin,max_occ
1568 !arrays
1569  ! MGTODO: This layout is not optimal (lots of cache thrashing, I will optimize it later on)
1570  real(dp_), intent(out) ::  tweight_tmp(nene, 4)
1571  real(dp_), intent(out) :: dtweightde_tmp(nene, 4)
1572  real(dp_),intent(in)  :: eigen_1tetra(4)
1573 
1574 !Local variables-------------------------------
1575 !  needed for gaussian replacement of Dirac functions
1576 !  the three coefficients of the DOS as quadratic form,
1577 !    in the interval [eig(ikpt-1), eig(ikpt)]
1578 !    for ikpt = 1 we add a point below eigen(1) which doesnt
1579 !    contribute to the DOS in any tetrahedron
1580 !scalars
1581  integer :: ieps,nn1,nn2,nn3,nn4
1582  real(dp_)  :: cc,cc1,cc2,cc3,dcc1de,dcc2de,dcc3de,dccde,deltaene,eps
1583  real(dp_)  :: epsilon21,epsilon31,epsilon32,epsilon41,epsilon42,epsilon43
1584  real(dp_)  :: gau_prefactor,gau_width,gau_width2,inv_epsilon21,inv_epsilon31,gval
1585  real(dp_)  :: inv_epsilon32,inv_epsilon41,inv_epsilon42,inv_epsilon43
1586  real(dp_)  :: deleps1, deleps2, deleps3, deleps4
1587  real(dp_)  :: invepsum, cc_pre, dccde_pre
1588  real(dp_)  :: cc1_pre, cc2_pre, cc3_pre
1589  real(dp_)  :: cc_tmp, dccde_tmp
1590  real(dp_)  :: dcc1de_pre, dcc2de_pre, dcc3de_pre
1591  real(dp_)  :: tmp,volconst,volconst_mult
1592 
1593 ! *********************************************************************
1594 
1595  volconst = tetra%vv/4.d0
1596 
1597  deltaene = (enemax-enemin) / (nene-1)
1598 
1599  ! This is output
1600  tweight_tmp = zero; dtweightde_tmp = zero
1601 
1602  volconst_mult = max_occ*volconst*float(tetra%tetra_mult(itetra))
1603 
1604  ! all notations are from Blochl PRB 49 16223 [[cite:Bloechl1994a]] Appendix B
1605  epsilon21 = eigen_1tetra(2)-eigen_1tetra(1)
1606  epsilon31 = eigen_1tetra(3)-eigen_1tetra(1)
1607  epsilon41 = eigen_1tetra(4)-eigen_1tetra(1)
1608  epsilon32 = eigen_1tetra(3)-eigen_1tetra(2)
1609  epsilon42 = eigen_1tetra(4)-eigen_1tetra(2)
1610  epsilon43 = eigen_1tetra(4)-eigen_1tetra(3)
1611  inv_epsilon21 = zero
1612  inv_epsilon31 = zero
1613  inv_epsilon41 = zero
1614  inv_epsilon32 = zero
1615  inv_epsilon42 = zero
1616  inv_epsilon43 = zero
1617  if (epsilon21 > tol6) inv_epsilon21 = 1.d0 / epsilon21
1618  if (epsilon31 > tol6) inv_epsilon31 = 1.d0 / epsilon31
1619  if (epsilon41 > tol6) inv_epsilon41 = 1.d0 / epsilon41
1620  if (epsilon32 > tol6) inv_epsilon32 = 1.d0 / epsilon32
1621  if (epsilon42 > tol6) inv_epsilon42 = 1.d0 / epsilon42
1622  if (epsilon43 > tol6) inv_epsilon43 = 1.d0 / epsilon43
1623  nn1 = int((eigen_1tetra(1)-enemin)/deltaene)+1
1624  nn2 = int((eigen_1tetra(2)-enemin)/deltaene)+1
1625  nn3 = int((eigen_1tetra(3)-enemin)/deltaene)+1
1626  nn4 = int((eigen_1tetra(4)-enemin)/deltaene)+1
1627 
1628  nn1 = max(1,nn1)
1629  nn1 = min(nn1,nene)
1630  nn2 = max(1,nn2)
1631  nn2 = min(nn2,nene)
1632  nn3 = max(1,nn3)
1633  nn3 = min(nn3,nene)
1634  nn4 = max(1,nn4)
1635  nn4 = min(nn4,nene)
1636 
1637  eps = enemin+nn1*deltaene
1638  !
1639  !interval enemin < eps < e1 nothing to do
1640  !
1641  !
1642  !interval e1 < eps < e2
1643  !
1644  deleps1 = eps-eigen_1tetra(1)
1645  cc_pre = volconst_mult*inv_epsilon21*inv_epsilon31*inv_epsilon41
1646  invepsum = inv_epsilon21+inv_epsilon31+inv_epsilon41
1647  dccde_pre = 3.d0*volconst_mult*inv_epsilon21*inv_epsilon31*inv_epsilon41
1648  do ieps=nn1+1,nn2
1649    cc = cc_pre * deleps1*deleps1*deleps1
1650    tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + cc*(4.d0-deleps1*invepsum)
1651    tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + cc*deleps1*inv_epsilon21
1652    tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + cc*deleps1*inv_epsilon31
1653    tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + cc*deleps1*inv_epsilon41
1654 
1655    dccde = dccde_pre * deleps1*deleps1
1656    dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + dccde*(4.d0 - deleps1*invepsum) -cc*invepsum
1657    dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + (dccde*deleps1 + cc) * inv_epsilon21
1658    dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + (dccde*deleps1 + cc) * inv_epsilon31
1659    dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + (dccde*deleps1 + cc) * inv_epsilon41
1660 
1661    if (bcorr == 1) then
1662      ! bxu, correction terms based on Bloechl's paper
1663      tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1664 &     4.d0*dccde_pre*deleps1*deleps1*(epsilon21+epsilon31+epsilon41)/40.d0
1665      tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1666 &     4.d0*dccde_pre*deleps1*deleps1*(-epsilon21+epsilon32+epsilon42)/40.d0
1667      tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1668 &     4.d0*dccde_pre*deleps1*deleps1*(-epsilon31-epsilon32+epsilon43)/40.d0
1669      tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1670 &     4.d0*dccde_pre*deleps1*deleps1*(-epsilon41-epsilon42-epsilon43)/40.d0
1671 
1672      dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + &
1673 &     8.d0*dccde_pre*deleps1*(epsilon21+epsilon31+epsilon41)/40.d0
1674      dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + &
1675 &     8.d0*dccde_pre*deleps1*(-epsilon21+epsilon32+epsilon42)/40.d0
1676      dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + &
1677 &     8.d0*dccde_pre*deleps1*(-epsilon31-epsilon32+epsilon43)/40.d0
1678      dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + &
1679 &     8.d0*dccde_pre*deleps1*(-epsilon41-epsilon42-epsilon43)/40.d0
1680    end if
1681 
1682    deleps1 = deleps1 + deltaene
1683  end do
1684 
1685  eps = eps + (nn2-nn1)*deltaene
1686  !
1687  !  interval e2 < eps < e3
1688  !
1689  deleps1 = eps-eigen_1tetra(1)
1690  deleps2 = eps-eigen_1tetra(2)
1691  deleps3 = eigen_1tetra(3)-eps
1692  deleps4 = eigen_1tetra(4)-eps
1693 
1694  cc1_pre = volconst_mult*inv_epsilon31*inv_epsilon41
1695  cc2_pre = volconst_mult*inv_epsilon41*inv_epsilon32*inv_epsilon31
1696  cc3_pre = volconst_mult*inv_epsilon42*inv_epsilon32*inv_epsilon41
1697 
1698  dcc1de_pre = 2.d0*cc1_pre
1699  dcc2de_pre = cc2_pre
1700  dcc3de_pre = cc3_pre
1701  do ieps=nn2+1,nn3
1702    cc1 = cc1_pre * deleps1*deleps1
1703    cc2 = cc2_pre * deleps1*deleps2*deleps3
1704    cc3 = cc3_pre * deleps2*deleps2*deleps4
1705 
1706    tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1707 &   cc1 + (cc1+cc2)*deleps3*inv_epsilon31 + (cc1+cc2+cc3)*deleps4*inv_epsilon41
1708    tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1709 &   cc1+cc2+cc3+(cc2+cc3)*deleps3*inv_epsilon32 + cc3*deleps4*inv_epsilon42
1710    tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1711 &   (cc1+cc2)*deleps1*inv_epsilon31 + (cc2+cc3)*deleps2*inv_epsilon32
1712    tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1713 &   (cc1+cc2+cc3)*deleps1*inv_epsilon41 + cc3*deleps2*inv_epsilon42
1714 
1715 
1716    dcc1de = dcc1de_pre * deleps1
1717    dcc2de = dcc2de_pre * (-deleps1*deleps2  +deleps1*deleps3  +deleps2*deleps3)
1718    dcc3de = dcc3de_pre * (2.d0*deleps2*deleps4  -deleps2*deleps2)
1719 
1720    dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) &
1721 &   + dcc1de &
1722 &   + ((dcc1de+dcc2de)*deleps3 -(cc1+cc2)) * inv_epsilon31 &
1723 &   + ((dcc1de+dcc2de+dcc3de)*deleps4 -(cc1+cc2+cc3)) * inv_epsilon41
1724    dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) &
1725 &   + dcc1de+dcc2de+dcc3de &
1726 &   + ((dcc2de+dcc3de)*deleps3 -(cc2+cc3) ) * inv_epsilon32 &
1727 &   + (dcc3de*deleps4  -cc3 ) * inv_epsilon42
1728    dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) &
1729 &   + ((dcc1de+dcc2de)*deleps1 + (cc1+cc2) ) * inv_epsilon31 &
1730 &   + ((dcc2de+dcc3de)*deleps2 + (cc2+cc3) ) * inv_epsilon32
1731    dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) &
1732 &   + ((dcc1de+dcc2de+dcc3de)*deleps1 + (cc1+cc2+cc3) ) * inv_epsilon41 &
1733 &   + (dcc3de*deleps2 + cc3) * inv_epsilon42
1734 
1735  if (bcorr == 1) then
1736    ! bxu, correction terms based on Bloechl's paper
1737    ! The correction terms may cause the dtweightde become negative
1738    tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1739 &   4.d0*cc1_pre* &
1740 &   (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1741 &   (epsilon21+epsilon31+epsilon41)/40.d0
1742    tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1743 &   4.d0*cc1_pre* &
1744 &   (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1745 &   (-epsilon21+epsilon32+epsilon42)/40.d0
1746    tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1747 &   4.d0*cc1_pre* &
1748 &   (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1749 &   (-epsilon31-epsilon32+epsilon43)/40.d0
1750    tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1751 &   4.d0*cc1_pre* &
1752 &   (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* &
1753 &   (-epsilon41-epsilon42-epsilon43)/40.d0
1754 
1755    dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + &
1756 &   4.d0*cc1_pre* &
1757 &   (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1758 &   (epsilon21+epsilon31+epsilon41)/40.d0
1759    dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + &
1760 &   4.d0*cc1_pre* &
1761 &   (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1762 &   (-epsilon21+epsilon32+epsilon42)/40.d0
1763    dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + &
1764 &   4.d0*cc1_pre* &
1765 &   (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1766 &   (-epsilon31-epsilon32+epsilon43)/40.d0
1767    dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + &
1768 &   4.d0*cc1_pre* &
1769 &   (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* &
1770 &   (-epsilon41-epsilon42-epsilon43)/40.d0
1771   end if
1772 
1773   deleps1 = deleps1 + deltaene
1774   deleps2 = deleps2 + deltaene
1775   deleps3 = deleps3 - deltaene
1776   deleps4 = deleps4 - deltaene
1777  end do
1778 
1779  eps = eps + (nn3-nn2)*deltaene
1780  !
1781  !  interval e3 < eps < e4
1782  !
1783  deleps4 = eigen_1tetra(4)-eps
1784  cc_pre = volconst_mult*inv_epsilon41*inv_epsilon42*inv_epsilon43
1785  invepsum = inv_epsilon41+inv_epsilon42+inv_epsilon43
1786  dccde_pre = -3.d0*cc_pre
1787  do ieps=nn3+1,nn4
1788    cc = cc_pre * deleps4*deleps4*deleps4
1789    cc_tmp = cc * deleps4
1790    tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + volconst_mult - cc_tmp*inv_epsilon41
1791    tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + volconst_mult - cc_tmp*inv_epsilon42
1792    tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + volconst_mult - cc_tmp*inv_epsilon43
1793    tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + volconst_mult - cc*4.d0 + cc_tmp*invepsum
1794 
1795 
1796    dccde = dccde_pre * deleps4*deleps4
1797    dccde_tmp = -dccde*deleps4 + cc
1798    dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + dccde_tmp * inv_epsilon41
1799    dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + dccde_tmp * inv_epsilon42
1800    dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + dccde_tmp * inv_epsilon43
1801    dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) - dccde*4.d0 - dccde_tmp*invepsum
1802 
1803    if (bcorr == 1) then
1804      ! bxu, correction terms based on Bloechl's paper
1805      ! The correction terms may cause the dtweightde become negative
1806      tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + &
1807 &     12.d0*cc_pre*deleps4*deleps4*(epsilon21+epsilon31+epsilon41)/40.d0
1808      tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + &
1809 &     12.d0*cc_pre*deleps4*deleps4*(-epsilon21+epsilon32+epsilon42)/40.d0
1810      tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + &
1811 &     12.d0*cc_pre*deleps4*deleps4*(-epsilon31-epsilon32+epsilon43)/40.d0
1812      tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + &
1813 &     12.d0*cc_pre*deleps4*deleps4*(-epsilon41-epsilon42-epsilon43)/40.d0
1814 
1815      dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) - &
1816 &     24.d0*cc_pre*deleps4*(epsilon21+epsilon31+epsilon41)/40.d0
1817      dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) - &
1818 &     24.d0*cc_pre*deleps4*(-epsilon21+epsilon32+epsilon42)/40.d0
1819      dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) - &
1820 &     24.d0*cc_pre*deleps4*(-epsilon31-epsilon32+epsilon43)/40.d0
1821      dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) - &
1822 &     24.d0*cc_pre*deleps4*(-epsilon41-epsilon42-epsilon43)/40.d0
1823    end if
1824 
1825    deleps4 = deleps4 - deltaene
1826  end do
1827  eps = eps + (nn4-nn3)*deltaene
1828  !
1829  !
1830  !  interval e4 < eps < enemax
1831  !
1832  do ieps=nn4+1,nene
1833    tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + volconst_mult
1834    tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + volconst_mult
1835    tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + volconst_mult
1836    tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + volconst_mult
1837    ! dtweightde unchanged by this tetrahedron
1838  end do
1839 
1840  !
1841  !  if we have a fully degenerate tetrahedron,
1842  !  1) the tweight is a Heaviside (step) function, which is correct above, but
1843  !  2) the dtweightde should contain a Dirac function: add a Gaussian here
1844  !
1845  if (epsilon41 < tol6) then
1846 
1847    !  to ensure the gaussian will integrate properly:
1848    !  WARNING: this smearing could be problematic if too large
1849    !  and doesnt integrate well if its too small
1850    gau_width = 10.0d0*deltaene
1851    gau_width2 = 1.0 / gau_width / gau_width
1852    gau_prefactor = volconst_mult / gau_width / sqrtpi
1853    !
1854    ! average position since bracket for epsilon41 is relatively large
1855    cc = (eigen_1tetra(1)+eigen_1tetra(2)+eigen_1tetra(3)+eigen_1tetra(4))/4.d0
1856    eps = enemin
1857    do ieps=1,nene
1858      tmp = eps - cc
1859      gval = gau_prefactor*exp(-tmp*tmp*gau_width2)
1860      !dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + gval
1861      !dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + gval
1862      !dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + gval
1863      dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + gval
1864      eps = eps + deltaene
1865    end do
1866  end if ! end degenerate tetrahedron if
1867 
1868 end subroutine get_onetetra_

m_tetrahedron/get_tetra_weight [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 get_tetra_weight

FUNCTION

 calculate integration weights and their derivatives from Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]

INPUTS

 eigen_in(nkpt)=eigenenergies for each k point
 enemin=minimal energy for DOS
 enemax=maximal energy for DOS
 max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)
 nene=number of energies for DOS
 nkpt=number of irreducible kpoints
 tetra<t_tetrahedron>
   %ntetra=number of tetrahedra
   %tetra_full(4,2,ntetra)=for each irred tetrahedron, the list of k point vertices
     1 -> irred kpoint   2 -> fullkpt
   %tetra_mult(ntetra)=for each irred tetrahedron, its multiplicity
   %vv = ratio of volume of one tetrahedron in reciprocal space to full BZ volume
 bcorr=1 to include Blochl correction else 0.
 comm=MPI communicator

OUTPUT

  tweight(nkpt,nene) = integration weights for each irred kpoint from all adjacent tetrahedra
  dtweightde(nkpt,nene) = derivative of tweight wrt energy

PARENTS

      ep_el_weights,ep_fs_weights,ep_ph_weights,thmeig

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

620 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
621 ! THIS FUNCTION IS DEPRECATED, USE tetra_blochl_weights
622 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
623 subroutine get_tetra_weight(eigen_in,enemin,enemax,max_occ,nene,nkpt,tetra,&
624   bcorr,tweight,dtweightde,comm)
625 
626 
627 !This section has been created automatically by the script Abilint (TD).
628 !Do not modify the following lines by hand.
629 #undef ABI_FUNC
630 #define ABI_FUNC 'get_tetra_weight'
631 !End of the abilint section
632 
633  implicit none
634 
635 !Arguments ------------------------------------
636 !scalars
637  integer,intent(in) :: nene,nkpt,bcorr,comm
638  type(t_tetrahedron), intent(in) :: tetra
639  real(dp_) ,intent(in) :: enemax,enemin,max_occ
640 !arrays
641  real(dp_) ,intent(in) :: eigen_in(nkpt)
642  real(dp_) ,intent(out) :: dtweightde(nkpt,nene),tweight(nkpt,nene)
643 
644 !Local variables-------------------------------
645  real(dp_) , allocatable :: dtweightde_ek(:, :), tweight_ek(:, :)
646 
647 ! *********************************************************************
648 
649  TETRA_ALLOCATE(dtweightde_ek, (nene, nkpt))
650  TETRA_ALLOCATE(tweight_ek, (nene, nkpt))
651 
652  call tetra_blochl_weights(tetra,eigen_in,enemin,enemax,max_occ,nene,nkpt,bcorr,tweight_ek,dtweightde_ek,comm)
653 
654  ! transpose: otherwise the data access is crap and the code slows by an order of magnitude
655  tweight    = transpose(tweight_ek)
656  dtweightde = transpose(dtweightde_ek)
657 
658  TETRA_DEALLOCATE(dtweightde_ek)
659  TETRA_DEALLOCATE(tweight_ek)
660 
661 end subroutine get_tetra_weight

m_tetrahedron/init_tetra [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 init_tetra

FUNCTION

 get tetrahedra characterized by apexes

INPUTS

  indkpt(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

OUTPUT

  tetra%tetra_full(4,2,ntetra)=for each tetrahedron,
     the different instances of the tetrahedron (fullbz kpoints)
  tetra%tetra_mult(ntetra) = store multiplicity of each irred tetrahedron
  tetra%tetra_wrap(3,4,ntetra) = store flag to wrap tetrahedron summit into IBZ
  tetra%ntetra = final number of irred tetra (dimensions of tetra_* remain larger)
  tetra%vv = tetrahedron volume divided by full BZ volume

PARENTS

      ep_el_weights,ep_fs_weights,ep_ph_weights,m_fstab,m_kpts,m_phonons
      thmeig

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

192 subroutine init_tetra (indkpt,gprimd,klatt,kpt_fullbz,nkpt_fullbz,tetra,ierr,errorstring)
193 
194 
195 !This section has been created automatically by the script Abilint (TD).
196 !Do not modify the following lines by hand.
197 #undef ABI_FUNC
198 #define ABI_FUNC 'init_tetra'
199 !End of the abilint section
200 
201  implicit none
202 
203 !Arguments ------------------------------------
204 !scalars
205  integer,intent(in) :: nkpt_fullbz
206 !arrays
207  integer,intent(in) :: indkpt(nkpt_fullbz)
208  real(dp_) ,intent(in) :: gprimd(3,3),klatt(3,3),kpt_fullbz(3,nkpt_fullbz)
209  integer, intent(out) :: ierr
210  character(len=80), intent(out) :: errorstring
211  type(t_tetrahedron),intent(out) :: tetra
212 
213 !Local variables-------------------------------
214 !scalars
215  integer :: ialltetra,ikpt2,ikpt_full,isummit,itetra,jalltetra,jsummit
216  integer :: symrankkpt,mtetra,itmp,ntetra_irred
217  real(dp_) :: shift1,shift2,shift3, rcvol,hashfactor
218  type(kptrank_type) :: kptrank_t
219 !arrays
220  integer :: tetra_shifts(3,4,6)  ! 3 dimensions, 4 summits, and 6 tetrahedra / kpoint box
221  real(dp_)  :: k1(3),k2(3),k3(3)
222  integer,allocatable :: tetra_full_(:,:,:)
223  integer,allocatable :: tetra_mult_(:)
224  integer,allocatable :: tetra_wrap_(:,:,:)
225  integer, allocatable :: reforder(:)
226  integer, allocatable :: irred_itetra(:)
227  real(dp_), allocatable :: tetra_hash(:)
228 
229 ! *********************************************************************
230 
231  ierr = 0
232  errorstring = ""
233 !jmb
234  shift1 = zero
235  shift2 = zero
236  shift3 = zero
237 
238  tetra%klatt = klatt
239 
240  mtetra = 6 * nkpt_fullbz
241  TETRA_ALLOCATE(tetra_full_, (4,2,mtetra))
242  TETRA_ALLOCATE(tetra_mult_, (mtetra))
243  TETRA_ALLOCATE(tetra_wrap_, (3,4,mtetra))
244 
245  tetra_mult_ = 1
246  tetra_full_ = 0
247  tetra_wrap_ = 0
248 
249 ! tetra_shifts(:,1,1) = (/0,0,0/)
250 ! tetra_shifts(:,2,1) = (/0,1,0/)
251 ! tetra_shifts(:,3,1) = (/0,1,1/)
252 ! tetra_shifts(:,4,1) = (/1,1,0/)
253 ! tetra_shifts(:,1,2) = (/0,0,0/)
254 ! tetra_shifts(:,2,2) = (/0,1,1/)
255 ! tetra_shifts(:,3,2) = (/1,1,0/)
256 ! tetra_shifts(:,4,2) = (/1,1,1/)
257 ! tetra_shifts(:,1,3) = (/0,0,0/)
258 ! tetra_shifts(:,2,3) = (/1,0,0/)
259 ! tetra_shifts(:,3,3) = (/1,1,0/)
260 ! tetra_shifts(:,4,3) = (/1,1,1/)
261 ! tetra_shifts(:,1,4) = (/0,0,0/)
262 ! tetra_shifts(:,2,4) = (/0,0,1/)
263 ! tetra_shifts(:,3,4) = (/1,0,0/)
264 ! tetra_shifts(:,4,4) = (/1,1,1/)
265 ! tetra_shifts(:,1,5) = (/0,0,1/)
266 ! tetra_shifts(:,2,5) = (/1,0,0/)
267 ! tetra_shifts(:,3,5) = (/1,0,1/)
268 ! tetra_shifts(:,4,5) = (/1,1,1/)
269 ! tetra_shifts(:,1,6) = (/0,0,0/)
270 ! tetra_shifts(:,2,6) = (/0,0,1/)
271 ! tetra_shifts(:,3,6) = (/0,1,1/)
272 ! tetra_shifts(:,4,6) = (/1,1,1/)
273 
274  ! bxu, the following division scheme is according to Bloechl's paper
275  tetra_shifts(:,1,1) = (/0,0,0/)
276  tetra_shifts(:,2,1) = (/1,0,0/)
277  tetra_shifts(:,3,1) = (/0,1,0/)
278  tetra_shifts(:,4,1) = (/1,0,1/)
279  tetra_shifts(:,1,2) = (/1,0,0/)
280  tetra_shifts(:,2,2) = (/1,1,0/)
281  tetra_shifts(:,3,2) = (/0,1,0/)
282  tetra_shifts(:,4,2) = (/1,0,1/)
283  tetra_shifts(:,1,3) = (/0,1,0/)
284  tetra_shifts(:,2,3) = (/1,1,0/)
285  tetra_shifts(:,3,3) = (/1,0,1/)
286  tetra_shifts(:,4,3) = (/1,1,1/)
287  tetra_shifts(:,1,4) = (/0,0,0/)
288  tetra_shifts(:,2,4) = (/0,1,0/)
289  tetra_shifts(:,3,4) = (/0,0,1/)
290  tetra_shifts(:,4,4) = (/1,0,1/)
291  tetra_shifts(:,1,5) = (/0,0,1/)
292  tetra_shifts(:,2,5) = (/1,0,1/)
293  tetra_shifts(:,3,5) = (/0,1,0/)
294  tetra_shifts(:,4,5) = (/0,1,1/)
295  tetra_shifts(:,1,6) = (/0,1,0/)
296  tetra_shifts(:,2,6) = (/1,0,1/)
297  tetra_shifts(:,3,6) = (/0,1,1/)
298  tetra_shifts(:,4,6) = (/1,1,1/)
299 
300  ! make full k-point rank arrays
301  call mkkptrank (kpt_fullbz,nkpt_fullbz,kptrank_t)
302 
303  ialltetra = 1
304  do ikpt_full=1,nkpt_fullbz
305    do itetra=1,6
306      do isummit=1,4
307 
308        k1(:) = kpt_fullbz(:,ikpt_full) &
309 &       + tetra_shifts(1,isummit,itetra)*klatt(:,1) &
310 &       + tetra_shifts(2,isummit,itetra)*klatt(:,2) &
311 &       + tetra_shifts(3,isummit,itetra)*klatt(:,3)
312 
313        ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full !
314        call get_rank_1kpt (k1,symrankkpt,kptrank_t)
315        ikpt2 = kptrank_t%invrank(symrankkpt)
316        if (ikpt2 < 1) then
317          errorstring='Error in ranking k-points - exiting with un-initialized tetrahedra.'
318          ierr = 2
319          call destroy_kptrank (kptrank_t)
320          TETRA_ALLOCATE(tetra%tetra_full, (4,2,1))
321          TETRA_ALLOCATE(tetra%tetra_mult, (1))
322          TETRA_ALLOCATE(tetra%tetra_wrap, (3,4,1))
323          TETRA_DEALLOCATE(tetra_full_)
324          TETRA_DEALLOCATE(tetra_mult_)
325          TETRA_DEALLOCATE(tetra_wrap_)
326          return
327        end if
328 
329        ! Store irreducible kpoint equivalent to kpt_fullbz(:,ikpt2)
330        tetra_full_(isummit,1,ialltetra) = indkpt(ikpt2)
331        tetra_full_(isummit,2,ialltetra) = ikpt2
332        shift1 = k1(1)-kpt_fullbz(1,ikpt2)
333        shift2 = k1(2)-kpt_fullbz(2,ikpt2)
334        shift3 = k1(3)-kpt_fullbz(3,ikpt2)
335        if (shift1>0.5d0) then
336          tetra_wrap_(1,isummit,ialltetra) = 1
337        else if (shift1<-0.5d0) then
338          tetra_wrap_(1,isummit,ialltetra) = -1
339        end if
340        if (shift2>0.5d0) then
341          tetra_wrap_(2,isummit,ialltetra) = 1
342        else if (shift2<-0.5d0) then
343          tetra_wrap_(2,isummit,ialltetra) = -1
344        end if
345        if (shift3>0.5d0) then
346          tetra_wrap_(3,isummit,ialltetra) = 1
347        else if (shift3<-0.5d0) then
348          tetra_wrap_(3,isummit,ialltetra) = -1
349        end if
350 
351        ! sort itetra summits
352        ! TODO: replace with sort_int
353        do jsummit=isummit,2,-1
354          if ( tetra_full_(jsummit,1,ialltetra)  <  tetra_full_(jsummit-1,1,ialltetra) ) then
355            itmp = tetra_full_(jsummit,1,ialltetra)
356            tetra_full_(jsummit,1,ialltetra) = tetra_full_(jsummit-1,1,ialltetra)
357            tetra_full_(jsummit-1,1,ialltetra) = itmp
358            itmp = tetra_full_(jsummit,2,ialltetra)
359            tetra_full_(jsummit,2,ialltetra) = tetra_full_(jsummit-1,2,ialltetra)
360            tetra_full_(jsummit-1,2,ialltetra) = itmp
361            ! keep fullbz_kpt tetrahedra points in same order
362            itmp = tetra_wrap_(1,jsummit,ialltetra)
363            tetra_wrap_(1,jsummit,ialltetra) = tetra_wrap_(1,jsummit-1,ialltetra)
364            tetra_wrap_(1,jsummit-1,ialltetra) = itmp
365            itmp = tetra_wrap_(2,jsummit,ialltetra)
366            tetra_wrap_(2,jsummit,ialltetra) = tetra_wrap_(2,jsummit-1,ialltetra)
367            tetra_wrap_(2,jsummit-1,ialltetra) = itmp
368            itmp = tetra_wrap_(1,jsummit,ialltetra)
369            tetra_wrap_(3,jsummit,ialltetra) = tetra_wrap_(3,jsummit-1,ialltetra)
370            tetra_wrap_(3,jsummit-1,ialltetra) = itmp
371          end if
372        end do ! jsummit
373 
374      end do ! isummit
375 
376      if (ialltetra > mtetra) then
377        write (errorstring, '(3a,i0,a,i0)' ) &
378         'init_tetra: BUG - ',&
379         ' ialltetra > mtetra ',&
380         ' ialltetra=  ',ialltetra,', mtetra= ',mtetra
381        ierr = 1
382        return
383      end if
384      ialltetra = ialltetra+1
385    end do ! itetra
386  end do ! ikpt_full
387 
388  call destroy_kptrank (kptrank_t)
389 
390  rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) &
391 & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) &
392 & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3)))
393 
394  ! Volume of all tetrahedra should be the same as that of tetra 1
395  ! this is the volume of 1 tetrahedron, should be coherent with notation in Lehmann & Taut
396  k1(:) = gprimd(:,1)*klatt(1,1) +  gprimd(:,2)*klatt(2,1) +  gprimd(:,3)*klatt(3,1)
397  k2(:) = gprimd(:,1)*klatt(1,2) +  gprimd(:,2)*klatt(2,2) +  gprimd(:,3)*klatt(3,2)
398  k3(:) = gprimd(:,1)*klatt(1,3) +  gprimd(:,2)*klatt(2,3) +  gprimd(:,3)*klatt(3,3)
399  tetra%vv  = abs (k1(1)*(k2(2)*k3(3)-k2(3)*k3(2)) &
400 & -k1(2)*(k2(1)*k3(3)-k2(3)*k3(1)) &
401 & +k1(3)*(k2(1)*k3(2)-k2(2)*k3(1))) / 6.d0 / rcvol
402 
403  ! eliminate equivalent tetrahedra by symmetry and account for them in multiplicity tetra_mult
404  tetra%ntetra = mtetra
405 
406  ! FIXME: could we replace this with a ranking algorithm to avoid the O(tetra%ntetra^2) step? For example:
407  ! get tetrahedron rank - problem too many combinations in principle = nkpt_irred^4 - only a few used in practice
408  ! sort ranks and keep indices
409 
410  ! make hash table = tetra_full_(1)*nkptirred**3+tetra_full_(2)*nkptirred**2+tetra_full_(3)*nkptirred**1+tetra_full_(4)
411 
412  hashfactor = 100.d0 ! *acos(-1.d0) ! 100 pi should be far from an integer...
413  TETRA_ALLOCATE(tetra_hash, (tetra%ntetra))
414  TETRA_ALLOCATE(reforder, (tetra%ntetra))
415 
416  !MG: In principle the order of the indices should not matter.
417  do ialltetra=1, tetra%ntetra
418    tetra_hash(ialltetra) = tetra_full_(1,1,ialltetra)*hashfactor**3+&
419 &      tetra_full_(2,1,ialltetra)*hashfactor**2+&
420 &      tetra_full_(3,1,ialltetra)*hashfactor**1+&
421 &      tetra_full_(4,1,ialltetra)
422    reforder(ialltetra) = ialltetra
423  end do
424 
425  call sort_tetra(tetra%ntetra, tetra_hash, reforder, tol6)
426 
427  ! determine number of tetra after reduction
428  TETRA_ALLOCATE(irred_itetra, (tetra%ntetra))
429  jalltetra = 1
430  irred_itetra(1) = 1
431  do ialltetra=2, tetra%ntetra
432    if (abs(tetra_hash(ialltetra)-tetra_hash(ialltetra-1)) > tol6) then
433      ! found a new series
434      jalltetra = jalltetra + 1
435    end if
436    irred_itetra(ialltetra) = jalltetra
437  end do
438 
439  ! reset number of tetra
440  ntetra_irred = jalltetra
441 
442  ! allocate definitive tetra arrays and transfer to new arrays
443  TETRA_ALLOCATE(tetra%tetra_full, (4,2,ntetra_irred))
444  TETRA_ALLOCATE(tetra%tetra_mult, (ntetra_irred))
445  TETRA_ALLOCATE(tetra%tetra_wrap, (3,4,ntetra_irred))
446 
447  ! eliminate equal rank tetrahedra and accumulate multiplicity into first one
448  tetra%tetra_full = 0
449  tetra%tetra_mult = 0
450  tetra%tetra_wrap = 0
451  jalltetra = 1
452  tetra%tetra_full(:,:,1) = tetra_full_(:,:,reforder(1))
453  tetra%tetra_mult(1) = 1
454  tetra%tetra_wrap(:,:,1) = tetra_wrap_(:,:,reforder(1))
455  do ialltetra=2, tetra%ntetra
456    ! TODO: check if tolerance is adapted
457    if (abs(tetra_hash(ialltetra)-tetra_hash(ialltetra-1)) > tol6) then
458      ! found a new series
459      jalltetra = jalltetra + 1
460      tetra%tetra_full(:,:,jalltetra) = tetra_full_(:,:,reforder(ialltetra))
461      tetra%tetra_wrap(:,:,jalltetra) = tetra_wrap_(:,:,reforder(ialltetra))
462      tetra%tetra_mult(jalltetra) = 1
463    else
464      ! TODO: add real check that the tetra are equivalent...
465      ! otherwise increment jalltetra here as well, generate new series?
466      tetra%tetra_mult(jalltetra) = tetra%tetra_mult(jalltetra) + tetra_mult_(reforder(ialltetra))
467      !tetra_mult_(reforder(ialltetra)) = 0
468    end if
469  end do
470 
471  ! reset of ntetra for final version after checks and debu
472  tetra%ntetra = ntetra_irred
473 
474  TETRA_DEALLOCATE(tetra_hash)
475  TETRA_DEALLOCATE(reforder)
476  TETRA_DEALLOCATE(irred_itetra)
477  TETRA_DEALLOCATE(tetra_full_)
478  TETRA_DEALLOCATE(tetra_mult_)
479  TETRA_DEALLOCATE(tetra_wrap_)
480 
481 end subroutine init_tetra

m_tetrahedron/sort_tetra [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

  sort_tetra

FUNCTION

  Sort double precision array list(n) into ascending numerical order using Heapsort
  algorithm, while making corresponding rearrangement of the integer
  array iperm. Consider that two double precision numbers
  within tolerance tol are equal.

INPUTS

  n        intent(in)    dimension of the list
  tol      intent(in)    numbers within tolerance are equal
  list(n)  intent(inout) list of double precision numbers to be sorted
  iperm(n) intent(inout) iperm(i)=i (very important)

OUTPUT

  list(n)  sorted list
  iperm(n) index of permutation given the right ascending order

PARENTS

      m_tetrahedron

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

1338 subroutine sort_tetra(n,list,iperm,tol)
1339 
1340 
1341 !This section has been created automatically by the script Abilint (TD).
1342 !Do not modify the following lines by hand.
1343 #undef ABI_FUNC
1344 #define ABI_FUNC 'sort_tetra'
1345 !End of the abilint section
1346 
1347  implicit none
1348 
1349  integer, intent(in) :: n
1350  integer, intent(inout) :: iperm(n)
1351  real(dp_), intent(inout) :: list(n)
1352  real(dp_), intent(in) :: tol
1353 
1354  integer :: l,ir,iap,i,j
1355  real(dp_) :: ap
1356  character(len=500) :: msg
1357 
1358  if (n==1) then
1359    ! Accomodate case of array of length 1: already sorted!
1360    return
1361  else if (n<1) then
1362   ! Should not call with n<1
1363   write(msg,1000) n
1364   1000  format(/,' sort_tetra has been called with array length n=',i12,/, &
1365 &  ' having a value less than 1.  This is not allowed.')
1366   TETRA_ERROR(msg)
1367 
1368  else ! n>1
1369 
1370   ! Conduct the usual sort
1371   l=n/2+1
1372   ir=n
1373 
1374   do   ! Infinite do-loop
1375    if (l>1) then
1376     l=l-1
1377     ap=list(l)
1378     iap=iperm(l)
1379 
1380    else ! l<=1
1381     ap=list(ir)
1382     iap=iperm(ir)
1383     list(ir)=list(1)
1384     iperm(ir)=iperm(1)
1385     ir=ir-1
1386 
1387     if (ir==1) then
1388      list(1)=ap
1389      iperm(1)=iap
1390      exit   ! This is the end of this algorithm
1391     end if
1392    end if ! l>1
1393 
1394    i=l
1395    j=l+l
1396 
1397    do while (j<=ir)
1398     if (j<ir) then
1399      if ( list(j)<list(j+1)-tol .or.  &
1400 &        (list(j)<list(j+1)+tol.and.iperm(j)<iperm(j+1))) j=j+1
1401     endif
1402     if (ap<list(j)-tol .or. (ap<list(j)+tol.and.iap<iperm(j))) then
1403      list(i)=list(j)
1404      iperm(i)=iperm(j)
1405      i=j
1406      j=j+j
1407     else
1408      j=ir+1
1409     end if
1410    enddo
1411 
1412    list(i)=ap
1413    iperm(i)=iap
1414 
1415   enddo ! End infinite do-loop
1416 
1417  end if ! n>1
1418 
1419 end subroutine sort_tetra

m_tetrahedron/split_work [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

  split_work

FUNCTION

  Splits the number of tasks, ntasks, among nprocs processors. Used for the MPI parallelization of simple loops.

INPUTS

  ntasks=number of tasks
  comm=MPI communicator.

OUTPUT

  nprocs=Number of MPI processes in the communicator.
  my_start,my_stop= indices defining the initial and final task for this processor
  ierr=Exit status.

NOTES

  If nprocs>ntasks then :
    my_start=ntasks+1
    my_stop=ntask

  In this particular case, loops of the form

  do ii=my_start,my_stop
   ...
  end do

  are not executed. Moreover allocation such as foo(my_start:my_stop) will generate a zero-sized array.

PARENTS

      m_tetrahedron

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

1493 subroutine split_work(ntasks,comm,nprocs,my_start,my_stop,ierr)
1494 
1495 
1496 !This section has been created automatically by the script Abilint (TD).
1497 !Do not modify the following lines by hand.
1498 #undef ABI_FUNC
1499 #define ABI_FUNC 'split_work'
1500 !End of the abilint section
1501 
1502  implicit none
1503 
1504 !Arguments ------------------------------------
1505  integer,intent(in)  :: ntasks,comm
1506  integer,intent(out) :: nprocs,my_start,my_stop,ierr
1507 
1508 !Local variables-------------------------------
1509  integer :: res,my_rank,block_p1,block,mpierr
1510 
1511 ! *************************************************************************
1512 
1513  nprocs = 1; my_start = 1; my_stop = ntasks; ierr = 1
1514 #ifdef HAVE_MPI
1515  call MPI_COMM_SIZE(comm,nprocs,mpierr); if (mpierr /= MPI_SUCCESS) return
1516  call MPI_COMM_RANK(comm,my_rank,mpierr); if (mpierr /= MPI_SUCCESS) return
1517 
1518  block   = ntasks/nprocs
1519  res     = MOD(ntasks,nprocs)
1520  block_p1= block+1
1521 
1522  if (my_rank<res) then
1523    my_start =  my_rank   *block_p1+1
1524    my_stop  = (my_rank+1)*block_p1
1525  else
1526    my_start = res*block_p1 + (my_rank-res  )*block + 1
1527    my_stop  = res*block_p1 + (my_rank-res+1)*block
1528  end if
1529 #endif
1530  ierr = 0
1531 
1532 end subroutine split_work

m_tetrahedron/t_tetrahedron [ Types ]

[ Top ] [ m_tetrahedron ] [ Types ]

NAME

 t_tetrahedron

FUNCTION

 tetrahedron geometry object

SOURCE

 70 type, public :: t_tetrahedron
 71 
 72   integer :: ntetra = 0
 73   ! Number of tetrahedra
 74 
 75   real(dp_)  :: vv
 76   ! volume of the tetrahedra
 77 
 78   real(dp_) :: klatt(3, 3)
 79   ! reciprocal of lattice vectors for full kpoint grid
 80 
 81   integer,allocatable :: tetra_full(:,:,:)
 82   !(4,2,ntetra)
 83   ! For each tetra
 84   !   (:,1,itetra) indices of the vertex in IBZ (symmetrical image)
 85   !   (:,2,itetra) indices of the vertexes in the BZ
 86 
 87   integer,allocatable :: tetra_mult(:)
 88   !(ntetra)
 89   ! multiplicity of each irred tetrahedron
 90 
 91   integer,allocatable :: tetra_wrap(:,:,:)
 92   !(3,4,ntetra)
 93   ! flag to wrap tetrahedron summit into IBZ
 94 
 95 end type t_tetrahedron
 96 
 97 public :: init_tetra           ! Initialize the object
 98                                ! (see also the high-level interface tetra_from_kptrlatt provided by m_kpts).
 99 public :: get_tetra_weight     ! Calculate integration weights and their derivatives. shape (nkpt,nene)
100 public :: tetra_blochl_weights ! Same as in get_tetra_weight but weights have shape (nene,nkpt)
101 public :: get_dbl_tetra_weight ! Calculate integration weights for double tetrahedron
102                                !  integration of delta functions (NB these correspond
103                                !  to the derivative terms in normal tetrahedron).
104 public :: destroy_tetra        ! Free memory.
105 
106 public :: tetra_write          ! Write text file with tetra info.
107 
108 public :: tetralib_has_mpi     ! Return True if the library has been compiled with MPI support
109 public :: tetra_get_onewk      ! Calculate integration weights and their derivatives for a single k-point in the IBZ.

m_tetrahedron/tetra_blochl_weights [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 tetra_blochl_weights

FUNCTION

 calculate integration weights and their derivatives from Blochl et al PRB 49 16223 [[cite:Bloechl1994a]]
 Same API as get_tetra_weight but weights here have shape (nene, nkpt)

PARENTS

      m_fstab,m_phonons,m_tetrahedron

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

682 subroutine tetra_blochl_weights(tetra,eigen_in,enemin,enemax,max_occ,nene,nkpt,&
683   bcorr,tweight_t,dtweightde_t,comm)
684 
685 
686 !This section has been created automatically by the script Abilint (TD).
687 !Do not modify the following lines by hand.
688 #undef ABI_FUNC
689 #define ABI_FUNC 'tetra_blochl_weights'
690 !End of the abilint section
691 
692  implicit none
693 
694 !Arguments ------------------------------------
695 !scalars
696  integer,intent(in) :: nene,nkpt,bcorr,comm
697  type(t_tetrahedron), intent(in) :: tetra
698  real(dp_) ,intent(in) :: enemax,enemin,max_occ
699 !arrays
700  real(dp_) ,intent(in) :: eigen_in(nkpt)
701  real(dp_) ,intent(out) :: dtweightde_t(nene,nkpt),tweight_t(nene,nkpt)
702 
703 !Local variables-------------------------------
704 !scalars
705  integer :: itetra,nprocs,my_start,my_stop,ierr
706  real(dp_) :: deltaene,volconst,volconst_mult
707 !arrays
708  integer :: ind_ibz(4)
709  real(dp_) :: eigen_1tetra(4)
710  real(dp_), allocatable :: tweight_tmp(:,:),dtweightde_tmp(:,:),buffer(:,:)
711 
712 ! *********************************************************************
713 
714  TETRA_ALLOCATE(tweight_tmp, (nene, 4))
715  TETRA_ALLOCATE(dtweightde_tmp, (nene, 4))
716  tweight_t = zero; dtweightde_t = zero
717 
718  volconst = tetra%vv/4.d0
719  if (nene <= 1) then
720    TETRA_ERROR('tetra_blochl_weights: nene must be at least 2')
721  else
722    deltaene = (enemax-enemin) / (nene-1)
723  end if
724 
725  call split_work(tetra%ntetra,comm,nprocs,my_start,my_stop,ierr)
726  if (ierr /= 0) TETRA_ERROR("Error in MPI layer")
727 
728  ! for each tetrahedron
729  do itetra=my_start,my_stop
730    tweight_tmp = zero
731    dtweightde_tmp = zero
732 
733    volconst_mult = max_occ*volconst*float(tetra%tetra_mult(itetra))
734 
735    ! Here we need the original ordering to reference the correct irred kpoints
736    ind_ibz(1) = tetra%tetra_full(1,1,itetra)
737    ind_ibz(2) = tetra%tetra_full(2,1,itetra)
738    ind_ibz(3) = tetra%tetra_full(3,1,itetra)
739    ind_ibz(4) = tetra%tetra_full(4,1,itetra)
740 
741    ! Sort energies before calling get_onetetra_
742    eigen_1tetra(1) = eigen_in(ind_ibz(1))
743    eigen_1tetra(2) = eigen_in(ind_ibz(2))
744    eigen_1tetra(3) = eigen_in(ind_ibz(3))
745    eigen_1tetra(4) = eigen_in(ind_ibz(4))
746    call sort_tetra(4,eigen_1tetra,ind_ibz,tol14)
747 
748    call get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,nene,bcorr,tweight_tmp,dtweightde_tmp)
749 
750    ! NOTE: the following blas calls are not working systematically, or do not give speed ups, strange...
751    !call daxpy (nene, 1.d0,    tweight_tmp(:,1), 1,    tweight_t(:,ind_ibz(1)), 1)
752    !call daxpy (nene, 1.d0,    tweight_tmp(:,2), 1,    tweight_t(:,ind_ibz(2)), 1)
753    !call daxpy (nene, 1.d0,    tweight_tmp(:,3), 1,    tweight_t(:,ind_ibz(3)), 1)
754    !call daxpy (nene, 1.d0,    tweight_tmp(:,4), 1,    tweight_t(:,ind_ibz(4)), 1)
755    !call daxpy (nene, 1.d0, dtweightde_tmp(:,1), 1, dtweightde_t(:,ind_ibz(1)), 1)
756    !call daxpy (nene, 1.d0, dtweightde_tmp(:,2), 1, dtweightde_t(:,ind_ibz(2)), 1)
757    !call daxpy (nene, 1.d0, dtweightde_tmp(:,3), 1, dtweightde_t(:,ind_ibz(3)), 1)
758    !call daxpy (nene, 1.d0, dtweightde_tmp(:,4), 1, dtweightde_t(:,ind_ibz(4)), 1)
759    tweight_t(:,ind_ibz(1)) = tweight_t(:,ind_ibz(1)) + tweight_tmp(:,1)
760    tweight_t(:,ind_ibz(2)) = tweight_t(:,ind_ibz(2)) + tweight_tmp(:,2)
761    tweight_t(:,ind_ibz(3)) = tweight_t(:,ind_ibz(3)) + tweight_tmp(:,3)
762    tweight_t(:,ind_ibz(4)) = tweight_t(:,ind_ibz(4)) + tweight_tmp(:,4)
763    dtweightde_t(:,ind_ibz(1)) = dtweightde_t(:,ind_ibz(1)) + dtweightde_tmp(:,1)
764    dtweightde_t(:,ind_ibz(2)) = dtweightde_t(:,ind_ibz(2)) + dtweightde_tmp(:,2)
765    dtweightde_t(:,ind_ibz(3)) = dtweightde_t(:,ind_ibz(3)) + dtweightde_tmp(:,3)
766    dtweightde_t(:,ind_ibz(4)) = dtweightde_t(:,ind_ibz(4)) + dtweightde_tmp(:,4)
767  end do ! itetra
768 
769  TETRA_DEALLOCATE(tweight_tmp)
770  TETRA_DEALLOCATE(dtweightde_tmp)
771 
772  if (nprocs > 1) then
773 #ifdef HAVE_MPI
774    TETRA_ALLOCATE(buffer, (nene, nkpt))
775    call MPI_ALLREDUCE(tweight_t,buffer,nene*nkpt,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
776    tweight_t = buffer
777 
778    call MPI_ALLREDUCE(dtweightde_t,buffer,nene*nkpt,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr)
779    dtweightde_t = buffer
780    TETRA_DEALLOCATE(buffer)
781 #endif
782  end if
783 
784 end subroutine tetra_blochl_weights

m_tetrahedron/tetra_get_onewk [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 tetra_get_onewk

FUNCTION

 Calculate integration weights and their derivatives for a single k-point in the IBZ.

INPUTS

 tetra<t_tetrahedron>=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.
 nene=number of energies for DOS
 nibz=number of irreducible kpoints
 eigen_ibz(nkibz)=eigenenergies for each k point
 enemin=minimal energy for DOS
 enemax=maximal energy for DOS
 max_occ=maximal occupation number (2 for nsppol=1, 1 for nsppol=2)

OUTPUT

  weights(nene,2) = integration weights for
    Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
    for a given (band, k-point, spin).

PARENTS

      m_ebands,m_epjdos,m_gruneisen,m_phgamma

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

1904 subroutine tetra_get_onewk(tetra,ik_ibz,bcorr,nene,nkibz,eig_ibz,enemin,enemax,max_occ,weights)
1905 
1906 
1907 !This section has been created automatically by the script Abilint (TD).
1908 !Do not modify the following lines by hand.
1909 #undef ABI_FUNC
1910 #define ABI_FUNC 'tetra_get_onewk'
1911 !End of the abilint section
1912 
1913  implicit none
1914 
1915 !Arguments ------------------------------------
1916 !scalars
1917  integer,intent(in) :: ik_ibz,nene,nkibz,bcorr
1918  type(t_tetrahedron), intent(in) :: tetra
1919  real(dp_) ,intent(in) :: enemin,enemax,max_occ
1920 !arrays
1921  real(dp_),intent(in) :: eig_ibz(nkibz)
1922  real(dp_),intent(out) :: weights(nene,2)
1923 
1924 !Local variables-------------------------------
1925 !scalars
1926  integer :: itetra,ii
1927  real(dp_) :: deltaene
1928 !arrays
1929  integer :: ind_ibz(4)
1930  real(dp_) :: tweight_tmp(nene,4),dtweightde_tmp(nene,4)
1931  real(dp_) :: eigen_1tetra(4)
1932 
1933 ! *********************************************************************
1934 
1935  if (nene <= 1) then
1936    TETRA_ERROR('tetra_blochl_weights: nene must be at least 2')
1937  else
1938    deltaene = (enemax-enemin) / (nene-1)
1939  end if
1940 
1941  weights = zero
1942 
1943  ! For each tetrahedron
1944  do itetra=1,tetra%ntetra
1945 
1946    ! Here we need the original ordering to reference the correct irred kpoints
1947    ind_ibz(1) = tetra%tetra_full(1,1,itetra)
1948    ind_ibz(2) = tetra%tetra_full(2,1,itetra)
1949    ind_ibz(3) = tetra%tetra_full(3,1,itetra)
1950    ind_ibz(4) = tetra%tetra_full(4,1,itetra)
1951    ! Cycle if this tetra does not contribute to this k-point.
1952    if (all(ind_ibz /= ik_ibz)) cycle
1953 
1954    ! Sort energies before calling get_onetetra_
1955    eigen_1tetra(1) = eig_ibz(ind_ibz(1))
1956    eigen_1tetra(2) = eig_ibz(ind_ibz(2))
1957    eigen_1tetra(3) = eig_ibz(ind_ibz(3))
1958    eigen_1tetra(4) = eig_ibz(ind_ibz(4))
1959    call sort_tetra(4,eigen_1tetra,ind_ibz,tol14)
1960 
1961    call get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,&
1962 &   nene,bcorr,tweight_tmp,dtweightde_tmp)
1963 
1964    ! Accumulate contributions to ik_ibz
1965    ! (there might be multiple vertexes that map onto ik_ibz)
1966    do ii=1,4
1967      if (ind_ibz(ii) == ik_ibz) then
1968        weights(:,1) = weights(:,1) + dtweightde_tmp(:,ii)
1969        weights(:,2) = weights(:,2) + tweight_tmp(:,ii)
1970      end if
1971    end do
1972  end do ! itetra
1973 
1974 end subroutine tetra_get_onewk

m_tetrahedron/tetra_write [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 tetra_write

FUNCTION

  Write text file with tetra info.

INPUTS

  tetra<t_tetrahedron>=tetrahedron geometry object
  nkibz=Number of k-points in the IBZ used to generate tetra
  kibz(3,nkibz)=Reduced coordinates of the IBZ
  path=Name of output file

OUTPUT

  Output is written to file.

PARENTS

      gstate,wfk_analyze

CHILDREN

      get_onetetra_,sort_tetra

SOURCE

510 subroutine tetra_write(tetra, nkibz, kibz, path)
511 
512 
513 !This section has been created automatically by the script Abilint (TD).
514 !Do not modify the following lines by hand.
515 #undef ABI_FUNC
516 #define ABI_FUNC 'tetra_write'
517 !End of the abilint section
518 
519  implicit none
520 
521 !Arguments ------------------------------------
522 !scalars
523  integer,intent(in) :: nkibz
524  character(len=*),intent(in) :: path
525  type(t_tetrahedron),intent(in) :: tetra
526 !arrays
527  real(dp_),intent(in) :: kibz(3,nkibz)
528 
529 !Local variables-------------------------------
530  integer,parameter :: version=1
531  integer :: ik,it,unt
532 #ifdef HAVE_LIBTETRA_ABINIT
533  character(len=500) :: msg
534 #endif
535 
536 ! *********************************************************************
537 
538 #ifdef HAVE_LIBTETRA_ABINIT
539  if (open_file(file=trim(path),iomsg=msg,newunit=unt,form="formatted",status="unknown",action="write")/=0) then
540    TETRA_ERROR(msg)
541  end if
542 #else
543  open(file=trim(path),newunit=unt,form="formatted",status="unknown",action="write")
544 #endif
545 
546  write(unt,*)version, " # version number"
547 
548  ! Write IBZ
549  write(unt,*)nkibz, " # number of k-points in the IBZ"
550  write(unt,"(a)")"<irreducible_zone>"
551  do ik=1,nkibz
552    write(unt,"(3es22.12)") kibz(:,ik)
553  end do
554  write(unt,"(a)")"</irreducible_zone>"
555 
556  ! Write tetra info
557  write(unt,"(i0,a)")tetra%ntetra, " # number of tetrahedra"
558  write(unt,"(es22.12,a)")tetra%vv, " # tetrahedron volume"
559 
560  write(unt,"(a)")"<tetra_full>"
561  do it=1,tetra%ntetra
562    write(unt,"(8(i0,1x))")tetra%tetra_full(:,:,it)
563  end do
564  write(unt,"(a)")"</tetra_full>"
565 
566  write(unt,"(a)")"<tetra_mult>"
567  do it=1,tetra%ntetra
568    write(unt,"(i0)")tetra%tetra_mult(it)
569  end do
570  write(unt,"(a)")"</tetra_mult>"
571 
572  write(unt,"(a)")"<tetra_wrap>"
573  do it=1,tetra%ntetra
574    write(unt,"(12(i0,1x))")tetra%tetra_wrap(:,:,it)
575  end do
576  write(unt,"(a)")"</tetra_wrap>"
577 
578  close(unt)
579 
580 end subroutine tetra_write

m_tetrahedron/tetralib_has_mpi [ Functions ]

[ Top ] [ m_tetrahedron ] [ Functions ]

NAME

 tetralib_has_mpi

FUNCTION

 Return True if library has been compiled with MPI support

PARENTS

CHILDREN

SOURCE

1437 logical function tetralib_has_mpi() result(ans)
1438 
1439 
1440 !This section has been created automatically by the script Abilint (TD).
1441 !Do not modify the following lines by hand.
1442 #undef ABI_FUNC
1443 #define ABI_FUNC 'tetralib_has_mpi'
1444 !End of the abilint section
1445 
1446   ans = .False.
1447 #ifdef HAVE_MPI
1448   ans = .True.
1449 #endif
1450 
1451 end function tetralib_has_mpi