TABLE OF CONTENTS
- ABINIT/m_kg
- ABINIT/mkkin_metdqdq
- ABINIT/mkpwind_k
- m_kg/getcut
- m_kg/getmpw
- m_kg/getph
- m_kg/kpgio
- m_kg/kpgstr
- m_kg/mkkin
- m_kg/mkkpg
- m_kg/mkkpgcart
- m_kg/ph1d3d
ABINIT/m_kg [ Modules ]
NAME
m_kg
FUNCTION
Low-level functions to operate of G-vectors.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (DCA, XG, GMR, MT, DRH, AR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_kg 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 use m_dtset 30 31 use defs_abitypes, only : MPI_type 32 use m_fftcore, only : kpgsph, bound 33 use m_mpinfo, only : proc_distrb_cycle 34 use m_time, only : timab 35 36 implicit none 37 38 private
ABINIT/mkkin_metdqdq [ Functions ]
NAME
mkkin_metdqdq
FUNCTION
Compute elements of the second q-gradient of the metric kinetic energy operator in reciprocal space at given k point wrt cartesian q components.
INPUTS
effmass=effective mass for electrons (1. in common case) gprimd(3,3)=reciprocal space dimensional primitive translations idir= strain perturbation direction kg(3,npw) = integer coordinates of planewaves in basis sphere. kpt(3) = reduced coordinates of k point npw = number of plane waves at kpt. qdir = direction of the first q-gradient
OUTPUT
dqdqkinpw(npw)=d/deps(istr) ( (1/2)*(2 pi)**2 * (k+G)**2 )
NOTES
**Since the 2nd derivative w.r.t q-vector is calculated along cartesian
directions, the 1/twopi**2 factor (that in the rest of the code is applied
in the reduced to cartesian derivative conversion process) is here
explicictly included in the formulas.
**Notice that idir=1-9, in contrast to the strain perturbation (idir=1-6),
because this term is not symmetric w.r.t permutations of the two strain
indices.
**A -i factor has been factorized out in all the contributions of the second
q-gradient of the metric Hamiltonian. This is lately included in the contribution
of the corresponing term (T4) to the flexoelectric tensor in dfpt_flexoout.F90
SOURCE
1218 subroutine mkkin_metdqdq(dqdqkinpw,effmass,gprimd,idir,kg,kpt,npw,qdir) 1219 1220 !Arguments ------------------------------- 1221 !scalars 1222 integer,intent(in) :: idir,npw,qdir 1223 real(dp),intent(in) :: effmass 1224 !arrays 1225 integer,intent(in) :: kg(3,npw) 1226 real(dp),intent(in) :: gprimd(3,3),kpt(3) 1227 real(dp),intent(out) :: dqdqkinpw(npw) 1228 1229 !Local variables ------------------------- 1230 !scalars 1231 integer :: beta,delta,gamma,ig,ka,kb 1232 real(dp) :: delbd,delbg,deldg 1233 real(dp) :: dkinetic,gpk1,gpk2,gpk3,htpi 1234 real(dp) :: gpkc(3) 1235 !arrays 1236 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 1237 1238 ! ********************************************************************* 1239 1240 !htpi is (1/2) (2 Pi): 1241 htpi=0.5_dp*two_pi 1242 1243 ka=idx(2*idir-1);kb=idx(2*idir) 1244 1245 !For easier formula implementation 1246 beta=ka 1247 delta=kb 1248 gamma=qdir 1249 1250 !Kronecker deltas 1251 delbd=0.0_dp; delbg=0.0_dp; deldg=0.0_dp 1252 if (beta==delta) delbd=1.0_dp 1253 if (beta==gamma) delbg=1.0_dp 1254 if (delta==gamma) deldg=1.0_dp 1255 1256 do ig=1,npw 1257 gpk1=dble(kg(1,ig))+kpt(1) 1258 gpk2=dble(kg(2,ig))+kpt(2) 1259 gpk3=dble(kg(3,ig))+kpt(3) 1260 1261 ! Obtain G in cartesian coordinates 1262 gpkc(1)=gprimd(1,1)*gpk1+gprimd(1,2)*gpk2+gprimd(1,3)*gpk3 1263 gpkc(2)=gprimd(2,1)*gpk1+gprimd(2,2)*gpk2+gprimd(2,3)*gpk3 1264 gpkc(3)=gprimd(3,1)*gpk1+gprimd(3,2)*gpk2+gprimd(3,3)*gpk3 1265 1266 dkinetic=htpi*(2.0_dp*deldg*gpkc(beta)+delbg*gpkc(delta)+ & 1267 & delbd*gpkc(gamma)) 1268 1269 dqdqkinpw(ig)=dkinetic/effmass 1270 end do 1271 1272 end subroutine mkkin_metdqdq
ABINIT/mkpwind_k [ Functions ]
NAME
mkpwind_k
FUNCTION
Make plane wave index at k point for basis at second k point,
needed to compute overlaps $\langle u_{k,n}|u_{k+b,n}\rangle$
as appear in Berry phase derived quantities
INPUTS
dk(3)=real vector difference of ket kpt - bra kpt dtset <type(dataset_type)>=all input variables in this dataset fnkpt=number of kpts in full BZ fkptns=kpts in full BZ gmet(3,3)=metric in reciprocal space indkk_f2ibz(fnkpt,6)=information on folding from FBZ to IBZ (see initberry or initorbmag) ikpt=index of bra k pt in FBZ ikpt1=index of neighbour ket k pt in FBZ mpi_enreg=information about MPI parallelization npwarr(dtset%nkpt)=npw at each kpt symrec(3,3,nsym) = symmetries in reciprocal space in terms of reciprocal space primitive translations
OUTPUT
pwind_k1(dtset%mpw)=output index of ikpt1 basis states refered to ikpt
SIDE EFFECTS
TODO
NOTES
SOURCE
986 subroutine mkpwind_k(dk,dtset,fnkpt,fkptns,gmet,indkk_f2ibz,ikpt,ikpt1,& 987 & mpi_enreg,npwarr,pwind_k1,symrec) 988 989 !Arguments ------------------------------------ 990 !scalars 991 integer,intent(in) :: fnkpt,ikpt,ikpt1 992 type(dataset_type),intent(in) :: dtset 993 type(MPI_type), intent(inout) :: mpi_enreg 994 995 !arrays 996 integer,intent(in) :: indkk_f2ibz(fnkpt,6) 997 integer,intent(in) :: npwarr(dtset%nkpt) 998 integer,intent(in) :: symrec(3,3,dtset%nsym) 999 integer,intent(out) :: pwind_k1(dtset%mpw) 1000 real(dp),intent(in) :: dk(3),fkptns(3,fnkpt),gmet(3,3) 1001 1002 !Local variables ------------------------- 1003 !scalars 1004 integer :: exchn2n3d,idum1,ikg1,ikpti,ikpt1i,ipw,istwf_k,isym,isym1,jpw,npw_k,npw_k1 1005 real(dp) :: ecut_eff 1006 1007 !arrays 1008 integer,allocatable :: kg_k(:,:),kg1_k(:,:) 1009 real(dp) :: dg(3),dum33(3,3),kpt(3),kpt1(3),iadum(3),iadum1(3) 1010 1011 ! *********************************************************************** 1012 1013 ikpti = indkk_f2ibz(ikpt,1) 1014 ikpt1i = indkk_f2ibz(ikpt1,1) 1015 1016 1017 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 1018 exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0 1019 1020 ! Build basis sphere of plane waves for the k-point 1021 ! we avoid using the global kg data because of difficulties in parallel-ism 1022 ABI_MALLOC(kg_k,(3,dtset%mpw)) 1023 kg_k(:,:) = 0 1024 kpt(:) = dtset%kptns(:,ikpti) 1025 call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg_k,kpt,1,mpi_enreg,dtset%mpw,npw_k) 1026 1027 ! Build basis sphere of plane waves for the nearest neighbour of the k-point 1028 ABI_MALLOC(kg1_k,(3,dtset%mpw)) 1029 kg1_k(:,:) = 0 1030 kpt1(:) = dtset%kptns(:,ikpt1i) 1031 call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,1,mpi_enreg,dtset%mpw,npw_k1) 1032 1033 ! Deal with symmetry transformations 1034 ! 1035 1036 ! bra k-point k(b) and IBZ k-point kIBZ(b) related by 1037 ! k(b) = alpha(b) S(b)^t kIBZ(b) + G(b) 1038 ! where alpha(b), S(b) and G(b) are given by indkk_f2ibz 1039 ! 1040 ! For the ket k-point: 1041 ! k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k) 1042 ! where GBZ(k) takes k(k) to the BZ 1043 ! 1044 1045 isym = indkk_f2ibz(ikpt,2) 1046 isym1 = indkk_f2ibz(ikpt1,2) 1047 1048 ! Construct transformed G vector that enters the matching condition: 1049 ! alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) ) 1050 1051 dg(:) = -indkk_f2ibz(ikpt,3:5) & 1052 & - nint(-fkptns(:,ikpt) - dk(:) - tol10 + fkptns(:,ikpt1)) & 1053 & + indkk_f2ibz(ikpt1,3:5) 1054 1055 iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:)) 1056 1057 dg(:) = iadum(:) 1058 1059 ! Construct S(k)^{t,-1} S(b)^{t} 1060 1061 dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym)) 1062 1063 ! Construct alpha(k) alpha(b) 1064 1065 pwind_k1(:) = 0 1066 npw_k = npwarr(ikpti) 1067 do ipw = 1, npw_k 1068 1069 ! NOTE: the bra G vector is taken for the sym-related IBZ k point, 1070 ! not for the FBZ k point 1071 1072 ! original code from initberry 1073 ! iadum(:) = kg(:,kgindex(ikpti) + ipw) 1074 1075 iadum(:) = kg_k(:,ipw) 1076 1077 ! to determine r.l.v. matchings, we transformed the bra vector 1078 ! Rotation 1079 iadum1(:)=0 1080 do idum1=1,3 1081 iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1) 1082 end do 1083 iadum(:)=iadum1(:) 1084 iadum(:) = iadum(:) + dg(:) 1085 1086 do jpw = 1, npw_k1 1087 iadum1(1:3) = kg1_k(1:3,jpw) 1088 if ( (iadum(1) == iadum1(1)).and. & 1089 & (iadum(2) == iadum1(2)).and. & 1090 & (iadum(3) == iadum1(3)) ) then 1091 pwind_k1(ipw) = jpw 1092 exit 1093 end if 1094 end do 1095 end do 1096 1097 ABI_FREE(kg_k) 1098 ABI_FREE(kg1_k) 1099 1100 end subroutine mkpwind_k
m_kg/getcut [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
getcut
FUNCTION
For input kpt, fft box dim ngfft(1:3), recip space metric gmet, and kinetic energy cutoff ecut, COMPUTES: if iboxcut==0: gsqcut: cut-off on G^2 for "large sphere" of radius double that of the basis sphere corresponding to ecut boxcut: where boxcut == gcut(box)/gcut(sphere). boxcut >=2 for no aliasing. boxcut < 1 is wrong and halts subroutine. if iboxcut==1: gsqcut: cut-off on G^2 for "large sphere" containing the whole fft box boxcut: no meaning (zero)
INPUTS
ecut=kinetic energy cutoff for planewave sphere (hartree) gmet(3,3)=reciprocal space metric (bohr^-2) iboxcut=0: compute gsqcut and boxcut with boxcut>=1 1: compute gsqcut for boxcut=1 (sphere_cutoff=box_cutoff) iout=unit number for output file kpt(3)=input k vector (reduced coordinates--in terms of reciprocal lattice primitive translations) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
OUTPUT
boxcut=defined above (dimensionless), ratio of basis sphere diameter to fft box length (smallest value) gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere--appropriate for charge density rho(G), Hartree potential, and pseudopotentials
NOTES
2*gcut arises from rho(g)=sum g prime (psi(g primt)*psi(g prime+g))
where psi(g) is only nonzero for |g| <= gcut).
ecut (currently in hartree) is proportional to gcut(sphere)**2.
SOURCE
96 subroutine getcut(boxcut, ecut, gmet, gsqcut, iboxcut, iout, kpt, ngfft) 97 98 !Arguments ------------------------------------ 99 !scalars 100 integer,intent(in) :: iboxcut,iout 101 real(dp),intent(in) :: ecut 102 real(dp),intent(out) :: boxcut,gsqcut 103 !arrays 104 integer,intent(in) :: ngfft(18) 105 real(dp),intent(in) :: gmet(3,3),kpt(3) 106 107 !Local variables------------------------------- 108 !scalars 109 integer :: plane 110 real(dp) :: boxsq,cutrad,ecut_pw,effcut,largesq,sphsq 111 character(len=1000) :: msg 112 !arrays 113 integer :: gbound(3) 114 115 ! ************************************************************************* 116 117 ! This is to treat the case in which ecut has not been initialized e.g. for wavelet computations. 118 ! The default for ecut is -1.0 , allowed only for wavelets calculations 119 ecut_pw=ecut 120 if(ecut<-tol8)ecut_pw=ten 121 122 !gcut(box)**2=boxsq; gcut(sphere)**2=sphsq 123 !get min. d**2 to boundary of fft box: 124 !(gmet sets dimensions: bohr**-2) 125 !ecut(sphere)=0.5*(2 pi)**2 * sphsq: 126 call bound(largesq,boxsq,gbound,gmet,kpt,ngfft,plane) 127 effcut=0.5_dp * (two_pi)**2 * boxsq 128 sphsq=2._dp*ecut_pw/two_pi**2 129 130 if (iboxcut/=0) then 131 boxcut=10._dp 132 gsqcut=(largesq/sphsq)*(2.0_dp*ecut)/two_pi**2 133 134 write(msg, '(a,a,3f8.4,a,3i4,a,a,f11.3,a,a)' ) ch10,& 135 ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,& 136 ' ecut(hartree)=',ecut_pw+tol8,ch10,'=> whole FFT box selected' 137 if(iout/=std_out) call wrtout(iout,msg) 138 call wrtout(std_out,msg) 139 else 140 141 ! Get G^2 cutoff for sphere of double radius of basis sphere 142 ! for selecting G s for rho(G), V_Hartree(G), and V_psp(G)-- 143 ! cut off at fft box boundary or double basis sphere radius, whichever 144 ! is smaller. If boxcut were 2, then relation would be 145 ! $ecut_eff = (1/2) * (2 Pi Gsmall)^2 and gsqcut=4*Gsmall^2$. 146 boxcut = sqrt(boxsq/sphsq) 147 cutrad = min(2.0_dp,boxcut) 148 gsqcut = (cutrad**2)*(2.0_dp*ecut_pw)/two_pi**2 149 150 if(ecut>-tol8)then 151 152 write(msg, '(a,a,3f8.4,a,3i4,a,a,f11.3,3x,a,f10.5)' ) ch10,& 153 ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,& 154 ' ecut(hartree)=',ecut+tol8,'=> boxcut(ratio)=',boxcut+tol8 155 if(iout/=std_out) call wrtout(iout,msg) 156 call wrtout(std_out,msg) 157 158 if (boxcut<1.0_dp) then 159 write(msg, '(9a,f12.6,6a)' )& 160 'Choice of acell, ngfft, and ecut',ch10,& 161 '===> basis sphere extends BEYOND fft box !',ch10,& 162 'Recall that boxcut=Gcut(box)/Gcut(sphere) must be > 1.',ch10,& 163 'Action: try larger ngfft or smaller ecut.',ch10,& 164 'Note that ecut=effcut/boxcut**2 and effcut=',effcut+tol8,ch10,& 165 'This situation might happen when optimizing the cell parameters.',ch10,& 166 'Your starting geometry might be crazy.',ch10,& 167 'See https://wiki.abinit.org/doku.php?id=howto:troubleshooting#incorrect_initial_geometry .' 168 if(iout/=std_out) call wrtout(iout,msg) 169 ABI_ERROR(msg) 170 end if 171 172 if (boxcut>2.2_dp) then 173 write(msg, '(a,a,a,a,a,a,a,a,a,a,a,f12.6,a,a)' ) ch10,& 174 ' getcut : COMMENT -',ch10,& 175 ' Note that boxcut > 2.2 ; recall that',' boxcut=Gcut(box)/Gcut(sphere) = 2',ch10,& 176 ' is sufficient for exact treatment of convolution.',ch10,& 177 ' Such a large boxcut is a waste : you could raise ecut',ch10,& 178 ' e.g. ecut=',effcut*0.25_dp+tol8,' Hartrees makes boxcut=2',ch10 179 if(iout/=std_out) call wrtout(iout,msg) 180 call wrtout(std_out,msg) 181 end if 182 183 if (boxcut<1.5_dp) then 184 write(msg, '(15a)' ) ch10,& 185 ' getcut : WARNING -',ch10,& 186 ' Note that boxcut < 1.5; this usually means',ch10,& 187 ' that the forces are being fairly strongly affected by',' the smallness of the fft box.',ch10,& 188 ' Be sure to test with larger ngfft(1:3) values.',ch10,& 189 ' This situation might happen when optimizing the cell parameters.',ch10,& 190 ' Your starting geometry might be crazy.',ch10,& 191 ' See https://wiki.abinit.org/doku.php?id=howto:troubleshooting#incorrect_initial_geometry .' 192 if(iout/=std_out) call wrtout(iout,msg) 193 call wrtout(std_out,msg) 194 end if 195 196 end if 197 198 end if ! iboxcut 199 200 end subroutine getcut
m_kg/getmpw [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
getmpw
FUNCTION
From input ecut, combined with ucvol and gmet, compute recommended mpw mpw is the maximum number of plane-waves in the wave-function basis for one processor of the WF group
INPUTS
ecut=plane-wave cutoff energy in Hartrees exchn2n3d=if n1, n2 and n3 are exchanged gmet(3,3)=reciprocal space metric (bohr**-2). istwfk(nkpt)=input option parameter that describes the storage of wfs kptns(3,nkpt)=real(dp) array for k points (normalisation is already taken into account) mpi_enreg=information about MPI parallelization nkpt=integer number of k points in the calculation
OUTPUT
mpw=maximal number of plane waves over all k points of the processor (for one processor of the WF group)
SOURCE
227 subroutine getmpw(ecut,exchn2n3d,gmet,istwfk,kptns,mpi_enreg,mpw,nkpt) 228 229 !Arguments ------------------------------------ 230 !scalars 231 integer,intent(in) :: exchn2n3d,nkpt 232 integer,intent(out) :: mpw 233 real(dp),intent(in) :: ecut 234 type(MPI_type),intent(inout) :: mpi_enreg 235 !arrays 236 integer,intent(in) :: istwfk(nkpt) 237 real(dp),intent(in) :: gmet(3,3),kptns(3,nkpt) 238 239 !Local variables------------------------------- 240 !scalars 241 integer :: ikpt,istwf_k,npw 242 ! integer :: npwwrk,pad=50 243 ! real(dp) :: scale=1.3_dp 244 character(len=500) :: msg 245 !arrays 246 integer,allocatable :: kg(:,:) 247 real(dp) :: kpoint(3) 248 249 ! ************************************************************************* 250 251 !An upper bound for mpw, might be obtained as follows 252 !the average number of plane-waves in the cutoff sphere is 253 !npwave = (2*ecut)**(3/2) * ucvol / (6*pi**2) 254 !the upper bound is calculated as 255 !npwwrk = int(scale * npwave) + pad 256 !rescale so an upper bound 257 !npwave=nint(ucvol*(2.0_dp*ecut)**1.5_dp/(6.0_dp*pi**2)) 258 !npwwrk=nint(dble(npwave)*scale)+pad 259 260 ABI_MALLOC(kg,(3,100)) 261 262 !set mpw to zero, as needed for only counting in kpgsph 263 mpw = 0 264 265 !Might be parallelized over k points ? ! 266 do ikpt = 1,nkpt 267 ! Do computation of G sphere, returning npw 268 kpoint(:)=kptns(:,ikpt) 269 istwf_k=istwfk(ikpt) 270 call kpgsph(ecut,exchn2n3d,gmet,0,ikpt,istwf_k,kg,kpoint,0,mpi_enreg,0,npw) 271 mpw = max(npw,mpw) 272 end do 273 274 write(msg,'(a,i0)') ' getmpw: optimal value of mpw= ',mpw 275 call wrtout(std_out,msg) 276 277 ABI_FREE(kg) 278 279 end subroutine getmpw
m_kg/getph [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
getph
FUNCTION
Compute three factors of one-dimensional structure factor phase for input atomic coordinates, for all planewaves which fit in fft box. The storage of these atomic factors is made according to the values provided by the index table atindx. This will save time in nonlop.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) natom=number of atoms in cell. n1,n2,n3=dimensions of fft box (ngfft(3)). xred(3,natom)=reduced atomic coordinates.
OUTPUT
ph1d(2,(2*n1+1)*natom+(2*n2+1)*natom+(2*n3+1)*natom)=exp(2Pi i G.xred) for integer vector G with components ranging from -nj <= G <= nj. Real and imag given in usual Fortran convention.
SOURCE
708 subroutine getph(atindx, natom, n1, n2, n3, ph1d, xred) 709 710 !Arguments ------------------------------------ 711 !scalars 712 integer,intent(in) :: n1,n2,n3,natom 713 !arrays 714 integer,intent(in) :: atindx(natom) 715 real(dp),intent(in) :: xred(3,natom) 716 real(dp),intent(out) :: ph1d(:,:) 717 718 !Local variables------------------------------- 719 !scalars 720 integer,parameter :: im=2,re=1 721 integer :: i1,i2,i3,ia,ii,ph1d_size1,ph1d_size2,ph1d_sizemin 722 !character(len=500) :: msg 723 real(dp) :: arg 724 725 ! ************************************************************************* 726 727 ph1d_size1=size(ph1d,1); ph1d_size2=size(ph1d,2) 728 ph1d_sizemin = (2*n1+1+2*n2+1+2*n3+1)*natom 729 if (ph1d_size1 /= 2 .or. ph1d_size2 < ph1d_sizemin) then 730 ABI_BUG('Wrong ph1d sizes!') 731 end if 732 733 do ia=1,natom 734 735 ! Store the phase factor of atom number ia in place atindx(ia) 736 i1=(atindx(ia)-1)*(2*n1+1) 737 i2=(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1) 738 i3=(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 739 740 do ii=1,2*n1+1 741 arg=two_pi*dble(ii-1-n1)*xred(1,ia) 742 ph1d(re,ii+i1)=dcos(arg) 743 ph1d(im,ii+i1)=dsin(arg) 744 end do 745 746 do ii=1,2*n2+1 747 arg=two_pi*dble(ii-1-n2)*xred(2,ia) 748 ph1d(re,ii+i2)=dcos(arg) 749 ph1d(im,ii+i2)=dsin(arg) 750 end do 751 752 do ii=1,2*n3+1 753 arg=two_pi*dble(ii-1-n3)*xred(3,ia) 754 ph1d(re,ii+i3)=dcos(arg) 755 ph1d(im,ii+i3)=dsin(arg) 756 end do 757 758 end do 759 760 ! This to avoid uninitialized ph1d values 761 if (ph1d_sizemin < ph1d_size2) ph1d(:,ph1d_sizemin+1:ph1d_size2)=zero 762 763 end subroutine getph
m_kg/kpgio [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
kpgio
FUNCTION
Do initialization of kg data.
INPUTS
ecut=kinetic energy planewave cutoff (hartree) exchn2n3d=if 1, n2 and n3 are exchanged gmet(3,3)=reciprocal space metric (bohr^-2) istwfk(nkpt)=input option parameter that describes the storage of wfs kptns(3,nkpt)=reduced coords of k points mkmem =number of k points treated by this node. character(len=4) : mode_paral=either 'COLL' or 'PERS', tells whether the loop over k points must be done by all processors or not, in case of parallel execution. mpi_enreg=information about MPI parallelization mpw=maximum number of planewaves as dimensioned in calling routine nband(nkpt*nsppol)=number of bands at each k point nkpt=number of k points nsppol=1 for unpolarized, 2 for polarized
OUTPUT
npwarr(nkpt)=array holding npw for each k point, taking into account the effect of istwfk, and the spreading over processors npwtot(nkpt)=array holding the total number of plane waves for each k point, kg(3,mpw*mkmem)=dimensionless coords of G vecs in basis sphere at k point
NOTES
Note that in case of band parallelism, the number of spin-up and spin-down bands must be equal at each k points
SOURCE
471 subroutine kpgio(ecut,exchn2n3d,gmet,istwfk,kg,kptns,mkmem,nband,nkpt,& 472 & mode_paral,mpi_enreg,mpw,npwarr,npwtot,nsppol) 473 474 !Arguments ------------------------------------ 475 !scalars 476 integer,intent(in) :: exchn2n3d,mkmem,mpw,nkpt,nsppol 477 real(dp),intent(in) :: ecut 478 character(len=4),intent(in) :: mode_paral 479 type(MPI_type),intent(inout) :: mpi_enreg 480 !arrays 481 integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol) 482 integer,intent(out) :: kg(3,mpw*mkmem),npwarr(nkpt),npwtot(nkpt) 483 real(dp),intent(in) :: gmet(3,3),kptns(3,nkpt) 484 485 !Local variables------------------------------- 486 !scalars 487 integer :: ierr,ikg,ikpt,istwf_k,me,nband_down,nband_k,npw1 488 logical :: test_npw 489 character(len=500) :: msg 490 !arrays 491 real(dp) :: kpoint(3) 492 493 ! ************************************************************************* 494 495 !Define me 496 me=mpi_enreg%me_kpt 497 498 if((mpi_enreg%paralbd==1) .and. (mode_paral=='PERS')) then 499 if(nsppol==2)then 500 do ikpt=1,nkpt 501 nband_k=nband(ikpt) 502 nband_down=nband(ikpt+nkpt) 503 if(nband_k/=nband_down)then 504 write(msg,'(a,a,a,a,a,a,a,a,i4,a,i4,a,a,a,i4,a,a,a)')ch10,& 505 ' kpgio : ERROR -',ch10,& 506 ' Band parallel case, one must have same number',ch10,& 507 ' of spin up and spin down bands, but input is :',ch10,& 508 ' nband(up)=',nband_k,', nband(down)=',nband_down,',',ch10,& 509 ' for ikpt=',ikpt,'.',ch10,& 510 ' Action: correct nband in your input file.' 511 ! MG: Tests v3(10,11,17) and v6(67) fail if this test is enabled 512 ! call wrtout(std_out,msg,mode_paral) 513 end if 514 end do 515 end if 516 end if 517 npwarr(:)=0 518 npwtot(:)=0 519 520 kg=0 521 ikg=0 522 !Find (k+G) sphere for each k. 523 524 do ikpt=1,nkpt 525 526 nband_k = nband(ikpt) 527 528 if(mode_paral=='PERS')then 529 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,me)) cycle 530 end if 531 532 kpoint(:)=kptns(:,ikpt) 533 istwf_k=istwfk(ikpt) 534 call kpgsph(ecut,exchn2n3d,gmet,ikg,ikpt,istwf_k,kg,kpoint,mkmem,mpi_enreg,mpw,npw1) 535 536 test_npw=.true. 537 if (xmpi_paral==1)then 538 if (mode_paral=='PERS')then 539 test_npw=(minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,1:nsppol))==me) 540 end if 541 end if 542 if (test_npw) npwarr(ikpt)=npw1 543 544 ! Make sure npw < nband never happens: 545 ! if (npw1<nband(ikpt)) then 546 ! write(msg, '(a,a,a,a,i5,a,3f8.4,a,a,i10,a,i10,a,a,a,a)' )ch10,& 547 ! & ' kpgio : ERROR -',ch10,& 548 ! & ' At k point number',ikpt,' k=',(kptns(ierr,ikpt),ierr=1,3),ch10,& 549 ! & ' npw=',npw1,' < nband=',nband(ikpt),ch10,& 550 ! & ' Indicates not enough planewaves for desired number of bands.',ch10,& 551 ! & ' Action: change either ecut or nband in input file.' 552 ! ABI_ERROR(msg) 553 ! end if 554 555 ! Find boundary of G sphere for efficient zero padding, 556 ! Shift to next section of each array kg 557 ikg=ikg+npw1 558 end do ! End of the loop over k points 559 560 ! TODO: this fails on some platforms if nproc > nkpt 561 if(mode_paral == 'PERS') then 562 call xmpi_sum(npwarr,mpi_enreg%comm_kpt,ierr) 563 end if 564 565 !if (mpi_enreg%nproc>1) call wrtout(std_out,' kpgio: loop on k-points done in parallel','COLL') 566 567 !XG030513 MPIWF : now, one should sum npwarr over all processors 568 !of the WF group, to get npwtot (to be spread on all procs of the WF group 569 npwtot(:)=npwarr(:) 570 571 !Taking into account istwfk 572 do ikpt=1,nkpt 573 if(istwfk(ikpt)>1)then 574 if(istwfk(ikpt)==2)then 575 npwtot(ikpt)=2*npwtot(ikpt)-1 576 else 577 npwtot(ikpt)=2*npwtot(ikpt) 578 end if 579 end if 580 end do 581 582 end subroutine kpgio
m_kg/kpgstr [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
kpgstr
FUNCTION
Compute elements of the derivative the kinetic energy operator in reciprocal space at given k point wrt a single cartesian strain component
INPUTS
ecut=cut-off energy for plane wave basis sphere (Ha) ecutsm=smearing energy for plane wave kinetic energy (Ha) effmass_free=effective mass for electrons (1. in common case) gmet(3,3) = reciprocal lattice metric tensor (Bohr**-2) gprimd(3,3)=reciprocal space dimensional primitive translations istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21 kg(3,npw) = integer coordinates of planewaves in basis sphere. kpt(3) = reduced coordinates of k point npw = number of plane waves at kpt.
OUTPUT
dkinpw(npw)=d/deps(istr) ( (1/2)*(2 pi)**2 * (k+G)**2 )
NOTES
Src_6response/kpg3.f
SOURCE
793 subroutine kpgstr(dkinpw,ecut,ecutsm,effmass_free,gmet,gprimd,istr,kg,kpt,npw) 794 795 !Arguments ------------------------------- 796 !scalars 797 integer,intent(in) :: istr,npw 798 real(dp),intent(in) :: ecut,ecutsm,effmass_free 799 !arrays 800 integer,intent(in) :: kg(3,npw) 801 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3) 802 real(dp),intent(out) :: dkinpw(npw) 803 804 !Local variables ------------------------- 805 !scalars 806 integer :: ig,ii,ka,kb 807 real(dp) :: dfsm,dkinetic,dkpg2,ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq 808 ! real(dp) :: d2fsm ! used in commented section below 809 real(dp) :: kpg2,xx 810 character(len=500) :: msg 811 !arrays 812 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 813 real(dp) :: dgmetds(3,3) 814 815 ! ********************************************************************* 816 817 !htpisq is (1/2) (2 Pi) **2: 818 htpisq=0.5_dp*(two_pi)**2 819 820 ecutsm_inv=0.0_dp 821 if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm 822 823 !Compute derivative of metric tensor wrt strain component istr 824 if(istr<1 .or. istr>6)then 825 write(msg, '(a,i10,a,a,a)' )& 826 'Input istr=',istr,' not allowed.',ch10,& 827 'Possible values are 1,2,3,4,5,6 only.' 828 ABI_BUG(msg) 829 end if 830 831 ka=idx(2*istr-1);kb=idx(2*istr) 832 do ii = 1,3 833 dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 834 end do 835 !For historical reasons: 836 dgmetds(:,:)=0.5_dp*dgmetds(:,:) 837 838 do ig=1,npw 839 gpk1=dble(kg(1,ig))+kpt(1) 840 gpk2=dble(kg(2,ig))+kpt(2) 841 gpk3=dble(kg(3,ig))+kpt(3) 842 kpg2=htpisq*& 843 & ( gmet(1,1)*gpk1**2+ & 844 & gmet(2,2)*gpk2**2+ & 845 & gmet(3,3)*gpk3**2 & 846 & +2.0_dp*(gpk1*gmet(1,2)*gpk2+ & 847 & gpk1*gmet(1,3)*gpk3+ & 848 & gpk2*gmet(2,3)*gpk3 ) ) 849 dkpg2=htpisq*2.0_dp*& 850 & (gpk1*(dgmetds(1,1)*gpk1+dgmetds(1,2)*gpk2+dgmetds(1,3)*gpk3)+ & 851 & gpk2*(dgmetds(2,1)*gpk1+dgmetds(2,2)*gpk2+dgmetds(2,3)*gpk3)+ & 852 & gpk3*(dgmetds(3,1)*gpk1+dgmetds(3,2)*gpk2+dgmetds(3,3)*gpk3) ) 853 dkinetic=dkpg2 854 if(kpg2>ecut-ecutsm)then 855 if(kpg2>ecut-tol12)then 856 ! The wavefunction has been filtered : no derivative 857 dkinetic=0.0_dp 858 else 859 xx=(ecut-kpg2)*ecutsm_inv 860 ! This kinetic cutoff smoothing function and its xx derivatives 861 ! were produced with Mathematica and the fortran code has been 862 ! numerically checked against Mathematica. 863 fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx)))) 864 dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2 865 ! d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*& 866 ! & (-144+45*xx))))))*fsm**3 867 dkinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm) 868 end if 869 end if 870 dkinpw(ig)=dkinetic/effmass_free 871 end do 872 873 end subroutine kpgstr
m_kg/mkkin [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
mkkin
FUNCTION
compute elements of kinetic energy operator in reciprocal space at a given k point
INPUTS
ecut=cut-off energy for plane wave basis sphere (Ha) ecutsm=smearing energy for plane wave kinetic energy (Ha) effmass_free=effective mass for electrons (1. in common case) gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$) idir1 = 1st direction of the derivative (if 1 <= idir1 <= 3, not used otherwise) idir2 = 2st direction of the derivative (if 1 <= idir1,idir2 <= 3, not used otherwise)) kg(3,npw)=integer coordinates of planewaves in basis sphere. kpt(3)=reduced coordinates of k point npw=number of plane waves at kpt.
OUTPUT
kinpw(npw)=(modified) kinetic energy (or derivative) for each plane wave (Hartree)
NOTES
Usually, the kinetic energy expression is $(1/2) (2 \pi)^2 (k+G)^2 $ However, the present implementation allows for a modification of this kinetic energy, in order to obtain smooth total energy curves with respect to the cut-off energy or the cell size and shape. Thus the usual expression is kept if it is lower then ecut-ecutsm, zero is returned beyond ecut, and in between, the kinetic energy is DIVIDED by a smearing factor (to make it infinite at the cut-off energy). The smearing factor is $x^2 (3-2x)$, where x = (ecut- unmodified energy)/ecutsm. This smearing factor is also used to derived a modified kinetic contribution to stress, in another routine (forstrnps.f) Also, in order to break slightly the symmetry between axes, that causes sometimes a degeneracy of eigenvalues and do not allow to obtain the same results on different machines, there is a modification by one part over 1.0e12 of the metric tensor elements (1,1) and (3,3)
SOURCE
323 subroutine mkkin (ecut,ecutsm,effmass_free,gmet,kg,kinpw,kpt,npw,idir1,idir2) 324 325 !Arguments ------------------------------------ 326 !scalars 327 integer,intent(in) :: npw 328 integer,intent(in) :: idir1,idir2 329 real(dp),intent(in) :: ecut,ecutsm,effmass_free 330 !arrays 331 integer,intent(in) :: kg(3,npw) 332 real(dp),intent(in) :: gmet(3,3),kpt(3) 333 real(dp),intent(out) :: kinpw(npw) 334 335 !Local variables------------------------------- 336 !scalars 337 integer :: ig,order 338 real(dp),parameter :: break_symm=1.0d-11 339 real(dp) :: ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq,kinetic,kpg2,dkpg2,xx 340 real(dp) :: d1kpg2,d2kpg2,ddfsm, dfsm 341 !arrays 342 real(dp) :: gmet_break(3,3) !, tsec(2) 343 344 ! ************************************************************************* 345 346 ! htpisq is (1/2) (2 Pi) **2: 347 htpisq=0.5_dp*(two_pi)**2 348 349 ecutsm_inv=0.0_dp 350 if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm 351 352 gmet_break(:,:)=gmet(:,:) 353 gmet_break(1,1)=(1.0_dp+break_symm)*gmet(1,1) 354 gmet_break(3,3)=(1.0_dp-break_symm)*gmet(3,3) 355 356 order=0 ! Compute the kinetic operator 357 if (idir1>0.and.idir1<4) then 358 order=1 ! Compute the 1st derivative of the kinetic operator 359 if (idir2>0.and.idir2<4) then 360 order=2 ! Compute the 2nd derivative of the kinetic operator 361 end if 362 end if 363 364 !$OMP PARALLEL DO PRIVATE(dkpg2,d1kpg2,d2kpg2,gpk1,gpk2,gpk3,ig,kinetic,kpg2,xx,fsm,dfsm,ddfsm) & 365 !$OMP SHARED(kinpw,ecut,ecutsm,ecutsm_inv) & 366 !$OMP SHARED(gmet_break,htpisq,idir1,idir2,kg,kpt,npw) 367 do ig=1,npw 368 gpk1=dble(kg(1,ig))+kpt(1) 369 gpk2=dble(kg(2,ig))+kpt(2) 370 gpk3=dble(kg(3,ig))+kpt(3) 371 kpg2=htpisq*& 372 & ( gmet_break(1,1)*gpk1**2+ & 373 & gmet_break(2,2)*gpk2**2+ & 374 & gmet_break(3,3)*gpk3**2 & 375 & +2.0_dp*(gpk1*gmet_break(1,2)*gpk2+ & 376 & gpk1*gmet_break(1,3)*gpk3+ & 377 & gpk2*gmet_break(2,3)*gpk3 ) ) 378 select case (order) 379 case(0) 380 kinetic=kpg2 381 case(1) 382 dkpg2=htpisq*2.0_dp*& 383 & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3) 384 kinetic=dkpg2 385 case(2) 386 dkpg2=htpisq*2.0_dp*gmet_break(idir1,idir2) 387 kinetic=dkpg2 388 end select 389 if(kpg2>ecut-ecutsm)then 390 if(kpg2>ecut-tol12)then 391 if(order==0) then 392 ! Will filter the wavefunction, based on this value, in cgwf.f, getghc.f and precon.f 393 kinetic=huge(zero)*1.d-10 394 else 395 ! The wavefunction has been filtered : no derivative 396 kinetic=0 397 end if 398 else 399 if(order==0) then 400 xx=max( (ecut-kpg2)*ecutsm_inv , 1.0d-20) 401 else 402 xx=(ecut-kpg2)*ecutsm_inv 403 end if 404 if(order==2) then 405 d1kpg2=htpisq*2.0_dp*& 406 & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3) 407 d2kpg2=htpisq*2.0_dp*& 408 & (gmet_break(idir2,1)*gpk1+gmet_break(idir2,2)*gpk2+gmet_break(idir2,3)*gpk3) 409 end if 410 ! This kinetic cutoff smoothing function and its xx derivatives 411 ! were produced with Mathematica and the fortran code has been 412 ! numerically checked against Mathematica. 413 fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx)))) 414 if(order>0) dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2 415 if(order>1) ddfsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*(-144+45*xx))))))*fsm**3 416 select case (order) 417 case(0) 418 kinetic=kpg2*fsm 419 case(1) 420 kinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm) 421 case(2) 422 kinetic=dkpg2*fsm& 423 & -2.0_dp*d1kpg2*dfsm*ecutsm_inv*d2kpg2& 424 & +kpg2*ddfsm*(ecutsm_inv**2)*d1kpg2*d2kpg2& 425 & -kpg2*dfsm*ecutsm_inv*dkpg2 426 end select 427 end if 428 end if 429 kinpw(ig)=kinetic/effmass_free 430 end do 431 !$OMP END PARALLEL DO 432 433 end subroutine mkkin
m_kg/mkkpg [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
mkkpg
FUNCTION
Compute all (k+G) vectors (dp, in reduced coordinates) for given k point, from integer coordinates of G vectors. Eventually compute related data.
INPUTS
kg(3,npw)=integer coords of planewaves in basis sphere kpt(3)=k point in terms of recip. translations nkpg=second dimension of array kpg npw=number of plane waves in reciprocal space
OUTPUT
kpg(npw,3)= (k+G) components
=== if nkpg==9 ===
kpg(npw,4:9)= [(k+G)_a].[(k+G)_b] quantities
SOURCE
899 subroutine mkkpg(kg, kpg, kpt, nkpg, npw) 900 901 !Arguments ------------------------------------ 902 !scalars 903 integer,intent(in) :: nkpg,npw 904 !arrays 905 integer,intent(in) :: kg(3,npw) 906 real(dp),intent(in) :: kpt(3) 907 real(dp),intent(out) :: kpg(npw,nkpg) 908 909 !Local variables------------------------------- 910 !scalars 911 integer :: ipw,mu,mua,mub 912 character(len=500) :: msg 913 !arrays 914 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 915 916 ! ************************************************************************* 917 918 if (nkpg==0) return 919 920 !-- Test nkpg -- 921 if (nkpg/=3.and.nkpg/=9) then 922 write(msg, '(a,i0)' )' Bad value for nkpg !',nkpg 923 ABI_BUG(msg) 924 end if 925 926 !-- Compute (k+G) -- 927 !$OMP PARALLEL DO COLLAPSE(2) & 928 !$OMP PRIVATE(mu,ipw) 929 do ipw=1,npw 930 do mu=1,3 931 kpg(ipw,mu)=kpt(mu)+dble(kg(mu,ipw)) 932 end do 933 end do 934 !$OMP END PARALLEL DO 935 936 !-- Compute [(k+G)_a].[(k+G)_b] -- 937 if (nkpg==9) then 938 !$OMP PARALLEL DO COLLAPSE(2) & 939 !$OMP PRIVATE(ipw,mu,mua,mub) 940 do ipw=1,npw 941 do mu=4,9 942 mua=alpha(mu-3);mub=beta(mu-3) 943 kpg(ipw,mu)=kpg(ipw,mua)*kpg(ipw,mub) 944 end do 945 end do 946 !$OMP END PARALLEL DO 947 end if 948 949 end subroutine mkkpg
m_kg/mkkpgcart [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
mkkpgcart
FUNCTION
Compute all (k+G) vectors (dp, in cartesian coordinates) for given k point, from integer coordinates of G vectors.
INPUTS
gprimd(3,3)=dimensional reciprocal space primitive translations kg(3,npw)=integer coords of planewaves in basis sphere kpt(3)=k point in terms of recip. translations nkpg=second dimension of array kpg npw=number of plane waves in reciprocal space
OUTPUT
kpg(npw,3)= (k+G) components
SOURCE
1124 subroutine mkkpgcart(gprimd,kg,kpgcar,kpt,nkpg,npw) 1125 1126 !Arguments ------------------------------------ 1127 !scalars 1128 integer,intent(in) :: nkpg,npw 1129 !arrays 1130 real(dp),intent(in) :: gprimd(3,3) 1131 integer,intent(in) :: kg(3,npw) 1132 real(dp),intent(in) :: kpt(3) 1133 real(dp),intent(out) :: kpgcar(npw,nkpg) 1134 1135 !Local variables------------------------------- 1136 !scalars 1137 integer :: ipw,mu 1138 character(len=500) :: msg 1139 !arrays 1140 real(dp),allocatable :: kpg(:,:) 1141 1142 ! ************************************************************************* 1143 1144 DBG_ENTER("COLL") 1145 1146 if (nkpg==0) return 1147 1148 !-- Test nkpg -- 1149 if (nkpg/=3) then 1150 write(msg, '(a,i0)' )' Bad value for nkpg !',nkpg 1151 ABI_BUG(msg) 1152 end if 1153 1154 !-- Compute (k+G) -- 1155 ABI_MALLOC(kpg,(npw,nkpg)) 1156 !$OMP PARALLEL DO COLLAPSE(2) & 1157 !$OMP PRIVATE(mu,ipw) 1158 do ipw=1,npw 1159 do mu=1,3 1160 kpg(ipw,mu)=kpt(mu)+dble(kg(mu,ipw)) 1161 end do 1162 end do 1163 !$OMP END PARALLEL DO 1164 1165 !$OMP PARALLEL DO & 1166 !$OMP PRIVATE(ipw) 1167 do ipw=1,npw 1168 kpgcar(ipw,1)=kpg(ipw,1)*gprimd(1,1)+kpg(ipw,2)*gprimd(1,2)+kpg(ipw,3)*gprimd(1,3) 1169 kpgcar(ipw,2)=kpg(ipw,1)*gprimd(2,1)+kpg(ipw,2)*gprimd(2,2)+kpg(ipw,3)*gprimd(2,3) 1170 kpgcar(ipw,3)=kpg(ipw,1)*gprimd(3,1)+kpg(ipw,2)*gprimd(3,2)+kpg(ipw,3)*gprimd(3,3) 1171 end do 1172 !$OMP END PARALLEL DO 1173 ABI_FREE(kpg) 1174 1175 DBG_EXIT("COLL") 1176 1177 end subroutine mkkpgcart
m_kg/ph1d3d [ Functions ]
[ Top ] [ m_kg ] [ Functions ]
NAME
ph1d3d
FUNCTION
Compute the three-dimensional phase factor $e^{i 2 \pi (k+G) cdot xred}$ from the three one-dimensional factors, the k point coordinates, and the atom positions, for all planewaves which fit in the fft box.
INPUTS
iatom, jatom= bounds of atom indices in ph1d for which ph3d has to be computed kg_k(3,npw_k)=reduced planewave coordinates. matblk= dimension of ph3d natom= dimension of ph1d npw=number of plane waves n1,n2,n3=dimensions of fft box (ngfft(3)). phkxred(2,natom)=phase factors exp(2 pi k.xred) ph1d(2,(2*n1+1)*natom+(2*n2+1)*natom+(2*n3+1)*natom)=exp(2Pi i G xred) for vectors (Gx,0,0), (0,Gy,0) and (0,0,Gz) with components ranging from -nj <= Gj <= nj
OUTPUT
ph3d(2,npw_k,matblk)=$e^{2 i \pi (k+G) cdot xred}$ for vectors (Gx,Gy,Gz),
and for atoms in the range iatom to jatom with respect to ph1d
SOURCE
612 subroutine ph1d3d(iatom, jatom, kg_k, matblk, natom, npw_k, n1, n2, n3, phkxred, ph1d, ph3d) 613 614 !Arguments ------------------------------------ 615 !scalars 616 integer,intent(in) :: iatom,jatom,matblk,n1,n2,n3,natom,npw_k 617 !arrays 618 integer,intent(in) :: kg_k(3,npw_k) 619 real(dp),intent(in) :: ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom) 620 real(dp),intent(in) :: phkxred(2,natom) 621 real(dp),intent(out) :: ph3d(2,npw_k,matblk) 622 623 !Local variables------------------------------- 624 !scalars 625 integer :: i1,ia,iatblk,ig,kg1,kg2,kg3,shift1,shift2,shift3 626 real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,phkxi,phkxr 627 character(len=500) :: msg 628 !arrays 629 real(dp),allocatable :: ph1kxred(:,:) 630 631 ! ************************************************************************* 632 633 if(matblk-1 < jatom-iatom)then 634 write(msg,'(a,a,a,a,a,i0,a,a,i0,a,i0,a)')& 635 'Input natom-1 must be larger or equal to jatom-iatom,',ch10,& 636 'while their value is : ',ch10,& 637 'natom-1 = ',natom-1,ch10,& 638 'jatom=',jatom,', iatom=',iatom,'.' 639 ABI_BUG(msg) 640 end if 641 642 ABI_MALLOC(ph1kxred,(2,-n1:n1)) 643 644 ! ia runs from iatom to jatom 645 do ia=iatom,jatom 646 647 ! iatblk runs from 1 to matblk 648 iatblk=ia-iatom+1 649 shift1=1+n1+(ia-1)*(2*n1+1) 650 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 651 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 652 ! Compute product of phkxred by phase for the first component of G vector 653 phkxr=phkxred(1,ia) 654 phkxi=phkxred(2,ia) 655 do i1=-n1,n1 656 ph1kxred(1,i1)=ph1d(1,i1+shift1)*phkxr-ph1d(2,i1+shift1)*phkxi 657 ph1kxred(2,i1)=ph1d(2,i1+shift1)*phkxr+ph1d(1,i1+shift1)*phkxi 658 end do 659 660 ! Compute tri-dimensional phase factor 661 !$OMP PARALLEL DO PRIVATE(ig,ph1r,ph1i,ph2r,ph2i,ph3r,ph3i,ph12r,ph12i,kg1,kg2,kg3) 662 do ig=1,npw_k 663 kg1=kg_k(1,ig) 664 kg2=kg_k(2,ig)+shift2 665 kg3=kg_k(3,ig)+shift3 666 ph1r=ph1kxred(1,kg1) 667 ph1i=ph1kxred(2,kg1) 668 ph2r=ph1d(1,kg2) 669 ph2i=ph1d(2,kg2) 670 ph3r=ph1d(1,kg3) 671 ph3i=ph1d(2,kg3) 672 ph12r=ph1r*ph2r-ph1i*ph2i 673 ph12i=ph1r*ph2i+ph1i*ph2r 674 ph3d(1,ig,iatblk)=ph12r*ph3r-ph12i*ph3i 675 ph3d(2,ig,iatblk)=ph12r*ph3i+ph12i*ph3r 676 end do 677 !$OMP END PARALLEL DO 678 end do 679 680 ABI_FREE(ph1kxred) 681 682 end subroutine ph1d3d