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_onetetra_wvals
- m_tetrahedron/tetra_get_onewk
- m_tetrahedron/tetra_get_onewk_wvals
- 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-2024 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 carefully the case of degenerate tetrahedron 2) Change API so that we can pass the energy mesh instead of omega_min and omega_max 3) Add table ik_ibz --> tetra_list to avoid cycling inside big loop over ntetra 4) Add options to get only delta and/or theta ?
SOURCE
23 #if defined HAVE_CONFIG_H 24 #include "config.h" 25 #endif 26 27 #include "libtetra.h" 28 29 module m_tetrahedron 30 31 ! make sure stdout is defined, as libtetra.h needs it 32 use, intrinsic :: iso_fortran_env, only : stdin=>input_unit, & 33 stdout=>output_unit, & 34 stderr=>error_unit 35 USE_MEMORY_PROFILING 36 USE_MSG_HANDLING 37 use m_krank 38 #ifdef HAVE_MPI2 39 use mpi 40 #endif 41 #ifdef HAVE_LIBTETRA_ABINIT 42 use m_io_tools, only : open_file 43 use m_xmpi 44 #endif 45 46 implicit none 47 48 #if defined HAVE_MPI1 49 include 'mpif.h' 50 #endif 51 52 private
m_tetrahedron/destroy_tetra [ Functions ]
[ Top ] [ m_tetrahedron ] [ Functions ]
NAME
destroy_tetra
FUNCTION
deallocate tetrahedra pointers if needed
SOURCE
134 subroutine destroy_tetra (tetra) 135 136 type(t_tetrahedron), intent(inout) :: tetra 137 138 if (allocated(tetra%tetra_full)) then 139 TETRA_DEALLOCATE(tetra%tetra_full) 140 end if 141 if (allocated(tetra%tetra_mult)) then 142 TETRA_DEALLOCATE(tetra%tetra_mult) 143 end if 144 if (allocated(tetra%tetra_wrap)) then 145 TETRA_DEALLOCATE(tetra%tetra_wrap) 146 end if 147 if (allocated(tetra%ibz_tetra_count)) then 148 TETRA_DEALLOCATE(tetra%ibz_tetra_count) 149 end if 150 if (allocated(tetra%ibz_tetra_mapping)) then 151 TETRA_DEALLOCATE(tetra%ibz_tetra_mapping) 152 end if 153 154 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
SOURCE
809 subroutine get_dbl_tetra_weight(eigen1_in,eigen2_in,enemin1,enemax1,enemin2,enemax2,& 810 & max_occ,nene1,nene2,nkpt,tetra,tweight,dtweightde, ierr) 811 812 !Arguments ------------------------------------ 813 !scalars 814 integer,intent(in) :: nene1,nene2,nkpt 815 integer,intent(out) :: ierr 816 type(t_tetrahedron), intent(in) :: tetra 817 real(dp),intent(in) :: enemax1,enemin1 818 real(dp),intent(in) :: enemax2,enemin2 819 real(dp),intent(in) :: max_occ 820 !arrays 821 real(dp),intent(in) :: eigen1_in(nkpt) 822 real(dp),intent(in) :: eigen2_in(nkpt) 823 real(dp),intent(out) :: dtweightde(nkpt,nene1,nene2),tweight(nkpt,nene1,nene2) 824 825 !Local variables------------------------------- 826 ! needed for gaussian replacement of Dirac functions 827 ! the three coefficients of the DOS as quadratic form, 828 ! in the interval [eig(ikpt-1), eig(ikpt)] 829 ! for ikpt = 1 we add a point below eigen(1) which doesnt 830 ! contribute to the DOS in any tetrahedron 831 !scalars 832 integer :: ieps1,ieps2,itetra 833 integer :: nn1_1,nn1_2,nn1_3,nn1_4 834 integer :: nn2_1,nn2_2,nn2_3 835 integer :: ind_a(3), ind_b(3), ind_c(3) 836 real(dp) :: deltaene1,eps1 837 real(dp) :: deltaene2,eps2 838 ! real(dp) :: gau_prefactor,gau_width,gau_width2 839 real(dp) :: epsilon1(4,4) 840 real(dp) :: epsilon2(4,4) 841 real(dp) :: inv_epsilon1(4,4) 842 real(dp) :: aa(3),bb(3),cc(3) 843 real(dp) :: delaa(3),delbb(3),delcc(3) 844 real(dp) :: delaa0,delbb0,delcc0 845 real(dp) :: inv_delaa(3),inv_delbb(3),inv_delcc(3) 846 real(dp) :: deleps1, deleps2 847 real(dp) :: inv_deleps1 848 real(dp) :: dccde1, dccde1_pre 849 real(dp) :: volconst,volconst_mult 850 real(dp) :: ii0, ii1, ii3 851 !arrays 852 integer :: ind_k(4) 853 real(dp), allocatable :: tweight_tmp(:,:,:) 854 real(dp), allocatable :: dtweightde_tmp(:,:,:) 855 real(dp) :: eigen1_1tetra(4) 856 real(dp) :: eigen2_1tetra(4) 857 858 ! ********************************************************************* 859 860 ierr = 0 861 if (nene1 <= 1 .or. nene2 <= 1) then 862 !'get_dbl_tetra_weight: nene must be at least 2' 863 ierr = 1 864 return 865 else 866 deltaene1 = (enemax1-enemin1) / (nene1-1) 867 deltaene2 = (enemax2-enemin2) / (nene2-1) 868 end if 869 870 TETRA_ALLOCATE(tweight_tmp, (4, nene2, nene1)) 871 TETRA_ALLOCATE(dtweightde_tmp, (4, nene2, nene1)) 872 873 !print *, "warning: for the moment, heaviside weights are 0. The delta function / DOS weights are the only ones calculated " 874 875 volconst = tetra%vv/4.d0 876 877 ! for each tetrahedron 878 do itetra=1,tetra%ntetra 879 ! these are for 1 tetrahedron only. 880 tweight_tmp = zero 881 dtweightde_tmp = zero 882 883 volconst_mult = max_occ*volconst*dble(tetra%tetra_mult(itetra)) 884 885 ! Here we need the original ordering to reference the correct irred kpoints 886 ! ind_k refers to the index in the full k list of the summits of the present tetrahedra 887 ! we can forget the order of the summits within the tetrahedron, because eigen1 fixes that 888 ! order with its increasing value 889 ind_k(1) = tetra%tetra_full(1,1,itetra) 890 ind_k(2) = tetra%tetra_full(2,1,itetra) 891 ind_k(3) = tetra%tetra_full(3,1,itetra) 892 ind_k(4) = tetra%tetra_full(4,1,itetra) 893 eigen1_1tetra(1) = eigen1_in(ind_k(1)) 894 eigen1_1tetra(2) = eigen1_in(ind_k(2)) 895 eigen1_1tetra(3) = eigen1_in(ind_k(3)) 896 eigen1_1tetra(4) = eigen1_in(ind_k(4)) 897 call sort_tetra(4,eigen1_1tetra,ind_k,tol14) 898 899 ! re-sort eigen2 values according to order chosen for eigen1. Eigen2 are _not_ in increasing order! 900 eigen2_1tetra(1) = eigen2_in(ind_k(1)) 901 eigen2_1tetra(2) = eigen2_in(ind_k(2)) 902 eigen2_1tetra(3) = eigen2_in(ind_k(3)) 903 eigen2_1tetra(4) = eigen2_in(ind_k(4)) 904 905 ! the epsilons are energy differences for the two eigenvalue sets 906 epsilon1 = zero 907 epsilon2 = zero 908 do ieps1 = 1, 4 909 do ieps2 = ieps1+1, 4 910 epsilon1(ieps1,ieps2) = eigen1_1tetra(ieps1)-eigen1_1tetra(ieps2) 911 epsilon1(ieps2,ieps1) = -epsilon1(ieps1,ieps2) 912 epsilon2(ieps1,ieps2) = eigen2_1tetra(ieps1)-eigen2_1tetra(ieps2) 913 epsilon2(ieps2,ieps1) = -epsilon2(ieps1,ieps2) 914 end do 915 end do 916 917 ! we precalculate the inverses to avoid doing tons of divisions in the energy loops below 918 ! Allen formulae only require the inverses of the differences of eigen1 + the a b c below 919 inv_epsilon1 = zero 920 do ieps1 = 1, 4 921 do ieps2 = ieps1+1, 4 922 if (abs(epsilon1(ieps1,ieps2)) > tol6) then 923 inv_epsilon1(ieps1,ieps2) = 1.d0 / epsilon1(ieps1,ieps2) 924 inv_epsilon1(ieps2,ieps1) = -inv_epsilon1(ieps1,ieps2) 925 end if 926 end do 927 end do 928 929 ! these bounds determine the intervals for Omega in Allen paper, and cases A, B, C 930 nn1_1 = int((eigen1_1tetra(1)-enemin1)/deltaene1)+1 931 nn1_2 = int((eigen1_1tetra(2)-enemin1)/deltaene1)+1 932 nn1_3 = int((eigen1_1tetra(3)-enemin1)/deltaene1)+1 933 nn1_4 = int((eigen1_1tetra(4)-enemin1)/deltaene1)+1 934 935 nn1_1 = max(1,nn1_1) 936 nn1_1 = min(nn1_1,nene1) 937 nn1_2 = max(1,nn1_2) 938 nn1_2 = min(nn1_2,nene1) 939 nn1_3 = max(1,nn1_3) 940 nn1_3 = min(nn1_3,nene1) 941 nn1_4 = max(1,nn1_4) 942 nn1_4 = min(nn1_4,nene1) 943 944 ! calculate Allen a_i b_i and c_i parameters 945 ! sort the a_i b_i c_i 946 ! 947 ! NOTE: indices here go from 1 to 4 instead of 0 to 3 as in Allen... 948 aa(1) = epsilon2(2,1) * inv_epsilon1(2,1) 949 aa(2) = epsilon2(3,1) * inv_epsilon1(3,1) 950 aa(3) = epsilon2(4,1) * inv_epsilon1(4,1) 951 ind_a = (/2,3,4/) 952 call sort_tetra(3,aa,ind_a,tol14) 953 ! aa are now in order a_s a_m a_l !!! Preserve the hash function ind_a to order the positions of k below 954 delaa(1) = aa(2)-aa(1) 955 delaa(2) = aa(3)-aa(1) 956 delaa(3) = aa(3)-aa(2) 957 inv_delaa = zero 958 if(delaa(1)> tol6) inv_delaa(1)= 1.0d0 / delaa(1) 959 if(delaa(2)> tol6) inv_delaa(2)= 1.0d0 / delaa(2) 960 if(delaa(3)> tol6) inv_delaa(3)= 1.0d0 / delaa(3) 961 962 bb(1) = epsilon2(1,2) * inv_epsilon1(1,2) 963 bb(2) = epsilon2(3,2) * inv_epsilon1(3,2) 964 bb(3) = epsilon2(4,2) * inv_epsilon1(4,2) 965 ind_b = (/1,3,4/) 966 call sort_tetra(3,bb,ind_b,tol14) 967 delbb(1) = bb(2)-bb(1) 968 delbb(2) = bb(3)-bb(1) 969 delbb(3) = bb(3)-bb(2) 970 inv_delbb = zero 971 if(delbb(1)> tol6) inv_delbb(1)= 1.0d0 / delbb(1) 972 if(delbb(2)> tol6) inv_delbb(2)= 1.0d0 / delbb(2) 973 if(delbb(3)> tol6) inv_delbb(3)= 1.0d0 / delbb(3) 974 975 cc(1) = epsilon2(1,4) * inv_epsilon1(1,4) 976 cc(2) = epsilon2(2,4) * inv_epsilon1(2,4) 977 cc(3) = epsilon2(3,4) * inv_epsilon1(3,4) 978 ind_c = (/1,2,3/) 979 call sort_tetra(3,cc,ind_c,tol14) 980 delcc(1) = cc(2)-cc(1) 981 delcc(2) = cc(3)-cc(1) 982 delcc(3) = cc(3)-cc(2) 983 inv_delcc = zero 984 if(delcc(1)> tol6) inv_delcc(1)= 1.0d0 / delcc(1) 985 if(delcc(2)> tol6) inv_delcc(2)= 1.0d0 / delcc(2) 986 if(delcc(3)> tol6) inv_delcc(3)= 1.0d0 / delcc(3) 987 988 !---------------------------------------------------------------------- 989 ! start main loop A B C over eps1 990 !---------------------------------------------------------------------- 991 992 ! 993 ! interval enemin1 < eps1 < e1 nothing to do 994 ! 995 ! 996 ! interval e1 < eps1 < e3 CASE A in Allen + first term in B 997 ! 998 ! NB: eps1 is not updated inside the loop, only between the loops 999 eps1 = enemin1+nn1_1*deltaene1 1000 deleps1 = eps1-eigen1_1tetra(1) ! this is Omega - omega_0 1001 dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(2,1)*inv_epsilon1(3,1)*inv_epsilon1(4,1) 1002 1003 ! note we go to nn1_3 1004 do ieps1=nn1_1+1,nn1_3 1005 1006 dccde1 = dccde1_pre * deleps1 ! this is f_0(Omega)*6*v 1007 1008 ! at fixed ieps1 we can find the pivot indices for the ieps2 loop 1009 nn2_1 = int((eigen2_1tetra(1)+deleps1*aa(1) -enemin2)/deltaene2)+1 1010 nn2_2 = int((eigen2_1tetra(1)+deleps1*aa(2) -enemin2)/deltaene2)+1 1011 nn2_3 = int((eigen2_1tetra(1)+deleps1*aa(3) -enemin2)/deltaene2)+1 1012 1013 nn2_1 = max(1,nn2_1) 1014 nn2_1 = min(nn2_1,nene2) 1015 nn2_2 = max(1,nn2_2) 1016 nn2_2 = min(nn2_2,nene2) 1017 nn2_3 = max(1,nn2_3) 1018 nn2_3 = min(nn2_3,nene2) 1019 1020 inv_deleps1 = 1.0d0 / deleps1 1021 1022 eps2 = enemin2+nn2_1*deltaene2 ! this is E 1023 deleps2 = eps2 - eigen2_1tetra(1) ! this is E-epsilon_0 1024 1025 !----------------------------------------------------------------------- 1026 ! This is case AI 1027 !----------------------------------------------------------------------- 1028 do ieps2 = nn2_1+1, nn2_2 1029 ! calculate running value of del "a" = a-a_s: first term should really mix eps1 and eps2 1030 delaa0 = deleps2*inv_deleps1 - aa(1) ! a - a_s 1031 1032 ii0 = dccde1*delaa0*inv_delaa(1)*inv_delaa(2) ! this is I_0(Omega E) 1033 1034 dtweightde_tmp(1,ieps2,ieps1) = dtweightde_tmp(1,ieps2,ieps1) + & 1035 & ii0*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_a(1),1)* & 1036 & (-2.0d0 + delaa0*inv_delaa(1)*epsilon1(ind_a(2),ind_a(1))*inv_epsilon1(ind_a(2),1) & 1037 & + delaa0*inv_delaa(2)*epsilon1(ind_a(3),ind_a(1))*inv_epsilon1(ind_a(3),1))) 1038 dtweightde_tmp(ind_a(1),ieps2,ieps1) = dtweightde_tmp(ind_a(1),ieps2,ieps1) + & 1039 & ii0*0.5d0*deleps1*inv_epsilon1(ind_a(1),1)*(2.0d0 - delaa0*inv_delaa(1) - delaa0*inv_delaa(2)) 1040 dtweightde_tmp(ind_a(2),ieps2,ieps1) = dtweightde_tmp(ind_a(2),ieps2,ieps1) + & 1041 & ii0*0.5d0*delaa0*inv_delaa(1)*deleps1*inv_epsilon1(ind_a(2),1) 1042 dtweightde_tmp(ind_a(3),ieps2,ieps1) = dtweightde_tmp(ind_a(3),ieps2,ieps1) + & 1043 & ii0*0.5d0*delaa0*inv_delaa(2)*deleps1*inv_epsilon1(ind_a(3),1) 1044 deleps2 = deleps2 + deltaene2 1045 end do 1046 1047 1048 eps2 = enemin2+nn2_2*deltaene2 ! this is E 1049 deleps2 = eps2 - eigen2_1tetra(1) ! E-E_0 1050 1051 !----------------------------------------------------------------------- 1052 ! This is case AII 1053 !----------------------------------------------------------------------- 1054 do ieps2 = nn2_2+1, nn2_3 1055 ! calculate running value of del "a" = a_l-a: first term should really mix eps1 and eps2 1056 delaa0 = aa(3) - deleps2*inv_deleps1 ! a_l - a 1057 1058 ii0 = dccde1*delaa0*inv_delaa(3)*inv_delaa(2) ! this is I_0(Omega E) 1059 1060 dtweightde_tmp(1,ieps2,ieps1) = dtweightde_tmp(1,ieps2,ieps1) + & 1061 & ii0*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_a(3),1)* & 1062 & (-2.0d0 + delaa0*inv_delaa(3)*epsilon1(ind_a(2),ind_a(3))*inv_epsilon1(ind_a(2),1) & 1063 & + delaa0*inv_delaa(2)*epsilon1(ind_a(1),ind_a(3))*inv_epsilon1(ind_a(1),1))) 1064 dtweightde_tmp(ind_a(3),ieps2,ieps1) = dtweightde_tmp(ind_a(3),ieps2,ieps1) + & 1065 & ii0*0.5d0*deleps1*inv_epsilon1(ind_a(3),1)*(2.0d0 - delaa0*inv_delaa(3) - delaa0*inv_delaa(2)) 1066 dtweightde_tmp(ind_a(2),ieps2,ieps1) = dtweightde_tmp(ind_a(2),ieps2,ieps1) + & 1067 & ii0*0.5d0*delaa0*inv_delaa(3)*deleps1*inv_epsilon1(ind_a(2),1) 1068 dtweightde_tmp(ind_a(1),ieps2,ieps1) = dtweightde_tmp(ind_a(1),ieps2,ieps1) + & 1069 & ii0*0.5d0*delaa0*inv_delaa(2)*deleps1*inv_epsilon1(ind_a(1),1) 1070 1071 deleps2 = deleps2 + deltaene2 1072 end do 1073 deleps1 = deleps1 + deltaene1 1074 end do 1075 ! 1076 ! interval e2 < eps < e3 1077 ! 1078 eps1 = eps1 + (nn1_2-nn1_1)*deltaene1 1079 1080 deleps1 = eps1-eigen1_1tetra(2) ! Omega - omega_1 1081 1082 dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(2,1)*inv_epsilon1(3,2)*inv_epsilon1(4,2) ! f1 function 1083 do ieps1=nn1_2+1,nn1_3 1084 1085 dccde1 = dccde1_pre * deleps1 ! f2(Omega) * 6 * v 1086 1087 ! at fixed ieps1 we can find the pivot indices for the ieps2 loop 1088 nn2_1 = int((eigen2_1tetra(2)+deleps1*bb(1) -enemin2)/deltaene2)+1 1089 nn2_2 = int((eigen2_1tetra(2)+deleps1*bb(2) -enemin2)/deltaene2)+1 1090 nn2_3 = int((eigen2_1tetra(2)+deleps1*bb(3) -enemin2)/deltaene2)+1 1091 1092 nn2_1 = max(1,nn2_1) 1093 nn2_1 = min(nn2_1,nene2) 1094 nn2_2 = max(1,nn2_2) 1095 nn2_2 = min(nn2_2,nene2) 1096 nn2_3 = max(1,nn2_3) 1097 nn2_3 = min(nn2_3,nene2) 1098 1099 inv_deleps1 = 1.0d0 / deleps1 1100 1101 eps2 = enemin2+nn2_1*deltaene2 ! starting value for E 1102 deleps2 = eps2 - eigen2_1tetra(2) ! E - epsilon_1 1103 1104 !----------------------------------------------------------------------- 1105 ! This is case BI 1106 !----------------------------------------------------------------------- 1107 do ieps2 = nn2_1+1, nn2_2 1108 ! calculate running value of del "b" = b-b_s: first term should really mix eps1 and eps2 1109 delbb0 = deleps2*inv_deleps1 - bb(1) 1110 1111 ii1 = dccde1*delbb0*inv_delbb(1)*inv_delbb(2) ! this is I_1(Omega E) 1112 1113 ! note negative sign here - we are correcting the I0 a0 term already calculated above 1114 dtweightde_tmp(2,ieps2,ieps1) = dtweightde_tmp(2,ieps2,ieps1) - & 1115 & ii1*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_b(1),2)* & 1116 & (-2.0d0 + delbb0*inv_delbb(1)*epsilon1(ind_b(2),ind_b(1))*inv_epsilon1(ind_b(2),2) & 1117 & + delbb0*inv_delbb(2)*epsilon1(ind_b(3),ind_b(1))*inv_epsilon1(ind_b(3),2))) 1118 dtweightde_tmp(ind_b(1),ieps2,ieps1) = dtweightde_tmp(ind_b(1),ieps2,ieps1) - & 1119 & ii1*0.5d0*deleps1*inv_epsilon1(ind_b(1),2)*(2.0d0 - delbb0*inv_delbb(1) - delbb0*inv_delbb(2)) 1120 dtweightde_tmp(ind_b(2),ieps2,ieps1) = dtweightde_tmp(ind_b(2),ieps2,ieps1) - & 1121 & ii1*0.5d0*delbb0*inv_delbb(1)*deleps1*inv_epsilon1(ind_b(2),2) 1122 dtweightde_tmp(ind_b(3),ieps2,ieps1) = dtweightde_tmp(ind_b(3),ieps2,ieps1) - & 1123 & ii1*0.5d0*delbb0*inv_delbb(2)*deleps1*inv_epsilon1(ind_b(3),2) 1124 deleps2 = deleps2 + deltaene2 1125 end do 1126 1127 eps2 = enemin2+nn2_2*deltaene2 1128 deleps2 = eps2 - eigen2_1tetra(2) 1129 1130 !----------------------------------------------------------------------- 1131 ! This is case BII 1132 !----------------------------------------------------------------------- 1133 do ieps2 = nn2_2+1, nn2_3 1134 ! calculate running value of del "b" = b_l-b: first term should really mix eps1 and eps2 1135 delbb0 = bb(3) - deleps2*inv_deleps1 1136 1137 ii1 = dccde1*delbb0*inv_delbb(3)*inv_delbb(2) ! this is I_1(Omega E) 1138 1139 ! note negative sign here - we are correcting the I0 a0 term already calculated above 1140 dtweightde_tmp(2,ieps2,ieps1) = dtweightde_tmp(2,ieps2,ieps1) - & 1141 & ii1*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_b(3),2)* & 1142 & (-2.0d0 + delbb0*inv_delbb(3)*epsilon1(ind_b(2),ind_b(3))*inv_epsilon1(ind_b(2),2) & 1143 & + delbb0*inv_delbb(2)*epsilon1(ind_b(1),ind_b(3))*inv_epsilon1(ind_b(1),2))) 1144 dtweightde_tmp(ind_b(3),ieps2,ieps1) = dtweightde_tmp(ind_b(3),ieps2,ieps1) - & 1145 & ii1*0.5d0*deleps1*inv_epsilon1(ind_b(3),2)*(2.0d0 - delbb0*inv_delbb(3) - delbb0*inv_delbb(2)) 1146 dtweightde_tmp(ind_b(2),ieps2,ieps1) = dtweightde_tmp(ind_b(2),ieps2,ieps1) - & 1147 & ii1*0.5d0*delbb0*inv_delbb(3)*deleps1*inv_epsilon1(ind_b(2),2) 1148 dtweightde_tmp(ind_b(1),ieps2,ieps1) = dtweightde_tmp(ind_b(1),ieps2,ieps1) - & 1149 & ii1*0.5d0*delbb0*inv_delbb(2)*deleps1*inv_epsilon1(ind_b(1),2) 1150 1151 deleps2 = deleps2 + deltaene2 1152 end do 1153 1154 deleps1 = deleps1 + deltaene1 1155 end do 1156 1157 ! 1158 ! interval e3 < eps < e4 1159 ! 1160 eps1 = eps1 + (nn1_3-nn1_2)*deltaene1 1161 deleps1 = eps1-eigen1_1tetra(4) 1162 dccde1_pre = 6.d0*volconst_mult*inv_epsilon1(4,1)*inv_epsilon1(4,2)*inv_epsilon1(4,3) 1163 do ieps1=nn1_3+1,nn1_4 1164 ! note - sign from definition of f3 1165 dccde1 = -dccde1_pre * deleps1 ! f3(Omega) * 6 * v 1166 1167 ! at fixed ieps1 we can find the pivot indices for the ieps2 loop 1168 ! NB: order is inverted for cc because deleps1 is defined negative (Omega is always less than omega_3) 1169 nn2_1 = int((eigen2_1tetra(4)+deleps1*cc(3) -enemin2)/deltaene2)+1 1170 nn2_2 = int((eigen2_1tetra(4)+deleps1*cc(2) -enemin2)/deltaene2)+1 1171 nn2_3 = int((eigen2_1tetra(4)+deleps1*cc(1) -enemin2)/deltaene2)+1 1172 1173 nn2_1 = max(1,nn2_1) 1174 nn2_1 = min(nn2_1,nene2) 1175 nn2_2 = max(1,nn2_2) 1176 nn2_2 = min(nn2_2,nene2) 1177 nn2_3 = max(1,nn2_3) 1178 nn2_3 = min(nn2_3,nene2) 1179 inv_deleps1 = 1.0d0 / deleps1 1180 1181 eps2 = enemin2+nn2_1*deltaene2 ! starting value for E 1182 deleps2 = eps2 - eigen2_1tetra(4) ! E - epsilon_3 1183 1184 !----------------------------------------------------------------------- 1185 ! This is case CII 1186 !----------------------------------------------------------------------- 1187 do ieps2 = nn2_1+1, nn2_2 1188 ! calculate running value of del "c" = c_l-c: first term should really mix eps1 and eps2 1189 delcc0 = cc(3) - deleps2*inv_deleps1 1190 1191 ii3 = dccde1*delcc0*inv_delcc(3)*inv_delcc(2) ! this is I_3(Omega E) 1192 1193 dtweightde_tmp(4,ieps2,ieps1) = dtweightde_tmp(4,ieps2,ieps1) + & 1194 & ii3*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_c(3),4)* & 1195 & (-2.0d0 + delcc0*inv_delcc(3)*epsilon1(ind_c(2),ind_c(3))*inv_epsilon1(ind_c(2),4) & 1196 & + delcc0*inv_delcc(2)*epsilon1(ind_c(1),ind_c(3))*inv_epsilon1(ind_c(1),4))) 1197 dtweightde_tmp(ind_c(3),ieps2,ieps1) = dtweightde_tmp(ind_c(3),ieps2,ieps1) + & 1198 & ii3*0.5d0*deleps1*inv_epsilon1(ind_c(3),4)*(2.0d0 - delcc0*inv_delcc(3) - delcc0*inv_delcc(2)) 1199 dtweightde_tmp(ind_c(2),ieps2,ieps1) = dtweightde_tmp(ind_c(2),ieps2,ieps1) + & 1200 & ii3*0.5d0*delcc0*inv_delcc(3)*deleps1*inv_epsilon1(ind_c(2),4) 1201 dtweightde_tmp(ind_c(1),ieps2,ieps1) = dtweightde_tmp(ind_c(1),ieps2,ieps1) + & 1202 & ii3*0.5d0*delcc0*inv_delcc(2)*deleps1*inv_epsilon1(ind_c(1),4) 1203 1204 deleps2 = deleps2 + deltaene2 1205 end do 1206 1207 1208 eps2 = enemin2+nn2_2*deltaene2 1209 deleps2 = eps2 - eigen2_1tetra(4) 1210 1211 !----------------------------------------------------------------------- 1212 ! This is case CI 1213 !----------------------------------------------------------------------- 1214 do ieps2 = nn2_2+1, nn2_3 1215 ! calculate running value of del "c" = c-c_s: first term should really mix eps1 and eps2 1216 delcc0 = deleps2*inv_deleps1 - cc(1) ! c - c_s 1217 1218 ii3 = dccde1*delcc0*inv_delcc(1)*inv_delcc(2) ! this is I_3(Omega E) 1219 1220 dtweightde_tmp(4,ieps2,ieps1) = dtweightde_tmp(4,ieps2,ieps1) + & 1221 & ii3*(1.d0 + 0.5d0*deleps1*inv_epsilon1(ind_c(1),4)* & 1222 & (-2.0d0 + delcc0*inv_delcc(1)*epsilon1(ind_c(2),ind_c(1))*inv_epsilon1(ind_c(2),4) & 1223 & + delcc0*inv_delcc(2)*epsilon1(ind_c(3),ind_c(1))*inv_epsilon1(ind_c(3),4))) 1224 dtweightde_tmp(ind_c(1),ieps2,ieps1) = dtweightde_tmp(ind_c(1),ieps2,ieps1) + & 1225 & ii3*0.5d0*deleps1*inv_epsilon1(ind_c(1),4)*(2.0d0 - delcc0*inv_delcc(1) - delcc0*inv_delcc(2)) 1226 dtweightde_tmp(ind_c(2),ieps2,ieps1) = dtweightde_tmp(ind_c(2),ieps2,ieps1) + & 1227 & ii3*0.5d0*delcc0*inv_delcc(1)*deleps1*inv_epsilon1(ind_c(2),4) 1228 dtweightde_tmp(ind_c(3),ieps2,ieps1) = dtweightde_tmp(ind_c(3),ieps2,ieps1) + & 1229 & ii3*0.5d0*delcc0*inv_delcc(2)*deleps1*inv_epsilon1(ind_c(3),4) 1230 deleps2 = deleps2 + deltaene2 1231 end do 1232 1233 deleps1 = deleps1 + deltaene1 1234 end do 1235 1236 eps1 = eps1 + (nn1_4-nn1_3)*deltaene1 1237 ! 1238 ! 1239 ! interval e4 < eps < enemax 1240 ! 1241 do ieps1=nn1_4+1,nene1 1242 ! dtweightde unchanged by this tetrahedron 1243 end do 1244 1245 ! if we have a fully degenerate tetrahedron, 1246 ! 1) the tweight is a Heaviside (step) function, which is correct above, but 1247 ! 2) the dtweightde should contain a Dirac function: add a Gaussian here 1248 1249 ! TODO: add treatment in double tetra case 1250 ! end degenerate tetrahedron if 1251 1252 ! the following blas calls are not working systematically, or do not give speed ups, strange... 1253 !call daxpy (nene, 1.d0, dtweightde_tmp(:,1), 1, dtweightde_t(:,ind_ibz(1)), 1) 1254 !call daxpy (nene, 1.d0, dtweightde_tmp(:,2), 1, dtweightde_t(:,ind_ibz(2)), 1) 1255 !call daxpy (nene, 1.d0, dtweightde_tmp(:,3), 1, dtweightde_t(:,ind_ibz(3)), 1) 1256 !call daxpy (nene, 1.d0, dtweightde_tmp(:,4), 1, dtweightde_t(:,ind_ibz(4)), 1) 1257 1258 do ieps2 = 1, nene2 1259 dtweightde(ind_k(1),:,ieps2) = dtweightde(ind_k(1),:,ieps2) + dtweightde_tmp(1,ieps2,:) 1260 dtweightde(ind_k(2),:,ieps2) = dtweightde(ind_k(2),:,ieps2) + dtweightde_tmp(2,ieps2,:) 1261 dtweightde(ind_k(3),:,ieps2) = dtweightde(ind_k(3),:,ieps2) + dtweightde_tmp(3,ieps2,:) 1262 dtweightde(ind_k(4),:,ieps2) = dtweightde(ind_k(4),:,ieps2) + dtweightde_tmp(4,ieps2,:) 1263 !tweight(nkpt,nene1,nene2) 1264 end do 1265 1266 end do ! itetra 1267 1268 ! transpose: otherwise the data access is crap and the code slows by an order of magnitude 1269 TETRA_DEALLOCATE(tweight_tmp) 1270 TETRA_DEALLOCATE(dtweightde_tmp) 1271 1272 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
SOURCE
1472 pure subroutine get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,nene,bcorr, & 1473 & tweight_tmp,dtweightde_tmp) 1474 1475 !Arguments ------------------------------------ 1476 !scalars 1477 integer,intent(in) :: nene,bcorr,itetra 1478 type(t_tetrahedron), intent(in) :: tetra 1479 real(dp) ,intent(in) :: enemax,enemin,max_occ 1480 !arrays 1481 ! MGTODO: This layout is not optimal (lots of cache thrashing, I will optimize it later on) 1482 real(dp), intent(out) :: tweight_tmp(nene, 4) 1483 real(dp), intent(out) :: dtweightde_tmp(nene, 4) 1484 real(dp),intent(in) :: eigen_1tetra(4) 1485 1486 !Local variables------------------------------- 1487 ! needed for gaussian replacement of Dirac functions 1488 ! the three coefficients of the DOS as quadratic form, 1489 ! in the interval [eig(ikpt-1), eig(ikpt)] 1490 ! for ikpt = 1 we add a point below eigen(1) which doesnt 1491 ! contribute to the DOS in any tetrahedron 1492 !scalars 1493 integer :: ieps,nn1,nn2,nn3,nn4 1494 real(dp) :: cc,cc1,cc2,cc3,dcc1de,dcc2de,dcc3de,dccde,deltaene,eps 1495 real(dp) :: epsilon21,epsilon31,epsilon32,epsilon41,epsilon42,epsilon43 1496 real(dp) :: gau_prefactor,gau_width,gau_width2,inv_epsilon21,inv_epsilon31,gval 1497 real(dp) :: inv_epsilon32,inv_epsilon41,inv_epsilon42,inv_epsilon43 1498 real(dp) :: deleps1, deleps2, deleps3, deleps4 1499 real(dp) :: invepsum, cc_pre, dccde_pre 1500 real(dp) :: cc1_pre, cc2_pre, cc3_pre 1501 real(dp) :: cc_tmp, dccde_tmp 1502 real(dp) :: dcc1de_pre, dcc2de_pre, dcc3de_pre 1503 real(dp) :: tmp,volconst,volconst_mult 1504 1505 ! ********************************************************************* 1506 1507 volconst = tetra%vv/4.d0 1508 1509 deltaene = (enemax-enemin) / (nene-1) 1510 1511 ! This is output 1512 tweight_tmp = zero; dtweightde_tmp = zero 1513 1514 volconst_mult = max_occ*volconst*dble(tetra%tetra_mult(itetra)) 1515 1516 ! all notations are from Blochl PRB 49 16223 [[cite:Bloechl1994a]] Appendix B 1517 epsilon21 = eigen_1tetra(2)-eigen_1tetra(1) 1518 epsilon31 = eigen_1tetra(3)-eigen_1tetra(1) 1519 epsilon41 = eigen_1tetra(4)-eigen_1tetra(1) 1520 epsilon32 = eigen_1tetra(3)-eigen_1tetra(2) 1521 epsilon42 = eigen_1tetra(4)-eigen_1tetra(2) 1522 epsilon43 = eigen_1tetra(4)-eigen_1tetra(3) 1523 inv_epsilon21 = zero; if (epsilon21 > tol6) inv_epsilon21 = 1.d0 / epsilon21 1524 inv_epsilon31 = zero; if (epsilon31 > tol6) inv_epsilon31 = 1.d0 / epsilon31 1525 inv_epsilon41 = zero; if (epsilon41 > tol6) inv_epsilon41 = 1.d0 / epsilon41 1526 inv_epsilon32 = zero; if (epsilon32 > tol6) inv_epsilon32 = 1.d0 / epsilon32 1527 inv_epsilon42 = zero; if (epsilon42 > tol6) inv_epsilon42 = 1.d0 / epsilon42 1528 inv_epsilon43 = zero; if (epsilon43 > tol6) inv_epsilon43 = 1.d0 / epsilon43 1529 1530 nn1 = int((eigen_1tetra(1)-enemin)/deltaene)+1 1531 nn2 = int((eigen_1tetra(2)-enemin)/deltaene)+1 1532 nn3 = int((eigen_1tetra(3)-enemin)/deltaene)+1 1533 nn4 = int((eigen_1tetra(4)-enemin)/deltaene)+1 1534 1535 nn1 = max(1,nn1) 1536 nn1 = min(nn1,nene) 1537 nn2 = max(1,nn2) 1538 nn2 = min(nn2,nene) 1539 nn3 = max(1,nn3) 1540 nn3 = min(nn3,nene) 1541 nn4 = max(1,nn4) 1542 nn4 = min(nn4,nene) 1543 1544 eps = enemin+nn1*deltaene 1545 ! 1546 !interval enemin < eps < e1 nothing to do 1547 ! 1548 ! 1549 !interval e1 < eps < e2 1550 ! 1551 deleps1 = eps-eigen_1tetra(1) 1552 cc_pre = volconst_mult*inv_epsilon21*inv_epsilon31*inv_epsilon41 1553 invepsum = inv_epsilon21+inv_epsilon31+inv_epsilon41 1554 dccde_pre = 3.d0*volconst_mult*inv_epsilon21*inv_epsilon31*inv_epsilon41 1555 do ieps=nn1+1,nn2 1556 cc = cc_pre * deleps1*deleps1*deleps1 1557 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + cc*(4.d0-deleps1*invepsum) 1558 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + cc*deleps1*inv_epsilon21 1559 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + cc*deleps1*inv_epsilon31 1560 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + cc*deleps1*inv_epsilon41 1561 1562 dccde = dccde_pre * deleps1*deleps1 1563 dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + dccde*(4.d0 - deleps1*invepsum) -cc*invepsum 1564 dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + (dccde*deleps1 + cc) * inv_epsilon21 1565 dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + (dccde*deleps1 + cc) * inv_epsilon31 1566 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + (dccde*deleps1 + cc) * inv_epsilon41 1567 1568 if (bcorr == 1) then 1569 ! bxu, correction terms based on Bloechl's paper 1570 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + & 1571 & 4.d0*dccde_pre*deleps1*deleps1*(epsilon21+epsilon31+epsilon41)/40.d0 1572 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + & 1573 & 4.d0*dccde_pre*deleps1*deleps1*(-epsilon21+epsilon32+epsilon42)/40.d0 1574 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + & 1575 & 4.d0*dccde_pre*deleps1*deleps1*(-epsilon31-epsilon32+epsilon43)/40.d0 1576 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + & 1577 & 4.d0*dccde_pre*deleps1*deleps1*(-epsilon41-epsilon42-epsilon43)/40.d0 1578 1579 dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + & 1580 & 8.d0*dccde_pre*deleps1*(epsilon21+epsilon31+epsilon41)/40.d0 1581 dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + & 1582 & 8.d0*dccde_pre*deleps1*(-epsilon21+epsilon32+epsilon42)/40.d0 1583 dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + & 1584 & 8.d0*dccde_pre*deleps1*(-epsilon31-epsilon32+epsilon43)/40.d0 1585 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + & 1586 & 8.d0*dccde_pre*deleps1*(-epsilon41-epsilon42-epsilon43)/40.d0 1587 end if 1588 1589 deleps1 = deleps1 + deltaene 1590 end do 1591 1592 eps = eps + (nn2-nn1)*deltaene 1593 ! 1594 ! interval e2 < eps < e3 1595 ! 1596 deleps1 = eps-eigen_1tetra(1) 1597 deleps2 = eps-eigen_1tetra(2) 1598 deleps3 = eigen_1tetra(3)-eps 1599 deleps4 = eigen_1tetra(4)-eps 1600 1601 cc1_pre = volconst_mult*inv_epsilon31*inv_epsilon41 1602 cc2_pre = volconst_mult*inv_epsilon41*inv_epsilon32*inv_epsilon31 1603 cc3_pre = volconst_mult*inv_epsilon42*inv_epsilon32*inv_epsilon41 1604 1605 dcc1de_pre = 2.d0*cc1_pre 1606 dcc2de_pre = cc2_pre 1607 dcc3de_pre = cc3_pre 1608 do ieps=nn2+1,nn3 1609 cc1 = cc1_pre * deleps1*deleps1 1610 cc2 = cc2_pre * deleps1*deleps2*deleps3 1611 cc3 = cc3_pre * deleps2*deleps2*deleps4 1612 1613 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + & 1614 & cc1 + (cc1+cc2)*deleps3*inv_epsilon31 + (cc1+cc2+cc3)*deleps4*inv_epsilon41 1615 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + & 1616 & cc1+cc2+cc3+(cc2+cc3)*deleps3*inv_epsilon32 + cc3*deleps4*inv_epsilon42 1617 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + & 1618 & (cc1+cc2)*deleps1*inv_epsilon31 + (cc2+cc3)*deleps2*inv_epsilon32 1619 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + & 1620 & (cc1+cc2+cc3)*deleps1*inv_epsilon41 + cc3*deleps2*inv_epsilon42 1621 1622 1623 dcc1de = dcc1de_pre * deleps1 1624 dcc2de = dcc2de_pre * (-deleps1*deleps2 +deleps1*deleps3 +deleps2*deleps3) 1625 dcc3de = dcc3de_pre * (2.d0*deleps2*deleps4 -deleps2*deleps2) 1626 1627 dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) & 1628 & + dcc1de & 1629 & + ((dcc1de+dcc2de)*deleps3 -(cc1+cc2)) * inv_epsilon31 & 1630 & + ((dcc1de+dcc2de+dcc3de)*deleps4 -(cc1+cc2+cc3)) * inv_epsilon41 1631 dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) & 1632 & + dcc1de+dcc2de+dcc3de & 1633 & + ((dcc2de+dcc3de)*deleps3 -(cc2+cc3) ) * inv_epsilon32 & 1634 & + (dcc3de*deleps4 -cc3 ) * inv_epsilon42 1635 dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) & 1636 & + ((dcc1de+dcc2de)*deleps1 + (cc1+cc2) ) * inv_epsilon31 & 1637 & + ((dcc2de+dcc3de)*deleps2 + (cc2+cc3) ) * inv_epsilon32 1638 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) & 1639 & + ((dcc1de+dcc2de+dcc3de)*deleps1 + (cc1+cc2+cc3) ) * inv_epsilon41 & 1640 & + (dcc3de*deleps2 + cc3) * inv_epsilon42 1641 1642 if (bcorr == 1) then 1643 ! bxu, correction terms based on Bloechl's paper 1644 ! The correction terms may cause the dtweightde become negative 1645 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + & 1646 & 4.d0*cc1_pre* & 1647 & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* & 1648 & (epsilon21+epsilon31+epsilon41)/40.d0 1649 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + & 1650 & 4.d0*cc1_pre* & 1651 & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* & 1652 & (-epsilon21+epsilon32+epsilon42)/40.d0 1653 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + & 1654 & 4.d0*cc1_pre* & 1655 & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* & 1656 & (-epsilon31-epsilon32+epsilon43)/40.d0 1657 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + & 1658 & 4.d0*cc1_pre* & 1659 & (3.d0*epsilon21+6.d0*deleps2-3.d0*(epsilon31+epsilon42)*deleps2**2.d0*inv_epsilon32*inv_epsilon42)* & 1660 & (-epsilon41-epsilon42-epsilon43)/40.d0 1661 1662 dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + & 1663 & 4.d0*cc1_pre* & 1664 & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* & 1665 & (epsilon21+epsilon31+epsilon41)/40.d0 1666 dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + & 1667 & 4.d0*cc1_pre* & 1668 & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* & 1669 & (-epsilon21+epsilon32+epsilon42)/40.d0 1670 dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + & 1671 & 4.d0*cc1_pre* & 1672 & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* & 1673 & (-epsilon31-epsilon32+epsilon43)/40.d0 1674 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + & 1675 & 4.d0*cc1_pre* & 1676 & (6.d0-6.d0*(epsilon31+epsilon42)*deleps2*inv_epsilon32*inv_epsilon42)* & 1677 & (-epsilon41-epsilon42-epsilon43)/40.d0 1678 end if 1679 1680 deleps1 = deleps1 + deltaene 1681 deleps2 = deleps2 + deltaene 1682 deleps3 = deleps3 - deltaene 1683 deleps4 = deleps4 - deltaene 1684 end do 1685 1686 eps = eps + (nn3-nn2)*deltaene 1687 ! 1688 ! interval e3 < eps < e4 1689 ! 1690 deleps4 = eigen_1tetra(4)-eps 1691 cc_pre = volconst_mult*inv_epsilon41*inv_epsilon42*inv_epsilon43 1692 invepsum = inv_epsilon41+inv_epsilon42+inv_epsilon43 1693 dccde_pre = -3.d0*cc_pre 1694 do ieps=nn3+1,nn4 1695 cc = cc_pre * deleps4*deleps4*deleps4 1696 cc_tmp = cc * deleps4 1697 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + volconst_mult - cc_tmp*inv_epsilon41 1698 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + volconst_mult - cc_tmp*inv_epsilon42 1699 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + volconst_mult - cc_tmp*inv_epsilon43 1700 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + volconst_mult - cc*4.d0 + cc_tmp*invepsum 1701 1702 dccde = dccde_pre * deleps4*deleps4 1703 dccde_tmp = -dccde*deleps4 + cc 1704 dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + dccde_tmp * inv_epsilon41 1705 dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + dccde_tmp * inv_epsilon42 1706 dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + dccde_tmp * inv_epsilon43 1707 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) - dccde*4.d0 - dccde_tmp*invepsum 1708 1709 if (bcorr == 1) then 1710 ! bxu, correction terms based on Bloechl's paper 1711 ! The correction terms may cause the dtweightde become negative 1712 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + & 1713 & 12.d0*cc_pre*deleps4*deleps4*(epsilon21+epsilon31+epsilon41)/40.d0 1714 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + & 1715 & 12.d0*cc_pre*deleps4*deleps4*(-epsilon21+epsilon32+epsilon42)/40.d0 1716 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + & 1717 & 12.d0*cc_pre*deleps4*deleps4*(-epsilon31-epsilon32+epsilon43)/40.d0 1718 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + & 1719 & 12.d0*cc_pre*deleps4*deleps4*(-epsilon41-epsilon42-epsilon43)/40.d0 1720 1721 dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) - & 1722 & 24.d0*cc_pre*deleps4*(epsilon21+epsilon31+epsilon41)/40.d0 1723 dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) - & 1724 & 24.d0*cc_pre*deleps4*(-epsilon21+epsilon32+epsilon42)/40.d0 1725 dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) - & 1726 & 24.d0*cc_pre*deleps4*(-epsilon31-epsilon32+epsilon43)/40.d0 1727 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) - & 1728 & 24.d0*cc_pre*deleps4*(-epsilon41-epsilon42-epsilon43)/40.d0 1729 end if 1730 1731 deleps4 = deleps4 - deltaene 1732 end do 1733 eps = eps + (nn4-nn3)*deltaene 1734 ! 1735 ! 1736 ! interval e4 < eps < enemax 1737 ! 1738 do ieps=nn4+1,nene 1739 tweight_tmp(ieps,1) = tweight_tmp(ieps,1) + volconst_mult 1740 tweight_tmp(ieps,2) = tweight_tmp(ieps,2) + volconst_mult 1741 tweight_tmp(ieps,3) = tweight_tmp(ieps,3) + volconst_mult 1742 tweight_tmp(ieps,4) = tweight_tmp(ieps,4) + volconst_mult 1743 ! dtweightde unchanged by this tetrahedron 1744 end do 1745 1746 ! 1747 ! if we have a fully degenerate tetrahedron, 1748 ! 1) the tweight is a Heaviside (step) function, which is correct above, but 1749 ! 2) the dtweightde should contain a Dirac function: add a Gaussian here 1750 ! 1751 if (epsilon41 < tol6) then 1752 1753 ! to ensure the gaussian will integrate properly: 1754 ! WARNING: this smearing could be problematic if too large 1755 ! and doesnt integrate well if its too small 1756 gau_width = 10.0d0*deltaene 1757 gau_width2 = 1.0 / gau_width / gau_width 1758 gau_prefactor = volconst_mult / gau_width / sqrtpi 1759 ! 1760 ! average position since bracket for epsilon41 is relatively large 1761 cc = (eigen_1tetra(1)+eigen_1tetra(2)+eigen_1tetra(3)+eigen_1tetra(4))/4.d0 1762 eps = enemin 1763 do ieps=1,nene 1764 tmp = eps - cc 1765 gval = gau_prefactor*exp(-tmp*tmp*gau_width2) 1766 ! MG TODO: I think this is not correct, because we have divided by 4 so 1767 ! the other points should be accumulated as well. 1768 ! There are however changes in the unit tests if I activate these lines... 1769 !dtweightde_tmp(ieps,1) = dtweightde_tmp(ieps,1) + gval 1770 !dtweightde_tmp(ieps,2) = dtweightde_tmp(ieps,2) + gval 1771 !dtweightde_tmp(ieps,3) = dtweightde_tmp(ieps,3) + gval 1772 dtweightde_tmp(ieps,4) = dtweightde_tmp(ieps,4) + gval 1773 eps = eps + deltaene 1774 end do 1775 end if ! end degenerate tetrahedron if 1776 1777 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
SOURCE
642 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 643 ! THIS FUNCTION IS DEPRECATED, USE tetra_blochl_weights 644 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 645 subroutine get_tetra_weight(eigen_in,enemin,enemax,max_occ,nene,nkpt,tetra,& 646 bcorr,tweight,dtweightde,comm) 647 648 !Arguments ------------------------------------ 649 !scalars 650 integer,intent(in) :: nene,nkpt,bcorr,comm 651 type(t_tetrahedron), intent(in) :: tetra 652 real(dp) ,intent(in) :: enemax,enemin,max_occ 653 !arrays 654 real(dp) ,intent(in) :: eigen_in(nkpt) 655 real(dp) ,intent(out) :: dtweightde(nkpt,nene),tweight(nkpt,nene) 656 657 !Local variables------------------------------- 658 real(dp), allocatable :: dtweightde_ek(:, :), tweight_ek(:, :) 659 660 ! ********************************************************************* 661 662 TETRA_ALLOCATE(dtweightde_ek, (nene, nkpt)) 663 TETRA_ALLOCATE(tweight_ek, (nene, nkpt)) 664 665 call tetra_blochl_weights(tetra,eigen_in,enemin,enemax,max_occ,nene,nkpt,bcorr,tweight_ek,dtweightde_ek,comm) 666 667 ! transpose: otherwise the data access is crap and the code slows by an order of magnitude 668 tweight = transpose(tweight_ek) 669 dtweightde = transpose(dtweightde_ek) 670 671 TETRA_DEALLOCATE(dtweightde_ek) 672 TETRA_DEALLOCATE(tweight_ek) 673 674 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 comm= MPI communicator
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
SOURCE
184 subroutine init_tetra(indkpt, gprimd, klatt, kpt_fullbz, nkpt_fullbz, tetra, ierr, errorstring, comm) 185 186 !Arguments ------------------------------------ 187 !scalars 188 integer,intent(in) :: nkpt_fullbz, comm 189 integer, intent(out) :: ierr 190 character(len=80), intent(out) :: errorstring 191 type(t_tetrahedron),intent(out) :: tetra 192 !arrays 193 integer,intent(in) :: indkpt(nkpt_fullbz) 194 real(dp) ,intent(in) :: gprimd(3,3),klatt(3,3),kpt_fullbz(3,nkpt_fullbz) 195 196 !Local variables------------------------------- 197 !scalars 198 integer :: ialltetra,ikpt2,ikpt_full,isummit,itetra,jalltetra,jsummit 199 integer :: ii,jj,ikibz,nkpt_ibz, my_rank, nprocs 200 integer :: symrankkpt,mtetra,itmp,ntetra_irred 201 real(dp) :: shift1,shift2,shift3, rcvol,hashfactor 202 !real :: cpu_start, cpu_stop 203 type(krank_t) :: krank 204 !arrays 205 integer :: ind_ibz(4), tetra_shifts(3,4,6) ! 3 dimensions, 4 summits, and 6 tetrahedra / kpoint box 206 real(dp) :: k1(3),k2(3),k3(3) 207 integer,allocatable :: tetra_full_(:,:,:) 208 integer,allocatable :: tetra_mult_(:) 209 integer,allocatable :: tetra_wrap_(:,:,:) 210 integer, allocatable :: reforder(:) 211 integer, allocatable :: irred_itetra(:) 212 real(dp), allocatable :: tetra_hash(:) 213 214 ! ********************************************************************* 215 216 !call cpu_time(cpu_start) 217 218 my_rank = 0; nprocs = 1 219 #ifdef HAVE_LIBTETRA_ABINIT 220 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 221 #endif 222 223 ierr = 0 224 errorstring = "" 225 !jmb 226 shift1 = zero 227 shift2 = zero 228 shift3 = zero 229 230 tetra%klatt = klatt 231 232 mtetra = 6 * nkpt_fullbz 233 TETRA_ALLOCATE(tetra_full_, (4,2,mtetra)) 234 TETRA_ALLOCATE(tetra_mult_, (mtetra)) 235 TETRA_ALLOCATE(tetra_wrap_, (3,4,mtetra)) 236 237 tetra_mult_ = 1 238 tetra_full_ = 0 239 tetra_wrap_ = 0 240 241 ! tetra_shifts(:,1,1) = (/0,0,0/) 242 ! tetra_shifts(:,2,1) = (/0,1,0/) 243 ! tetra_shifts(:,3,1) = (/0,1,1/) 244 ! tetra_shifts(:,4,1) = (/1,1,0/) 245 ! tetra_shifts(:,1,2) = (/0,0,0/) 246 ! tetra_shifts(:,2,2) = (/0,1,1/) 247 ! tetra_shifts(:,3,2) = (/1,1,0/) 248 ! tetra_shifts(:,4,2) = (/1,1,1/) 249 ! tetra_shifts(:,1,3) = (/0,0,0/) 250 ! tetra_shifts(:,2,3) = (/1,0,0/) 251 ! tetra_shifts(:,3,3) = (/1,1,0/) 252 ! tetra_shifts(:,4,3) = (/1,1,1/) 253 ! tetra_shifts(:,1,4) = (/0,0,0/) 254 ! tetra_shifts(:,2,4) = (/0,0,1/) 255 ! tetra_shifts(:,3,4) = (/1,0,0/) 256 ! tetra_shifts(:,4,4) = (/1,1,1/) 257 ! tetra_shifts(:,1,5) = (/0,0,1/) 258 ! tetra_shifts(:,2,5) = (/1,0,0/) 259 ! tetra_shifts(:,3,5) = (/1,0,1/) 260 ! tetra_shifts(:,4,5) = (/1,1,1/) 261 ! tetra_shifts(:,1,6) = (/0,0,0/) 262 ! tetra_shifts(:,2,6) = (/0,0,1/) 263 ! tetra_shifts(:,3,6) = (/0,1,1/) 264 ! tetra_shifts(:,4,6) = (/1,1,1/) 265 266 ! bxu, the following division scheme is according to Bloechl's paper 267 tetra_shifts(:,1,1) = (/0,0,0/) 268 tetra_shifts(:,2,1) = (/1,0,0/) 269 tetra_shifts(:,3,1) = (/0,1,0/) 270 tetra_shifts(:,4,1) = (/1,0,1/) 271 tetra_shifts(:,1,2) = (/1,0,0/) 272 tetra_shifts(:,2,2) = (/1,1,0/) 273 tetra_shifts(:,3,2) = (/0,1,0/) 274 tetra_shifts(:,4,2) = (/1,0,1/) 275 tetra_shifts(:,1,3) = (/0,1,0/) 276 tetra_shifts(:,2,3) = (/1,1,0/) 277 tetra_shifts(:,3,3) = (/1,0,1/) 278 tetra_shifts(:,4,3) = (/1,1,1/) 279 tetra_shifts(:,1,4) = (/0,0,0/) 280 tetra_shifts(:,2,4) = (/0,1,0/) 281 tetra_shifts(:,3,4) = (/0,0,1/) 282 tetra_shifts(:,4,4) = (/1,0,1/) 283 tetra_shifts(:,1,5) = (/0,0,1/) 284 tetra_shifts(:,2,5) = (/1,0,1/) 285 tetra_shifts(:,3,5) = (/0,1,0/) 286 tetra_shifts(:,4,5) = (/0,1,1/) 287 tetra_shifts(:,1,6) = (/0,1,0/) 288 tetra_shifts(:,2,6) = (/1,0,1/) 289 tetra_shifts(:,3,6) = (/0,1,1/) 290 tetra_shifts(:,4,6) = (/1,1,1/) 291 292 ! Make full k-point rank arrays 293 ! TODO: Lot of memory allocated here if dense mesh e.g ~ 300 ** 3 294 krank = krank_new(nkpt_fullbz, kpt_fullbz) 295 296 ialltetra = 1 297 do ikpt_full=1,nkpt_fullbz 298 do itetra=1,6 299 !ialltetra = itetra + (ikpt_full -1) * 6 300 !if (mod(ialltetra, nprocs) /= my_rank) cycle ! MPI parallelism. 301 do isummit=1,4 302 k1(:) = kpt_fullbz(:,ikpt_full) & 303 + tetra_shifts(1,isummit,itetra)*klatt(:,1) & 304 + tetra_shifts(2,isummit,itetra)*klatt(:,2) & 305 + tetra_shifts(3,isummit,itetra)*klatt(:,3) 306 307 ! Find full kpoint which is summit isummit of tetrahedron itetra around full kpt ikpt_full ! 308 symrankkpt = krank%get_rank(k1) 309 ikpt2 = krank%invrank(symrankkpt) 310 if (ikpt2 < 1) then 311 errorstring='Error in ranking k-points - exiting with un-initialized tetrahedra.' 312 ierr = 2 313 call krank%free() 314 TETRA_ALLOCATE(tetra%tetra_full, (4,2,1)) 315 TETRA_ALLOCATE(tetra%tetra_mult, (1)) 316 TETRA_ALLOCATE(tetra%tetra_wrap, (3,4,1)) 317 TETRA_DEALLOCATE(tetra_full_) 318 TETRA_DEALLOCATE(tetra_mult_) 319 TETRA_DEALLOCATE(tetra_wrap_) 320 return 321 end if 322 323 ! Store irreducible kpoint equivalent to kpt_fullbz(:,ikpt2) 324 tetra_full_(isummit,1,ialltetra) = indkpt(ikpt2) 325 tetra_full_(isummit,2,ialltetra) = ikpt2 326 shift1 = k1(1)-kpt_fullbz(1,ikpt2) 327 shift2 = k1(2)-kpt_fullbz(2,ikpt2) 328 shift3 = k1(3)-kpt_fullbz(3,ikpt2) 329 if (shift1>0.5d0) then 330 tetra_wrap_(1,isummit,ialltetra) = 1 331 else if (shift1<-0.5d0) then 332 tetra_wrap_(1,isummit,ialltetra) = -1 333 end if 334 if (shift2>0.5d0) then 335 tetra_wrap_(2,isummit,ialltetra) = 1 336 else if (shift2<-0.5d0) then 337 tetra_wrap_(2,isummit,ialltetra) = -1 338 end if 339 if (shift3>0.5d0) then 340 tetra_wrap_(3,isummit,ialltetra) = 1 341 else if (shift3<-0.5d0) then 342 tetra_wrap_(3,isummit,ialltetra) = -1 343 end if 344 345 ! sort itetra summits 346 ! TODO: replace with sort_int 347 do jsummit=isummit,2,-1 348 if ( tetra_full_(jsummit,1,ialltetra) < tetra_full_(jsummit-1,1,ialltetra) ) then 349 itmp = tetra_full_(jsummit,1,ialltetra) 350 tetra_full_(jsummit,1,ialltetra) = tetra_full_(jsummit-1,1,ialltetra) 351 tetra_full_(jsummit-1,1,ialltetra) = itmp 352 itmp = tetra_full_(jsummit,2,ialltetra) 353 tetra_full_(jsummit,2,ialltetra) = tetra_full_(jsummit-1,2,ialltetra) 354 tetra_full_(jsummit-1,2,ialltetra) = itmp 355 ! keep fullbz_kpt tetrahedra points in same order 356 itmp = tetra_wrap_(1,jsummit,ialltetra) 357 tetra_wrap_(1,jsummit,ialltetra) = tetra_wrap_(1,jsummit-1,ialltetra) 358 tetra_wrap_(1,jsummit-1,ialltetra) = itmp 359 itmp = tetra_wrap_(2,jsummit,ialltetra) 360 tetra_wrap_(2,jsummit,ialltetra) = tetra_wrap_(2,jsummit-1,ialltetra) 361 tetra_wrap_(2,jsummit-1,ialltetra) = itmp 362 itmp = tetra_wrap_(1,jsummit,ialltetra) 363 tetra_wrap_(3,jsummit,ialltetra) = tetra_wrap_(3,jsummit-1,ialltetra) 364 tetra_wrap_(3,jsummit-1,ialltetra) = itmp 365 end if 366 end do ! jsummit 367 368 end do ! isummit 369 370 if (ialltetra > mtetra) then 371 write (errorstring, '(3a,i0,a,i0)' ) & 372 'init_tetra: BUG - ',& 373 ' ialltetra > mtetra ',& 374 ' ialltetra= ',ialltetra,', mtetra= ',mtetra 375 ierr = 1 376 return 377 end if 378 ialltetra = ialltetra+1 379 end do ! itetra 380 end do ! ikpt_full 381 382 !call cpu_time(cpu_stop) 383 !write(*,*)"tetra_init ikpt_loop:", cpu_stop - cpu_start 384 !cpu_start = cpu_stop 385 386 call krank%free() 387 388 rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) & 389 & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) & 390 & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3))) 391 392 ! Volume of all tetrahedra should be the same as that of tetra 1 393 ! this is the volume of 1 tetrahedron, should be coherent with notation in Lehmann & Taut 394 k1(:) = gprimd(:,1)*klatt(1,1) + gprimd(:,2)*klatt(2,1) + gprimd(:,3)*klatt(3,1) 395 k2(:) = gprimd(:,1)*klatt(1,2) + gprimd(:,2)*klatt(2,2) + gprimd(:,3)*klatt(3,2) 396 k3(:) = gprimd(:,1)*klatt(1,3) + gprimd(:,2)*klatt(2,3) + gprimd(:,3)*klatt(3,3) 397 tetra%vv = abs (k1(1)*(k2(2)*k3(3)-k2(3)*k3(2)) & 398 & -k1(2)*(k2(1)*k3(3)-k2(3)*k3(1)) & 399 & +k1(3)*(k2(1)*k3(2)-k2(2)*k3(1))) / 6.d0 / rcvol 400 401 ! eliminate equivalent tetrahedra by symmetry and account for them in multiplicity tetra_mult 402 tetra%ntetra = mtetra 403 404 ! FIXME: could we replace this with a ranking algorithm to avoid the O(tetra%ntetra^2) step? For example: 405 ! get tetrahedron rank - problem too many combinations in principle = nkpt_irred^4 - only a few used in practice 406 ! sort ranks and keep indices 407 408 ! make hash table = tetra_full_(1)*nkptirred**3+tetra_full_(2)*nkptirred**2+tetra_full_(3)*nkptirred**1+tetra_full_(4) 409 410 hashfactor = 100.d0 ! *acos(-1.d0) ! 100 pi should be far from an integer... 411 TETRA_ALLOCATE(tetra_hash, (tetra%ntetra)) 412 TETRA_ALLOCATE(reforder, (tetra%ntetra)) 413 414 !MG: In principle the order of the indices should not matter. 415 do ialltetra=1, tetra%ntetra 416 tetra_hash(ialltetra) = tetra_full_(1,1,ialltetra)*hashfactor**3+& 417 & tetra_full_(2,1,ialltetra)*hashfactor**2+& 418 & tetra_full_(3,1,ialltetra)*hashfactor**1+& 419 & tetra_full_(4,1,ialltetra) 420 reforder(ialltetra) = ialltetra 421 end do 422 423 call sort_tetra(tetra%ntetra, tetra_hash, reforder, tol6) 424 ! Most of the wall-time is spent in the preamble of this routine (up to this point). 425 ! sort_tetra is not easy to parallelize... 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 ! Create mapping between the irreducible k-points 482 ! and all the tetrahedron contributing with some weight 483 nkpt_ibz = maxval(indkpt) 484 485 ! 1. First we count what is the maximum number of distinct tetrahedra that each k-point contains 486 TETRA_ALLOCATE(tetra%ibz_tetra_count,(nkpt_ibz)) 487 tetra%ibz_tetra_count(:) = 0 488 489 ! Count max tetra contributing 490 do ii=1,tetra%ntetra 491 ! Here we need the original ordering to reference the correct irred kpoints 492 ind_ibz(:) = tetra%tetra_full(:,1,ii) 493 ! count max tetra contributing 494 do jj=1,4 495 ikibz = ind_ibz(jj) 496 if (ikibz > nkpt_ibz) cycle 497 tetra%ibz_tetra_count(ikibz) = tetra%ibz_tetra_count(ikibz) + 1 498 end do 499 end do 500 501 ! 2. Then we build mapping of ikbz to tetra 502 TETRA_ALLOCATE(tetra%ibz_tetra_mapping,(nkpt_ibz,maxval(tetra%ibz_tetra_count))) 503 tetra%ibz_tetra_count(:) = 0 504 do ii=1,tetra%ntetra 505 ! Here we need the original ordering to reference the correct irred kpoints 506 ind_ibz(:) = tetra%tetra_full(:,1,ii) 507 ! Use the counter to move pointer and then fill index 508 do jj=1,4 509 ikibz = ind_ibz(jj) 510 if (ikibz > nkpt_ibz) cycle 511 ! avoid putting the same index twice 512 if (tetra%ibz_tetra_count(ikibz) > 0) then 513 if (tetra%ibz_tetra_mapping(ikibz,tetra%ibz_tetra_count(ikibz)) == ii) cycle 514 end if 515 tetra%ibz_tetra_count(ikibz) = tetra%ibz_tetra_count(ikibz) + 1 516 tetra%ibz_tetra_mapping(ikibz,tetra%ibz_tetra_count(ikibz)) = ii 517 end do 518 end do 519 520 !call cpu_time(cpu_stop) 521 !write(*,*)"tetra_init 2nd part:", cpu_stop - cpu_start 522 !cpu_start = cpu_stop 523 524 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
SOURCE
1298 subroutine sort_tetra(n,list,iperm,tol) 1299 1300 integer, intent(in) :: n 1301 integer, intent(inout) :: iperm(n) 1302 real(dp), intent(inout) :: list(n) 1303 real(dp), intent(in) :: tol 1304 1305 integer :: l,ir,iap,i,j 1306 real(dp) :: ap 1307 character(len=500) :: msg 1308 1309 if (n==1) then 1310 ! Accomodate case of array of length 1: already sorted! 1311 return 1312 else if (n<1) then 1313 ! Should not call with n<1 1314 write(msg,1000) n 1315 1000 format(/,' sort_tetra has been called with array length n=',i12,/, & 1316 & ' having a value less than 1. This is not allowed.') 1317 TETRA_ERROR(msg) 1318 1319 else ! n>1 1320 1321 ! Conduct the usual sort 1322 l=n/2+1 1323 ir=n 1324 1325 do ! Infinite do-loop 1326 if (l>1) then 1327 l=l-1 1328 ap=list(l) 1329 iap=iperm(l) 1330 1331 else ! l<=1 1332 ap=list(ir) 1333 iap=iperm(ir) 1334 list(ir)=list(1) 1335 iperm(ir)=iperm(1) 1336 ir=ir-1 1337 1338 if (ir==1) then 1339 list(1)=ap 1340 iperm(1)=iap 1341 exit ! This is the end of this algorithm 1342 end if 1343 end if ! l>1 1344 1345 i=l 1346 j=l+l 1347 1348 do while (j<=ir) 1349 if (j<ir) then 1350 if ( list(j)<list(j+1)-tol .or. & 1351 & (list(j)<list(j+1)+tol.and.iperm(j)<iperm(j+1))) j=j+1 1352 endif 1353 if (ap<list(j)-tol .or. (ap<list(j)+tol.and.iap<iperm(j))) then 1354 list(i)=list(j) 1355 iperm(i)=iperm(j) 1356 i=j 1357 j=j+j 1358 else 1359 j=ir+1 1360 end if 1361 enddo 1362 1363 list(i)=ap 1364 iperm(i)=iap 1365 1366 enddo ! End infinite do-loop 1367 1368 end if ! n>1 1369 1370 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.
SOURCE
1427 subroutine split_work(ntasks,comm,nprocs,my_start,my_stop,ierr) 1428 1429 !Arguments ------------------------------------ 1430 integer,intent(in) :: ntasks,comm 1431 integer,intent(out) :: nprocs,my_start,my_stop,ierr 1432 1433 !Local variables------------------------------- 1434 integer :: res,my_rank,block_p1,block,mpierr 1435 1436 ! ************************************************************************* 1437 1438 nprocs = 1; my_start = 1; my_stop = ntasks; ierr = 1 1439 #ifdef HAVE_MPI 1440 call MPI_COMM_SIZE(comm,nprocs,mpierr); if (mpierr /= MPI_SUCCESS) return 1441 call MPI_COMM_RANK(comm,my_rank,mpierr); if (mpierr /= MPI_SUCCESS) return 1442 1443 block = ntasks/nprocs 1444 res = MOD(ntasks,nprocs) 1445 block_p1= block+1 1446 1447 if (my_rank<res) then 1448 my_start = my_rank *block_p1+1 1449 my_stop = (my_rank+1)*block_p1 1450 else 1451 my_start = res*block_p1 + (my_rank-res )*block + 1 1452 my_stop = res*block_p1 + (my_rank-res+1)*block 1453 end if 1454 #endif 1455 ierr = 0 1456 1457 end subroutine split_work
m_tetrahedron/t_tetrahedron [ Types ]
[ Top ] [ m_tetrahedron ] [ Types ]
NAME
t_tetrahedron
FUNCTION
tetrahedron geometry object
SOURCE
71 type, public :: t_tetrahedron 72 73 integer :: ntetra = 0 74 ! Number of tetrahedra 75 76 real(dp) :: vv 77 ! volume of the tetrahedra 78 79 real(dp) :: klatt(3, 3) 80 ! reciprocal of lattice vectors for full kpoint grid 81 82 integer,allocatable :: tetra_full(:,:,:) 83 !(4,2,ntetra) 84 ! For each tetra 85 ! (:,1,itetra) indices of the vertex in IBZ (symmetrical image) 86 ! (:,2,itetra) indices of the vertexes in the BZ 87 88 integer,allocatable :: tetra_mult(:) 89 !(ntetra) 90 ! multiplicity of each irred tetrahedron 91 92 integer,allocatable :: tetra_wrap(:,:,:) 93 !(3,4,ntetra) 94 ! flag to wrap tetrahedron summit into IBZ 95 96 integer,allocatable :: ibz_tetra_count(:) 97 ! ibz_tetra_mapping(nkpt_ibz) 98 ! Number of tetrahedra associated to a point in the IBZ. 99 100 integer,allocatable :: ibz_tetra_mapping(:,:) 101 ! ibz_tetra_mapping(nkpt_ibz, maxval(tetra%ibz_tetra_count))) 102 ! map ikbz to tetra index. 103 104 end type t_tetrahedron 105 106 public :: init_tetra ! Initialize the object 107 ! See also the high-level interface tetra_from_kptrlatt provided by m_kpts. 108 public :: get_tetra_weight ! Calculate integration weights and their derivatives. shape (nkpt, nene). 109 public :: tetra_blochl_weights ! Same as in get_tetra_weight but weights have shape (nene, nkpt). 110 public :: get_dbl_tetra_weight ! Calculate integration weights for double tetrahedron integration of delta functions. 111 ! (NB these correspond to the derivative terms in normal tetrahedron). 112 public :: destroy_tetra ! Free memory. 113 public :: tetra_write ! Write text file (XML format) with tetra info. 114 public :: tetralib_has_mpi ! Return True if the library has been compiled with MPI support. 115 public :: tetra_get_onewk ! Calculate integration weights and their derivatives for a single k-point in the IBZ. 116 public :: tetra_get_onewk_wvals ! Similar to tetra_get_onewk_wvalsa but receives arbitrary list of frequency points. 117 public :: tetra_get_onetetra_wvals ! Get weights for one tetrahedra with arbitrary list of frequency points
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 output weights here have shape (nene, nkpt)
SOURCE
689 subroutine tetra_blochl_weights(tetra,eigen_in,enemin,enemax,max_occ,nene,nkpt,& 690 bcorr,tweight_t,dtweightde_t,comm) 691 692 !Arguments ------------------------------------ 693 !scalars 694 integer,intent(in) :: nene,nkpt,bcorr,comm 695 type(t_tetrahedron), intent(in) :: tetra 696 real(dp) ,intent(in) :: enemax,enemin,max_occ 697 !arrays 698 real(dp) ,intent(in) :: eigen_in(nkpt) 699 real(dp) ,intent(out) :: dtweightde_t(nene,nkpt),tweight_t(nene,nkpt) 700 701 !Local variables------------------------------- 702 !scalars 703 integer :: itetra,nprocs,my_start,my_stop,ierr,ii 704 !arrays 705 integer :: ind_ibz(4) 706 real(dp) :: eigen_1tetra(4) 707 real(dp), allocatable :: tweight_tmp(:,:),dtweightde_tmp(:,:),buffer(:,:) 708 709 ! ********************************************************************* 710 711 TETRA_ALLOCATE(tweight_tmp, (nene, 4)) 712 TETRA_ALLOCATE(dtweightde_tmp, (nene, 4)) 713 tweight_t = zero; dtweightde_t = zero 714 715 call split_work(tetra%ntetra, comm, nprocs, my_start, my_stop, ierr) 716 if (ierr /= 0) TETRA_ERROR("Error in MPI layer") 717 718 ! for each tetrahedron 719 do itetra=my_start,my_stop 720 tweight_tmp = zero 721 dtweightde_tmp = zero 722 723 ! Here we need the original ordering to reference the correct irred kpoints 724 ind_ibz(:) = tetra%tetra_full(:,1,itetra) 725 726 ! Sort energies before calling get_onetetra_ 727 eigen_1tetra(:) = eigen_in(ind_ibz(:)) 728 call sort_tetra(4, eigen_1tetra, ind_ibz, tol14) 729 730 call get_onetetra_(tetra,itetra,eigen_1tetra,enemin,enemax,max_occ,nene,bcorr,tweight_tmp,dtweightde_tmp) 731 732 ! NOTE: the following blas calls are not working systematically, or do not give speed ups, strange... 733 !if (nene > 100) then 734 ! do ii=1,4 735 ! call daxpy (nene, 1.d0, tweight_tmp(:,ii), 1, tweight_t(:,ind_ibz(ii)), 1) 736 ! end do 737 ! do ii=1,4 738 ! call daxpy (nene, 1.d0, dtweightde_tmp(:,ii), 1, dtweightde_t(:,ind_ibz(ii)), 1) 739 ! end do 740 !else 741 do ii=1,4 742 tweight_t(:,ind_ibz(ii)) = tweight_t(:,ind_ibz(ii)) + tweight_tmp(:,ii) 743 end do 744 do ii=1,4 745 dtweightde_t(:,ind_ibz(ii)) = dtweightde_t(:,ind_ibz(ii)) + dtweightde_tmp(:,ii) 746 end do 747 !end if 748 end do ! itetra 749 750 TETRA_DEALLOCATE(tweight_tmp) 751 TETRA_DEALLOCATE(dtweightde_tmp) 752 753 if (nprocs > 1) then 754 #ifdef HAVE_MPI 755 TETRA_ALLOCATE(buffer, (nene, nkpt)) 756 call MPI_ALLREDUCE(tweight_t,buffer,nene*nkpt,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) 757 tweight_t = buffer 758 759 call MPI_ALLREDUCE(dtweightde_t,buffer,nene*nkpt,MPI_DOUBLE_PRECISION,MPI_SUM,comm,ierr) 760 dtweightde_t = buffer 761 TETRA_DEALLOCATE(buffer) 762 #endif 763 end if 764 765 end subroutine tetra_blochl_weights
m_tetrahedron/tetra_get_onetetra_wvals [ Functions ]
[ Top ] [ m_tetrahedron ] [ Functions ]
NAME
tetra_get_onetetra_wvals
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. nw=number of energies in wvals nibz=number of irreducible kpoints wvals(nw)=Frequency points. eigen_ibz(nkibz)=eigenenergies for each k point [wtol]: If present, frequency points that differ by less that wtol are treated as equivalent. and the tetrahedron integration is performed only once per frequency point.
OUTPUT
weights(nw,2) = integration weights for
Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
for a given (band, k-point, spin).
SOURCE
1974 subroutine tetra_get_onetetra_wvals(tetra, itetra, eigen_1tetra, bcorr, nw, wvals, weights, wtol) 1975 1976 !Arguments ------------------------------------ 1977 !scalars 1978 integer,intent(in) :: nw,bcorr 1979 real(dp), optional, intent(in) :: wtol 1980 type(t_tetrahedron), intent(in) :: tetra 1981 !arrays 1982 real(dp),intent(in) :: wvals(nw) 1983 real(dp),intent(out) :: weights(nw, 2, 4) 1984 1985 !Local variables------------------------------- 1986 !scalars 1987 !integer,save :: done = 0 1988 integer,parameter :: nene3=3 1989 integer :: itetra,ii,idx,iw,ie 1990 integer :: ind(4) 1991 logical :: samew 1992 real(dp),parameter :: max_occ1 = one 1993 real(dp) :: enemin, enemax 1994 !arrays 1995 real(dp) :: theta_tmp(nene3,4), delta_tmp(nene3,4), eigen_1tetra(4) 1996 1997 ! ********************************************************************* 1998 1999 ind = [1,2,3,4] 2000 call sort_tetra(4, eigen_1tetra, ind, tol14) 2001 weights = 0 2002 2003 !for all the frequencies 2004 do iw=1,nw 2005 samew = .False. 2006 if (present(wtol)) then 2007 if (iw > 1) samew = abs(wvals(iw) - wvals(iw - 1)) < wtol 2008 end if 2009 if (.not. samew) then 2010 enemin = wvals(iw) - 0.01 2011 enemax = wvals(iw) + 0.01 2012 ie = nene3 / 2 + 1 2013 call get_onetetra_(tetra, itetra, eigen_1tetra, enemin, enemax, max_occ1, nene3, bcorr, & 2014 theta_tmp, delta_tmp) 2015 end if 2016 2017 ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz) 2018 do ii=1,4 2019 idx = ind(ii) 2020 weights(iw, 1, idx) = weights(iw, 1, idx) + delta_tmp(ie, ii) 2021 weights(iw, 2, idx) = weights(iw, 2, idx) + theta_tmp(ie, ii) 2022 end do 2023 end do !iw 2024 2025 end subroutine tetra_get_onetetra_wvals
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).
SOURCE
1807 subroutine tetra_get_onewk(tetra,ik_ibz,bcorr,nene,nkibz,eig_ibz,enemin,enemax,max_occ,weights) 1808 1809 !Arguments ------------------------------------ 1810 !scalars 1811 integer,intent(in) :: ik_ibz,nene,nkibz,bcorr 1812 type(t_tetrahedron), intent(in) :: tetra 1813 real(dp) ,intent(in) :: enemin,enemax,max_occ 1814 !arrays 1815 real(dp),intent(in) :: eig_ibz(nkibz) 1816 real(dp),intent(out) :: weights(nene,2) 1817 1818 !Local variables------------------------------- 1819 !scalars 1820 integer :: itetra,ii 1821 !arrays 1822 integer :: ind_ibz(4) 1823 real(dp) :: tweight_tmp(nene,4),dtweightde_tmp(nene,4),eigen_1tetra(4) 1824 1825 ! ********************************************************************* 1826 1827 weights = zero 1828 1829 ! For each tetrahedron 1830 do itetra=1,tetra%ntetra 1831 1832 ! Here we need the original ordering to reference the correct irred kpoints 1833 ind_ibz(:) = tetra%tetra_full(:,1,itetra) 1834 ! Cycle if this tetra does not contribute to this k-point. 1835 if (all(ind_ibz /= ik_ibz)) cycle 1836 1837 ! Sort energies before calling get_onetetra_ 1838 eigen_1tetra(:) = eig_ibz(ind_ibz(:)) 1839 call sort_tetra(4, eigen_1tetra, ind_ibz, tol14) 1840 1841 call get_onetetra_(tetra, itetra, eigen_1tetra, enemin, enemax, max_occ, nene, bcorr, & 1842 tweight_tmp, dtweightde_tmp) 1843 1844 ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz) 1845 do ii=1,4 1846 if (ind_ibz(ii) == ik_ibz) then 1847 weights(:,1) = weights(:,1) + dtweightde_tmp(:,ii) 1848 weights(:,2) = weights(:,2) + tweight_tmp(:,ii) 1849 end if 1850 end do 1851 end do ! itetra 1852 1853 end subroutine tetra_get_onewk
m_tetrahedron/tetra_get_onewk_wvals [ Functions ]
[ Top ] [ m_tetrahedron ] [ Functions ]
NAME
tetra_get_onewk_wvals
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. nw=number of energies in wvals nibz=number of irreducible kpoints wvals(nw)=Frequency points. eigen_ibz(nkibz)=eigenenergies for each k point [wtol]: If present, frequency points that differ by less that wtol are treated as equivalent. and the tetrahedron integration is performed only once per frequency point.
OUTPUT
weights(nw,2) = integration weights for
Dirac delta (derivative of theta wrt energy) and Theta (Heaviside function)
for a given (band, k-point, spin).
SOURCE
1883 subroutine tetra_get_onewk_wvals(tetra, ik_ibz, bcorr, nw, wvals, nkibz, eig_ibz, weights, wtol) 1884 1885 !Arguments ------------------------------------ 1886 !scalars 1887 integer,intent(in) :: ik_ibz,nw,nkibz,bcorr 1888 real(dp), optional, intent(in) :: wtol 1889 type(t_tetrahedron), intent(in) :: tetra 1890 !arrays 1891 real(dp),intent(in) :: wvals(nw) 1892 real(dp),intent(in) :: eig_ibz(nkibz) 1893 real(dp),intent(out) :: weights(nw, 2) 1894 1895 !Local variables------------------------------- 1896 !scalars 1897 !integer,save :: done = 0 1898 integer,parameter :: nene=3 1899 integer :: itetra,ii,jj,iw,ie 1900 logical :: samew 1901 real(dp),parameter :: max_occ1 = one 1902 real(dp) :: enemin, enemax 1903 !arrays 1904 integer :: ind_ibz(4) 1905 real(dp) :: theta_tmp(nene,4), delta_tmp(nene,4), eigen_1tetra(4) 1906 1907 ! ********************************************************************* 1908 1909 weights = zero 1910 1911 ! For each tetrahedron 1912 do jj=1,tetra%ibz_tetra_count(ik_ibz) 1913 itetra = tetra%ibz_tetra_mapping(ik_ibz,jj) 1914 1915 ! Here we need the original ordering to reference the correct irred kpoints 1916 ind_ibz(:) = tetra%tetra_full(:,1,itetra) 1917 1918 ! Sort energies before calling get_onetetra_ 1919 eigen_1tetra(:) = eig_ibz(ind_ibz(:)) 1920 call sort_tetra(4, eigen_1tetra, ind_ibz, tol14) 1921 1922 do iw=1,nw 1923 samew = .False. 1924 if (present(wtol)) then 1925 if (iw > 1) samew = abs(wvals(iw) - wvals(iw - 1)) < wtol 1926 end if 1927 if (.not. samew) then 1928 enemin = wvals(iw) - 0.01; enemax = wvals(iw) + 0.01 1929 ie = nene / 2 + 1 1930 call get_onetetra_(tetra, itetra, eigen_1tetra, enemin, enemax, max_occ1, nene, bcorr, & 1931 theta_tmp, delta_tmp) 1932 end if 1933 1934 ! Accumulate contributions to ik_ibz (there might be multiple vertexes that map onto ik_ibz) 1935 do ii=1,4 1936 if (ind_ibz(ii) == ik_ibz) then 1937 weights(iw, 1) = weights(iw, 1) + delta_tmp(ie, ii) 1938 weights(iw, 2) = weights(iw, 2) + theta_tmp(ie, ii) 1939 end if 1940 end do 1941 end do ! iw 1942 end do ! itetra 1943 1944 end subroutine tetra_get_onewk_wvals
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.
SOURCE
547 subroutine tetra_write(tetra, nkibz, kibz, path) 548 549 !Arguments ------------------------------------ 550 !scalars 551 integer,intent(in) :: nkibz 552 character(len=*),intent(in) :: path 553 type(t_tetrahedron),intent(in) :: tetra 554 !arrays 555 real(dp),intent(in) :: kibz(3,nkibz) 556 557 !Local variables------------------------------- 558 integer,parameter :: version=1 559 integer :: ik,it,unt 560 #ifdef HAVE_LIBTETRA_ABINIT 561 character(len=500) :: msg 562 #endif 563 564 ! ********************************************************************* 565 566 #ifdef HAVE_LIBTETRA_ABINIT 567 if (open_file(file=trim(path),iomsg=msg,newunit=unt,form="formatted",status="unknown",action="write")/=0) then 568 TETRA_ERROR(msg) 569 end if 570 #else 571 open(file=trim(path),newunit=unt,form="formatted",status="unknown",action="write") 572 #endif 573 574 write(unt,*)version, " # version number" 575 576 ! Write IBZ 577 write(unt,*)nkibz, " # number of k-points in the IBZ" 578 write(unt,"(a)")"<irreducible_zone>" 579 do ik=1,nkibz 580 write(unt,"(3es22.12)") kibz(:,ik) 581 end do 582 write(unt,"(a)")"</irreducible_zone>" 583 584 ! Write tetra info 585 write(unt,"(i0,a)")tetra%ntetra, " # number of tetrahedra" 586 write(unt,"(es22.12,a)")tetra%vv, " # tetrahedron volume" 587 588 write(unt,"(a)")"<tetra_full>" 589 do it=1,tetra%ntetra 590 write(unt,"(8(i0,1x))")tetra%tetra_full(:,:,it) 591 end do 592 write(unt,"(a)")"</tetra_full>" 593 594 write(unt,"(a)")"<tetra_mult>" 595 do it=1,tetra%ntetra 596 write(unt,"(i0)")tetra%tetra_mult(it) 597 end do 598 write(unt,"(a)")"</tetra_mult>" 599 600 write(unt,"(a)")"<tetra_wrap>" 601 do it=1,tetra%ntetra 602 write(unt,"(12(i0,1x))")tetra%tetra_wrap(:,:,it) 603 end do 604 write(unt,"(a)")"</tetra_wrap>" 605 606 close(unt) 607 608 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
SOURCE
1384 logical function tetralib_has_mpi() result(ans) 1385 1386 ans = .False. 1387 #ifdef HAVE_MPI 1388 ans = .True. 1389 #endif 1390 1391 end function tetralib_has_mpi