TABLE OF CONTENTS
- ABINIT/m_tetrahedron
- m_tetrahedron/destroy_tetra
- m_tetrahedron/get_dbl_tetra_weight
- m_tetrahedron/get_onetetra_
- m_tetrahedron/get_tetra_weight
- m_tetrahedron/init_tetra
- m_tetrahedron/sort_tetra
- m_tetrahedron/split_work
- m_tetrahedron/t_tetrahedron
- m_tetrahedron/tetra_blochl_weights
- m_tetrahedron/tetra_get_onewk
- m_tetrahedron/tetra_write
- m_tetrahedron/tetralib_has_mpi
ABINIT/m_tetrahedron [ 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