TABLE OF CONTENTS
ABINIT/corrmetalwf1 [ Functions ]
NAME
corrmetalwf1
FUNCTION
Response function calculation only: Correct 1st-order wave-function, taking into account "metallic" occupations. 1st-order WF orthogonal to C_n,k+q, restore the "active space" content of the first-order WF. receives a single band at k as input, and works on all bands at k+q
INPUTS
cgq(2,mcgq)=planewave coefficients of wavefunctions at k+q
cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
cwaveprj(natom,nspinor)= 1st-order wave-function before correction
projected on NL projectors (PAW)
cycle_bands(nband)=array of logicals for bands we have on this cpu
eig1(2*nband**2)=first-order eigenvalues (hartree)
fermie1=derivative of fermi energy wrt (strain) perturbation
ghc(2,npw1*nspinor)=<G|H0-eig0_k.I|C1 band,k> (NCPP) or <G|H0-eig0_k.S0|C1 band,k> (PAW)
(C1 before correction)
iband=index of current band
ibgq=shift to be applied on the location of data in the array cprjq
icgq=shift to be applied on the location of data in the array cgq
istwf_k=option parameter that describes the storage of wfs
mcgq=second dimension of the cgq array
mcprjq=second dimension of the cprjq array
mpi_enreg=information about MPI parallelization
natom=number of atoms in cell
nband=number of bands
npw1=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
occ(nband)=occupation number for each band for each k.
rocceig(nband,nband)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
if this ratio has been attributed to the band n (second argument), zero otherwise
timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
usepaw=flag for PAW
OUTPUT
cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
projected on NL projectors (PAW)
edocc(nband)=correction to 2nd-order total energy coming from changes of occupations
wf_corrected=flag put to 1 if input cwave1 is effectively different from output cwavef
NOTES
Was part of dfpt_vtowfk before.
SOURCE
932 subroutine corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,edocc,eig1,fermie1,ghc,iband, & 933 & ibgq,icgq,istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,nspinor,occ,rocceig,timcount,& 934 & usepaw,wf_corrected) 935 936 !Arguments ------------------------------------ 937 !scalars 938 integer,intent(in) :: iband,ibgq,icgq,istwf_k,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw 939 integer,intent(out) :: wf_corrected 940 real(dp),intent(in) :: fermie1 941 type(MPI_type),intent(in) :: mpi_enreg 942 !arrays 943 logical,intent(in) :: cycle_bands(nband) 944 real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor) 945 real(dp),intent(in) :: eig1(2*nband**2),ghc(2,npw1*nspinor),occ(nband),rocceig(nband,nband) 946 real(dp),intent(out) :: cwave1(2,npw1*nspinor),edocc(nband) 947 type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw) 948 type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw) !vz_i 949 950 !Local variables------------------------------- 951 !scalars 952 integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii 953 integer :: ibandkq_me, ierr, iband_ 954 integer :: wf_corrected_ 955 real(dp) :: facti,factr,invocc 956 real(dp) :: edocc_tmp 957 !arrays 958 integer :: bands_treated_now(nband) 959 real(dp) :: tsec(2) 960 real(dp),allocatable :: cwcorr(:,:) 961 type(pawcprj_type) :: cwaveprj1_corr(natom,nspinor*usepaw) 962 963 ! ********************************************************************* 964 965 DBG_ENTER("COLL") 966 967 call timab(214+timcount,1,tsec) 968 969 bands_treated_now = 0 970 bands_treated_now(iband) = 1 971 call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr) 972 973 !At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf). 974 !Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 . 975 976 ABI_MALLOC(cwcorr,(2,npw1*nspinor)) 977 978 wf_corrected=0 979 980 if(usepaw==1) then 981 call pawcprj_alloc(cwaveprj1_corr, cwaveprj1(1,1)%ncpgr, cwaveprj1(:,1)%nlmn) 982 end if 983 984 ! loop iband_ over all bands being treated for the moment 985 ! all procs in pool of bands should be working on the same iband_ at a given time 986 ! I will save in _my_ array cwave1, if iband==iband_ 987 do iband_ = 1, nband 988 if (bands_treated_now(iband_) == 0) cycle 989 990 cwcorr = zero 991 if (usepaw==1) then 992 call pawcprj_set_zero(cwaveprj1_corr) 993 end if 994 edocc_tmp = zero 995 wf_corrected_=0 996 997 !Correct WF only for occupied states 998 if (abs(occ(iband_)) > tol8) then 999 invocc=one/occ(iband_) 1000 1001 1002 ! Loop over WF at k+q subspace 1003 ibandkq_me = 0 1004 do ibandkq=1,nband 1005 if(cycle_bands(ibandkq)) cycle 1006 ibandkq_me = ibandkq_me + 1 1007 1008 ! Select bands with variable occupation 1009 if (abs(rocceig(ibandkq,iband_))>tol8) then 1010 1011 wf_corrected_=1 1012 1013 index_eig1=2*ibandkq-1+(iband_-1)*2*nband 1014 index_cgq=npw1*nspinor*(ibandkq_me-1)+icgq 1015 1016 if(ibandkq==iband_) then 1017 factr=rocceig(ibandkq,iband_)*invocc*(eig1(index_eig1)-fermie1) 1018 else 1019 factr=rocceig(ibandkq,iband_)*invocc*eig1(index_eig1) 1020 end if 1021 facti = rocceig(ibandkq,iband_)*invocc*eig1(index_eig1+1) 1022 1023 ! Apply correction to 1st-order WF 1024 !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwcorr,facti,factr,index_cgq,npw1,nspinor) 1025 do ii=1,npw1*nspinor 1026 cwcorr(1,ii)=cwcorr(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq)) 1027 cwcorr(2,ii)=cwcorr(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq)) 1028 end do 1029 1030 ! In the PAW case, also apply correction to projected WF 1031 if (usepaw==1) then 1032 index_cprjq=nspinor*(ibandkq_me-1)+ibgq 1033 call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1_corr) 1034 end if 1035 1036 ! The factor of two is needed because we compute the 2DTE, and not E(2) 1037 edocc_tmp = edocc_tmp-two*(factr*eig1(index_eig1)+facti*eig1(index_eig1+1)) 1038 1039 end if ! Variable occupations 1040 end do ! Loop over k+q subspace 1041 end if ! if occupied states 1042 1043 ! 3) reduce over bands to get all contributions to correction 1044 ! need MPI reduce over band communicator only 1045 call xmpi_sum(cwcorr, mpi_enreg%comm_band, ierr) 1046 if (usepaw==1) then 1047 call pawcprj_mpi_sum(cwaveprj1_corr, mpi_enreg%comm_band, ierr) 1048 end if 1049 1050 ! this sums over the k+q contributions to the present iband_ 1051 call xmpi_sum(edocc_tmp, mpi_enreg%comm_band, ierr) 1052 call xmpi_sum(wf_corrected_, mpi_enreg%comm_band, ierr) 1053 1054 1055 ! 4) add correction to the cwave1 1056 ! if I have iband_, correct my cwave1 1057 if (iband_==iband) then 1058 if (wf_corrected_ > 0) wf_corrected = 1 1059 edocc(iband) = edocc_tmp 1060 cwave1 = cwcorr 1061 call cg_zaxpy(npw1*nspinor,(/one,zero/),cwavef,cwave1) 1062 !Idem for cprj 1063 if (usepaw==1) then 1064 call pawcprj_copy(cwaveprj,cwaveprj1) 1065 call pawcprj_zaxpby((/one,zero/),(/one,zero/),cwaveprj1_corr,cwaveprj1) 1066 end if 1067 end if 1068 end do ! loop over all bands presently running in parallel 1069 1070 1071 !In the PAW case, compute <Psi^(1)_ortho|H-Eig0_k.S|Psi^(1)_parallel> contribution to 2DTE 1072 if (usepaw==1.and.wf_corrected==1) then 1073 !$OMP WORKSHARE 1074 cwcorr(:,:)=cwave1(:,:)-cwavef(:,:) 1075 !$OMP END WORKSHARE 1076 call dotprod_g(factr,facti,istwf_k,npw1*nspinor,1,cwcorr,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1077 edocc(iband)=edocc(iband)+four*factr 1078 end if 1079 1080 1081 ABI_FREE(cwcorr) 1082 if (usepaw==1) then 1083 call pawcprj_free(cwaveprj1_corr) 1084 end if 1085 1086 call timab(214+timcount,2,tsec) 1087 1088 DBG_EXIT("COLL") 1089 1090 end subroutine corrmetalwf1
ABINIT/dfpt_vtowfk [ Functions ]
NAME
dfpt_vtowfk
FUNCTION
This routine compute the partial density at a given k-point, for a given spin-polarization, from a fixed potential (vlocal1).
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions
cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
cplex=1 if rhoaug1 is real, 2 if rhoaug1 is complex
cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
projected with non-local projectors: cprj=<p_i|Cnk>
cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q>
dim_eig2rf = dimension for the second order eigenvalues
dtfil <type(datafiles_type)>=variables related to files
dtset <type(dataset_type)>=all input variables for this dataset
eig0_k(nband_k)=GS eigenvalues at k (hartree)
eig0_kq(nband_k)=GS eigenvalues at k+Q (hartree)
fermie1=derivative of fermi energy wrt (strain) perturbation
grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
ibg=shift to be applied on the location of data in the array cprj
ibgq=shift to be applied on the location of data in the array cprjq
ibg1=shift to be applied on the location of data in the array cprj1
icg=shift to be applied on the location of data in the array cg
icgq=shift to be applied on the location of data in the array cgq
icg1=shift to be applied on the location of data in the array cg1
idir=direction of the current perturbation
ikpt=k-point index number
ipert=type of the perturbation
isppol=1 index of current spin component
mband=maximum number of bands
mband_mem=maximum number of bands on this cpu
mcgq=second dimension of the cgq array
mcprjq=second dimension of the cprjq array
mkmem =number of k points trated by this node (GS data).
mk1mem =number of k points treated by this node (RF data)
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw or wfs at k
mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
natom=number of atoms in cell.
nband_k=number of bands at this k point for that spin polarization
ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
nnsclo_now=number of non-self-consistent loops for the current vtrial
(often 1 for SCF calculation, =nstep for non-SCF calculations)
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
n4,n5,n6 used for dimensioning real space arrays
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
prtvol=control print volume and debugging output
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
rf_hamk_dir2 <type(rf_hamiltonian_type)>= (used only when ipert=natom+11, so q=0)
same as rf_hamkq, but the direction of the perturbation is different
rhoaug1(cplex*n4,n5,n6,nspden)= density in electrons/bohr**3,
on the augmented fft grid. (cumulative, so input as well as output)
rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
if this ratio has been attributed to the band n (second argument), zero otherwise
ddk<wfk_t>=struct info for DDK file.
wtk_k=weight assigned to the k point.
OUTPUT
cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF
wavefunctions at k,q. They are orthogonalized to the occupied states.
cg1_active(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)=pw coefficients of RF
wavefunctions at k,q. They are orthogonalized to the active. Only needed for ieigrf/=0
edocc_k(nband_k)=correction to 2nd-order total energy coming
from changes of occupation
eeig0_k(nband_k)=zero-order eigenvalues contribution to 2nd-order total
energy from all bands at this k point.
eig1_k(2*nband_k**2)=first-order eigenvalues (hartree)
ek0_k(nband_k)=0-order kinetic energy contribution to 2nd-order total
energy from all bands at this k point.
ek1_k(nband_k)=1st-order kinetic energy contribution to 2nd-order total
energy from all bands at this k point.
eloc0_k(nband_k)=zero-order local contribution to 2nd-order total energy
from all bands at this k point.
end0_k(nband_k)=0-order nuclear dipole energy contribution to 2nd-order total
energy from all bands at this k point.
end1_k(nband_k)=1st-order nuclear dipole energy contribution to 2nd-order total
energy from all bands at this k point.
enl0_k(nband_k)=zero-order non-local contribution to 2nd-order total energy
from all bands at this k point.
enl1_k(nband_k)=first-order non-local contribution to 2nd-order total energy
from all bands at this k point.
evxctau0_k(nband_k)=0-order vxctau energy contribution to 2nd-order total
energy from all bands at this k point.
evxctau1_k(nband_k)=1-order vxctau energy contribution to 2nd-order total
energy from all bands at this k point.
gh1c_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
gh0c1_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}k+q-eig^{(0)}nk|\Psi^{(1)}kq>
The wavefunction is orthogonal to the active space (for metals). It is not coherent with cg1.
resid_k(nband_k)=residuals for each band over all k points,
rhoaug1(cplex*n4,n5,n6,nspden)= density in electrons/bohr**3,
on the augmented fft grid. (cumulative, so input as well as output).
==== if (gs_hamkq%usepaw==1) ====
cprj1(natom,nspinor*mband*mk1mem*nsppol*usecprj)=
1st-order wave functions at k,q projected with non-local projectors:
cprj1=<p_i|C1nk,q> where p_i is a non-local projector
pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
(cumulative, so input as well as output)
SOURCE
169 subroutine dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,& 170 & dim_eig2rf,dtfil,dtset,& 171 & edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,& 172 & ek0_k,ek1_k,eloc0_k,end0_k,end1_k,enl0_k,enl1_k,evxctau0_k,evxctau1_k,& 173 & fermie1,ffnl1,ffnl1_test,gh0c1_set,gh1c_set,grad_berry,gs_hamkq,& 174 & ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,& 175 & isppol,mband,mband_mem,mcgq,mcprjq,mkmem,mk1mem,& 176 & mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,& 177 & nnsclo_now,npw_k,npw1_k,nspinor,nsppol,& 178 & n4,n5,n6,occ_k,pawrhoij1,prtvol,psps,resid_k,rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,& 179 & ddk_f,wtk_k,nlines_done,cg1_out) 180 181 !Arguments ------------------------------------ 182 !scalars 183 integer,intent(in) :: cplex,dim_eig2rf,ibg 184 integer,intent(in) :: ibg1,ibgq,icg,icg1,icgq,idir,ikpt,ipert,isppol 185 integer,intent(in) :: mband,mcgq,mcprjq,mk1mem,mkmem 186 integer,intent(in) :: mband_mem 187 integer,intent(in) :: mpw,mpw1,n4,n5,n6,natom,ncpgr 188 integer,intent(in) :: nnsclo_now,nspinor,nsppol,prtvol 189 integer,optional,intent(in) :: cg1_out 190 integer,intent(in) :: nband_k,npw1_k,npw_k 191 integer,intent(inout) :: nlines_done 192 real(dp),intent(in) :: fermie1,wtk_k 193 type(MPI_type),intent(in) :: mpi_enreg 194 type(datafiles_type),intent(in) :: dtfil 195 type(dataset_type),intent(in) :: dtset 196 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 197 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq,rf_hamk_dir2 198 type(pseudopotential_type),intent(in) :: psps 199 !arrays 200 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol),cgq(2,mcgq) 201 real(dp),intent(in) :: eig0_k(nband_k),eig0_kq(nband_k) 202 real(dp),intent(in) :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:) 203 real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband_k) 204 real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k) 205 real(dp),intent(inout) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) 206 real(dp),intent(inout) :: rhoaug1(cplex*n4,n5,n6,gs_hamkq%nvloc) 207 real(dp),intent(inout) :: cg1_active(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf) 208 real(dp),intent(inout) :: gh1c_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf) 209 real(dp),intent(inout) :: gh0c1_set(2,mpw1*nspinor*mband_mem*mk1mem*nsppol*dim_eig2rf) 210 real(dp),intent(inout) :: edocc_k(nband_k),eeig0_k(nband_k),eig1_k(2*nband_k**2) 211 real(dp),intent(out) :: ek0_k(nband_k),eloc0_k(nband_k) 212 real(dp),intent(inout) :: ek1_k(nband_k) 213 real(dp),intent(out) :: end0_k(nband_k),end1_k(nband_k),enl0_k(nband_k),enl1_k(nband_k) 214 real(dp),intent(out) :: evxctau0_k(nband_k),evxctau1_k(nband_k) 215 real(dp),intent(out) :: resid_k(nband_k) 216 !TODO: PAW distrib bands mband_mem 217 type(pawcprj_type),intent(in) :: cprj(natom,nspinor*mband_mem*mkmem*nsppol*gs_hamkq%usecprj) 218 type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq) 219 type(pawcprj_type),intent(inout) :: cprj1(natom,nspinor*mband_mem*mk1mem*nsppol*gs_hamkq%usecprj) 220 type(pawrhoij_type),intent(inout) :: pawrhoij1(natom*gs_hamkq%usepaw) 221 type(wfk_t),intent(inout) :: ddk_f(4) 222 223 !Local variables------------------------------- 224 !scalars 225 integer,parameter :: level=14,tim_fourwf=5 226 integer,save :: nskip=0 227 integer :: iband,idir0,ierr,igs,igscq,ii,dim_dcwf,inonsc 228 integer :: iband_me,nband_me !, unit_me 229 integer :: iorder_cprj,iorder_cprj1,ipw,iscf_mod,ispinor,me,mgscq,nkpt_max 230 integer :: option,opt_gvnlx1,quit,test_ddk 231 integer :: tocceig,usedcwavef,ptr,shift_band 232 real(dp) :: aa,ai,ar,eig0nk,resid,residk,scprod,energy_factor 233 character(len=500) :: message 234 type(rf2_t) :: rf2 235 !arrays 236 logical,allocatable :: cycle_bands(:) 237 integer :: rank_band(nband_k), bands_treated_now(nband_k) 238 real(dp) :: tsec(2) 239 real(dp),allocatable :: cwave0(:,:),cwave1(:,:),cwavef(:,:) 240 real(dp),allocatable :: dcwavef(:,:),gh1c_n(:,:),gh0c1(:,:),ghc_vectornd(:,:) 241 real(dp),allocatable :: ghc_vxctau(:,:) 242 real(dp),allocatable :: gsc(:,:),gscq(:,:),gvnlx1(:,:),gvnlxc(:,:) 243 real(dp),pointer :: kinpw1(:) 244 type(pawcprj_type),allocatable :: cwaveprj(:,:),cwaveprj0(:,:),cwaveprj1(:,:) 245 246 ! ********************************************************************* 247 248 DBG_ENTER('COLL') 249 250 !Keep track of total time spent in dfpt_vtowfk 251 call timab(128,1,tsec) 252 253 nkpt_max=50; if (xmpi_paral==1) nkpt_max=-1 254 255 if(prtvol>2 .or. ikpt<=nkpt_max)then 256 write(message,'(2a,i5,2x,a,3f9.5,2x,a)')ch10,' Non-SCF iterations; k pt #',ikpt,'k=',gs_hamkq%kpt_k(:),'band residuals:' 257 call wrtout(std_out,message) 258 end if 259 260 !Initializations and allocations 261 me=mpi_enreg%me_kpt 262 quit=0 263 264 !The value of iscf must be modified if ddk perturbation 265 iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3 266 267 kinpw1 => gs_hamkq%kinpw_kp 268 ABI_MALLOC(gh0c1,(2,npw1_k*nspinor)) 269 ABI_MALLOC(gvnlxc,(2,npw1_k*nspinor)) 270 ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor)) 271 ABI_MALLOC(cwave0,(2,npw_k*nspinor)) 272 ABI_MALLOC(cwavef,(2,npw1_k*nspinor)) 273 ABI_MALLOC(cwave1,(2,npw1_k*nspinor)) 274 ABI_MALLOC(gh1c_n,(2,npw1_k*nspinor)) 275 if (gs_hamkq%usepaw==1) then 276 ABI_MALLOC(gsc,(2,npw1_k*nspinor)) 277 else 278 ABI_MALLOC(gsc,(0,0)) 279 end if 280 281 !Read the npw and kg records of wf files 282 test_ddk=0 283 if ((ipert==natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and.& 284 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 285 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.& 286 & ipert==natom+10.or.ipert==natom+11) then 287 test_ddk=1 288 if(ipert==natom+10.or.ipert==natom+11) test_ddk=0 289 end if 290 291 !Additional stuff for PAW 292 ABI_MALLOC(cwaveprj0,(0,0)) 293 if (gs_hamkq%usepaw==1) then 294 ! 1-Compute all <g|S|Cnk+q> 295 igscq=0 296 !TODO MJV: PAW mband_mem 297 mgscq=mpw1*nspinor*mband_mem 298 ABI_MALLOC_OR_DIE(gscq,(2,mgscq), ierr) 299 300 call getgsc(cgq,cprjq,gs_hamkq,gscq,ibgq,icgq,igscq,ikpt,isppol,mcgq,mcprjq,& 301 & mgscq,mpi_enreg,natom,nband_k,npw1_k,dtset%nspinor,select_k=KPRIME_H_KPRIME) 302 ! 2-Initialize additional scalars/arrays 303 iorder_cprj=0;iorder_cprj1=0 304 dim_dcwf=npw1_k*nspinor;if (ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) dim_dcwf=0 305 ABI_MALLOC(dcwavef,(2,dim_dcwf)) 306 if (gs_hamkq%usecprj==1) then 307 ABI_FREE(cwaveprj0) 308 ABI_MALLOC(cwaveprj0,(natom,nspinor)) 309 call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj) 310 end if 311 ABI_MALLOC(cwaveprj,(natom,nspinor)) 312 ABI_MALLOC(cwaveprj1,(natom,nspinor)) 313 call pawcprj_alloc(cwaveprj ,0,gs_hamkq%dimcprj) 314 call pawcprj_alloc(cwaveprj1,0,gs_hamkq%dimcprj) 315 else 316 igscq=0;mgscq=0;dim_dcwf=0 317 ABI_MALLOC(gscq,(0,0)) 318 ABI_MALLOC(dcwavef,(0,0)) 319 ABI_MALLOC(cwaveprj,(0,0)) 320 ABI_MALLOC(cwaveprj1,(0,0)) 321 end if 322 323 energy_factor=two 324 if(ipert==natom+10.or.ipert==natom+11) energy_factor=six 325 326 !For rf2 perturbation : 327 if(ipert==natom+10.or.ipert==natom+11) then 328 call rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,ffnl1,ffnl1_test,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,& 329 mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f) 330 end if 331 332 call timab(139,1,tsec) 333 334 !====================================================================== 335 !================== LOOP OVER BANDS ================================== 336 !====================================================================== 337 338 call proc_distrb_band(rank_band,mpi_enreg%proc_distrb,ikpt,isppol,mband,& 339 & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band) 340 341 iband_me = 0 342 do iband=1,nband_k 343 344 ! Skip bands not treated by current proc 345 if( (mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me)) cycle 346 iband_me = iband_me + 1 347 348 !unit_me = 300+iband 349 !Get ground-state wavefunctions 350 ptr = 1+(iband_me-1)*npw_k*nspinor+icg 351 call cg_zcopy(npw_k*nspinor,cg(1,ptr),cwave0) 352 353 ! Get PAW ground state projected WF (cprj) 354 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.ipert/=natom+10.and.ipert/=natom+11) then 355 idir0 = idir 356 if(ipert==natom+3.or.ipert==natom+4) idir0 =1 357 ! PAW distributes cprj by band and k 358 call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,natom,iband_me,ibg,ikpt,iorder_cprj,& 359 & isppol,mband_mem,mkmem,natom,1,nband_me,nspinor,nsppol,dtfil%unpaw,& 360 !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,& 361 & icpgr=idir0,ncpgr=ncpgr) 362 end if 363 364 ! Get first-order wavefunctions 365 ptr = 1+(iband_me-1)*npw1_k*nspinor+icg1 366 call cg_zcopy(npw1_k*nspinor,cg1(1,ptr),cwavef) 367 368 ! Read PAW projected 1st-order WF (cprj) 369 ! Unuseful for the time being (will be recomputed in dfpt_cgwf) 370 ! if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 371 !TODO MJV: PAW 372 ! call pawcprj_get(gs_hamkq%atindx1,cwaveprj,cprj1,natom,iband,ibg1,ikpt,iorder_cprj1,& 373 ! & isppol,mband,mk1mem,natom,1,nband_k,nspinor,nsppol,dtfil%unpaw1, 374 ! & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 375 ! end if 376 377 ! Filter the wavefunctions for large modified kinetic energy 378 ! The GS wavefunctions should already be non-zero 379 do ispinor=1,nspinor 380 igs=(ispinor-1)*npw1_k 381 do ipw=1+igs,npw1_k+igs 382 if(kinpw1(ipw-igs)>huge(zero)*1.d-20)then 383 cwavef(1,ipw)=zero 384 cwavef(2,ipw)=zero 385 end if 386 end do 387 end do 388 389 390 ! If electric field, the derivative of the wf should be read, and multiplied by i. 391 if(test_ddk==1) then 392 ii = ddk_f(1)%findk(gs_hamkq%kpt_k) 393 ABI_CHECK(ii == ikpt, "ii != ikpt, something is wrong with k-point, check kptopt/ngkpt, etc") 394 !TODO MJV: check if this iband should be _me 395 call ddk_f(1)%read_bks(iband, ikpt, isppol, xmpio_single, cg_bks=gvnlx1) 396 397 ! Multiplication by -i 398 ! MVeithen 021212 : use + i instead, 399 ! See X. Gonze, Phys. Rev. B 55, 10337 (1997) [[cite:Gonze1997]] Eq. (79) 400 ! the operator used to compute the first-order derivative 401 ! of the wavefunctions with respect to an electric field 402 ! is $+i \frac{d}{dk}$ 403 ! This change will affect the computation of the 2dtes from non 404 ! stationary expressions, see dfpt_nstdy.f and dfpt_nstwf.f 405 do ipw=1,npw1_k*nspinor 406 ! aa=gvnlx1(1,ipw) 407 ! gvnlx1(1,ipw)=gvnlx1(2,ipw) 408 ! gvnlx1(2,ipw)=-aa 409 aa=gvnlx1(1,ipw) 410 gvnlx1(1,ipw)=-gvnlx1(2,ipw) 411 gvnlx1(2,ipw)=aa 412 end do 413 end if 414 415 ! Unlike in GS calculations, the inonsc loop is inside the band loop 416 ! nnsclo_now=number of non-self-consistent loops for the current vtrial 417 ! (often 1 for SCF calculation, =nstep for non-SCF calculations) 418 do inonsc=1,nnsclo_now 419 420 ! Note that the following translation occurs in the called routine : 421 ! iband->band, nband_k->nband, npw_k->npw, npw1_k->npw1 422 eig0nk=eig0_k(iband) 423 usedcwavef=gs_hamkq%usepaw;if (dim_dcwf==0) usedcwavef=0 424 if (inonsc==1) usedcwavef=2*usedcwavef 425 opt_gvnlx1=0;if (ipert==natom+2) opt_gvnlx1=1 426 if (ipert==natom+2.and.gs_hamkq%usepaw==1.and.inonsc==1) opt_gvnlx1=2 427 428 if ( (ipert/=natom+10 .and. ipert/=natom+11) .or. abs(occ_k(iband))>tol8 ) then 429 nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 430 431 bands_treated_now = 0 432 bands_treated_now(iband) = 1 433 call xmpi_sum(bands_treated_now,mpi_enreg%comm_band,ierr) 434 435 if (dtset%rf2_dkdk==2 .and. (idir==1 .or. idir==2 .or. idir==3)) then 436 eig1_k = zero 437 resid = zero 438 else 439 call dfpt_cgwf(iband,iband_me,rank_band,bands_treated_now,dtset%berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,& 440 & rf2,dcwavef,& 441 & eig0_k,eig0_kq,eig1_k,gh0c1,gh1c_n,grad_berry,gsc,gscq,gs_hamkq,gvnlxc,gvnlx1,icgq,& 442 & idir,ipert,igscq,mcgq,mgscq,mpi_enreg,mpw1,natom,nband_k,nband_me,dtset%nbdbuf,dtset%nline,& 443 & npw_k,npw1_k,nspinor,opt_gvnlx1,prtvol,quit,resid,rf_hamkq,dtset%dfpt_sciss,dtset%tolrde,& 444 & dtset%tolwfr,usedcwavef,dtset%wfoptalg,nlines_done) 445 end if 446 resid_k(iband)=resid 447 448 else 449 resid_k(iband)=zero 450 end if 451 452 if (ipert/=natom+10 .and. ipert/= natom+11) then 453 ! At this stage, the 1st order function cwavef is orthogonal to cgq (unlike 454 ! when it is input to dfpt_cgwf). Here, restore the "active space" content 455 ! of the first-order wavefunction, to give cwave1. 456 ! PAW: note that dcwavef (1st-order change of WF due to overlap change) 457 ! remains in the subspace orthogonal to cgq 458 call proc_distrb_cycle_bands(cycle_bands, mpi_enreg%proc_distrb,ikpt,isppol,me) 459 if (dtset%prtfull1wf>0) then 460 call full_active_wf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,eig1_k,fermie1,& 461 & eig0nk,eig0_kq,dtset%elph2_imagden,iband,ibgq,icgq,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,& 462 & 0,gs_hamkq%usepaw) 463 edocc_k=zero 464 tocceig=1 465 else 466 call corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,edocc_k,eig1_k,fermie1,gh0c1,& 467 & iband,ibgq,icgq,gs_hamkq%istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,& 468 & occ_k,rocceig,0,gs_hamkq%usepaw,tocceig) 469 end if 470 ABI_FREE (cycle_bands) 471 else 472 tocceig=0 473 call cg_zcopy(npw1_k*nspinor,cwavef,cwave1) 474 if (gs_hamkq%usepaw==1) then 475 call pawcprj_copy(cwaveprj,cwaveprj1) 476 end if 477 end if 478 479 if (abs(occ_k(iband))<= tol8) then 480 ek0_k(iband)=zero 481 ek1_k(iband)=zero 482 eeig0_k(iband)=zero 483 end0_k(iband)=zero 484 end1_k(iband)=zero 485 enl0_k(iband)=zero 486 enl1_k(iband)=zero 487 eloc0_k(iband)=zero 488 evxctau0_k(iband)=zero 489 evxctau1_k(iband)=zero 490 nskip=nskip+1 491 else 492 ! Compute the 0-order kinetic operator contribution (with cwavef) 493 call meanvalue_g(ar,kinpw1,0,gs_hamkq%istwf_k,mpi_enreg,npw1_k,nspinor,cwavef,cwavef,0) 494 ! There is an additional factor of 2 with respect to the bare matrix element 495 ek0_k(iband)=energy_factor*ar 496 ! Compute the 1-order kinetic operator contribution (with cwave1 and cwave0), if needed. 497 ! Note that this is called only for ddk or strain, so that npw1_k=npw_k 498 if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4)then 499 call matrixelmt_g(ai,ar,rf_hamkq%dkinpw_k,gs_hamkq%istwf_k,0,npw_k,nspinor,cwave1,cwave0,& 500 & mpi_enreg%me_g0, mpi_enreg%comm_fft) 501 ! There is an additional factor of 4 with respect to the bare matrix element 502 ek1_k(iband)=two*energy_factor*ar 503 end if 504 505 ! Compute the 0-order nuclear dipole contribution (with cwavef) 506 ! only relevant for DDK 507 if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(gs_hamkq%vectornd)) ) then 508 ABI_MALLOC(ghc_vectornd,(2,npw_k*nspinor)) 509 ! ndat hard-coded as 1 510 call getghc_nucdip(cwavef,ghc_vectornd,gs_hamkq%gbound_k,gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kpt_k,& 511 & gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw_k,gs_hamkq%nvloc,& 512 & gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,nspinor,gs_hamkq%vectornd,gs_hamkq%gpu_option) 513 ! There is an additional factor of 2 with respect to the bare matrix element 514 end0_k(iband)=energy_factor*(DOT_PRODUCT(cwavef(1,1:npw_k*nspinor),ghc_vectornd(1,1:npw_k*nspinor))+& 515 & DOT_PRODUCT(cwavef(2,1:npw_k*nspinor),ghc_vectornd(2,1:npw_k*nspinor))) 516 ABI_FREE(ghc_vectornd) 517 end if 518 519 ! Compute the 0-order vxctau contribution (with cwavef) 520 ! only relevant for DDK 521 if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(gs_hamkq%vxctaulocal)) ) then 522 ABI_MALLOC(ghc_vxctau,(2,npw_k*nspinor)) 523 ! ndat hard-coded as 1 524 call getghc_mGGA(cwavef,ghc_vxctau,gs_hamkq%gbound_k,gs_hamkq%gprimd,gs_hamkq%istwf_k,& 525 & gs_hamkq%kg_k,gs_hamkq%kpt_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw_k,& 526 & gs_hamkq%nvloc,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,nspinor,gs_hamkq%vxctaulocal,& 527 & gs_hamkq%gpu_option) 528 ! There is an additional factor of 2 with respect to the bare matrix element 529 evxctau0_k(iband)=energy_factor*(DOT_PRODUCT(cwavef(1,1:npw_k*nspinor),ghc_vxctau(1,1:npw_k*nspinor))+& 530 & DOT_PRODUCT(cwavef(2,1:npw_k*nspinor),ghc_vxctau(2,1:npw_k*nspinor))) 531 ABI_FREE(ghc_vxctau) 532 end if 533 534 ! Compute the 1-order nuclear dipole contribution (with cwave1 and cwave0), if needed. 535 ! only relevant for DDK 536 if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(rf_hamkq%vectornd)) ) then 537 ABI_MALLOC(ghc_vectornd,(2,npw1_k*nspinor)) 538 ! ndat hard-coded as 1 539 call getgh1ndc(cwave1,ghc_vectornd,gs_hamkq%gbound_k,gs_hamkq%istwf_k,gs_hamkq%kg_k,& 540 & gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw1_k,gs_hamkq%nvloc,& 541 & gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,nspinor,rf_hamkq%vectornd,gs_hamkq%gpu_option) 542 ! There is an additional factor of 4 with respect to the bare matrix element 543 end1_k(iband)=two*energy_factor*(DOT_PRODUCT(cwave0(1,1:npw_k*nspinor),ghc_vectornd(1,1:npw_k*nspinor))+& 544 & DOT_PRODUCT(cwave0(2,1:npw_k*nspinor),ghc_vectornd(2,1:npw_k*nspinor))) 545 ABI_FREE(ghc_vectornd) 546 end if 547 548 ! Compute the 1-order vxctau contribution (with cwave1 and cwave0), if needed. 549 ! only relevant for DDK 550 if( (ipert .EQ. natom+1) .AND. (ASSOCIATED(rf_hamkq%vxctaulocal)) ) then 551 ABI_MALLOC(ghc_vxctau,(2,npw1_k*nspinor)) 552 ! ndat hard-coded as 1 553 call getgh1c_mGGA(cwave1,rf_hamkq%dkinpw_k,gs_hamkq%gbound_k,ghc_vxctau,gs_hamkq%gprimd,idir,gs_hamkq%istwf_k,& 554 & gs_hamkq%kg_k,kinpw1,gs_hamkq%mgfft,mpi_enreg,nspinor,gs_hamkq%n4,gs_hamkq%n5,& 555 & gs_hamkq%n6,1,gs_hamkq%ngfft,npw_k,gs_hamkq%nvloc,rf_hamkq%vxctaulocal,gs_hamkq%gpu_option) 556 ! There is an additional factor of 4 with respect to the bare matrix element 557 evxctau1_k(iband)=two*energy_factor*(DOT_PRODUCT(cwave0(1,1:npw_k*nspinor),ghc_vxctau(1,1:npw_k*nspinor))+& 558 & DOT_PRODUCT(cwave0(2,1:npw_k*nspinor),ghc_vxctau(2,1:npw_k*nspinor))) 559 ABI_FREE(ghc_vxctau) 560 end if 561 ! 562 ! Compute eigenvalue part of total energy (with cwavef) 563 if (gs_hamkq%usepaw==1) then 564 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gsc,mpi_enreg%me_g0,& 565 & mpi_enreg%comm_spinorfft) 566 else 567 call sqnorm_g(scprod,gs_hamkq%istwf_k,npw1_k*nspinor,cwavef,mpi_enreg%me_g0,& 568 & mpi_enreg%comm_fft) 569 end if 570 eeig0_k(iband)=-energy_factor*(eig0_k(iband)- (dtset%dfpt_sciss) )*scprod 571 572 ! Compute nonlocal psp contributions to nonlocal energy: 573 ! <G|Vnl+VFockACE|C1nk(perp)> is contained in gvnlxc (with cwavef) 574 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gvnlxc,mpi_enreg%me_g0,& 575 & mpi_enreg%comm_spinorfft) 576 enl0_k(iband)=energy_factor*scprod 577 578 if(ipert/=natom+10.and.ipert/=natom+11) then 579 ! <G|Vnl1|Cnk> is contained in gvnlx1 (with cwave1) 580 ! gvnlx1 contains at this stage first order kinetic energy, first order nuclear dipole, 581 ! first order vxctau1 582 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,gvnlx1,mpi_enreg%me_g0,& 583 & mpi_enreg%comm_spinorfft) 584 enl1_k(iband)=two*energy_factor*scprod 585 end if 586 587 ! Removal of the 1st-order kinetic energy from the 1st-order non-local part. 588 if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4) then 589 enl1_k(iband)=enl1_k(iband)-ek1_k(iband) 590 end if 591 ! enl1_k still contains first order nuclear dipole, first order vxctau1 in addition to 592 ! first order nonlocal 593 594 ! Accumulate 1st-order density (only at the last inonsc) 595 ! Accumulate zero-order potential part of the 2nd-order total energy 596 ! BUGFIX from Max Stengel: need to initialize eloc at each inonsc iteration, in case nnonsc > 1 597 eloc0_k(iband) = zero 598 option=2;if (iscf_mod>0.and.inonsc==nnsclo_now) option=3 599 call dfpt_accrho(cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,eloc0_k(iband),& 600 & gs_hamkq,iband,idir,ipert,isppol,dtset%kptopt,mpi_enreg,natom,nband_k,ncpgr,& 601 & npw_k,npw1_k,nspinor,occ_k,option,pawrhoij1,rhoaug1,tim_fourwf,tocceig,wtk_k) 602 if(ipert==natom+10.or.ipert==natom+11) eloc0_k(iband)=energy_factor*eloc0_k(iband)/two 603 604 if(ipert==natom+10.or.ipert==natom+11) then 605 shift_band=(iband-1)*npw1_k*nspinor 606 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,& 607 & rf2%RHS_Stern(:,1+shift_band:npw1_k*nspinor+shift_band),mpi_enreg%me_g0, mpi_enreg%comm_spinorfft) 608 ek1_k(iband)=two*energy_factor*scprod 609 end if 610 611 end if ! End of non-zero occupation 612 613 ! Exit loop over inonsc if converged and non-self-consistent 614 if (iscf_mod<0 .and. resid<dtset%tolwfr) exit 615 616 end do ! End loop over inonsc 617 618 ! Get first-order eigenvalues and wavefunctions 619 ptr = 1+(iband_me-1)*npw1_k*nspinor+icg1 620 if (.not. present(cg1_out)) then 621 call cg_zcopy(npw1_k*nspinor,cwave1,cg1(1,ptr)) 622 end if 623 if(dim_eig2rf > 0) then 624 if (.not. present(cg1_out)) then 625 cg1_active(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)=cwavef(:,:) 626 end if 627 gh1c_set(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)=gh1c_n(:,:) 628 gh0c1_set(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband_me*npw1_k*nspinor+icg1)=gh0c1(:,:) 629 end if 630 631 ! PAW: write first-order projected wavefunctions 632 if (psps%usepaw==1.and.gs_hamkq%usecprj==1) then 633 call pawcprj_put(gs_hamkq%atindx,cwaveprj,cprj1,natom,iband_me,ibg1,ikpt,iorder_cprj1,isppol,& 634 & mband_mem,mk1mem,natom,1,nband_me,gs_hamkq%dimcprj,nspinor,nsppol,dtfil%unpaw1) 635 !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,to_be_gathered=.true.) 636 end if 637 638 end do 639 640 !====================================================================== 641 !================== END LOOP OVER BANDS ============================== 642 !====================================================================== 643 644 ! select eig1 matrix for only my slice of bands at present k-point 645 ! this enables global xmpi_sum in dfpt_vtorho without double counting 646 ii = 0 647 do iband=1, nband_k 648 if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) then 649 eig1_k(ii+1:ii+2*nband_k) = zero 650 end if 651 ii = ii + 2*nband_k 652 end do 653 654 ! call xmpi_sum(eig1_k,mpi_enreg%comm_band,ierr) 655 656 ! NB: no need to sum eXX_k over band communicator here, as it is a sub-comm of kpt, 657 ! and full mpi_sum will be done higher up. 658 659 !For rf2 perturbation 660 if(ipert==natom+10.or.ipert==natom+11) call rf2_destroy(rf2) 661 662 !Find largest resid over bands at this k point 663 residk=maxval(resid_k(:)) 664 if (prtvol>2 .or. ikpt<=nkpt_max) then 665 do ii=0,(nband_k-1)/8 666 write(message,'(1p,8e10.2)')(resid_k(iband),iband=1+ii*8,min(nband_k,8+ii*8)) 667 call wrtout(std_out,message) 668 end do 669 end if 670 671 call timab(139,2,tsec) 672 call timab(130,1,tsec) 673 674 ABI_FREE(cwave0) 675 ABI_FREE(cwavef) 676 ABI_FREE(cwave1) 677 ABI_FREE(gh0c1) 678 ABI_FREE(gvnlxc) 679 ABI_FREE(gvnlx1) 680 ABI_FREE(gh1c_n) 681 682 if (gs_hamkq%usepaw==1) then 683 call pawcprj_free(cwaveprj) 684 call pawcprj_free(cwaveprj1) 685 if (gs_hamkq%usecprj==1) then 686 call pawcprj_free(cwaveprj0) 687 end if 688 end if 689 ABI_FREE(dcwavef) 690 ABI_FREE(gscq) 691 ABI_FREE(gsc) 692 ABI_FREE(cwaveprj0) 693 ABI_FREE(cwaveprj) 694 ABI_FREE(cwaveprj1) 695 696 697 !################################################################### 698 699 ! Write the number of one-way 3D ffts skipped until now (in case of fixed occupation numbers) 700 call xmpi_sum(nskip,mpi_enreg%comm_band,ierr) 701 if (iscf_mod>0 .and. (prtvol>2 .or. ikpt<=nkpt_max)) then 702 write(message,'(a,i0)')' dfpt_vtowfk : number of one-way 3D ffts skipped in vtowfk3 until now =',nskip 703 call wrtout(std_out,message) 704 end if 705 706 if (prtvol<=2 .and. ikpt==nkpt_max+1) then 707 write(message,'(3a)') ch10,' dfpt_vtowfk : prtvol=0, 1 or 2, do not print more k-points.',ch10 708 call wrtout(std_out,message) 709 end if 710 711 if (residk>dtset%tolwfr .and. iscf_mod<=0 .and. iscf_mod/=-3) then 712 write(message,'(a,2i0,a,es13.5)')'Wavefunctions not converged for nnsclo,ikpt=',nnsclo_now,ikpt,' max resid=',residk 713 ABI_WARNING(message) 714 end if 715 716 call timab(130,2,tsec) 717 call timab(128,2,tsec) 718 719 DBG_EXIT('COLL') 720 721 end subroutine dfpt_vtowfk
ABINIT/full_active_wf1 [ Functions ]
NAME
full_active_wf1
FUNCTION
Response function calculation only: Restore the full "active space" contribution to the 1st-order wavefunctions. The 1st-order WF corrected in this way will no longer be orthogonal to the other occupied states. This routine will be only used in a non self-consistent calculation of the 1st-order WF for post-processing purposes. Therefore, it does not compute the contribution of the 2DTE coming from the change of occupations.
INPUTS
cgq(2,mcgq)=planewave coefficients of wavefunctions at k+q
cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
cwaveprj(natom,nspinor)= 1st-order wave-function before correction
projected on NL projectors (PAW)
cycle_bands(nband)=array of logicals for bands we have on this cpu
eig1(2*nband**2)=first-order eigenvalues (hartree)
fermie1=derivative of fermi energy wrt (strain) perturbation
eig0nk=energy of the band at k being corrected
eig0_kq(nband)=energies of the bands at k+q
elph2_imagden=imaginary parameter to broaden the energy denominators
iband=index of current band
ibgq=shift to be applied on the location of data in the array cprjq
icgq=shift to be applied on the location of data in the array cgq
mcgq=second dimension of the cgq array
mcprjq=second dimension of the cprjq array
mpi_enreg=information about MPI parallelization
natom=number of atoms in cell
nband=number of bands
npw1=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
usepaw=flag for PAW
OUTPUT
cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
projected on NL projectors (PAW)
SOURCE
769 subroutine full_active_wf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,cycle_bands,eig1,& 770 & fermie1,eig0nk,eig0_kq,elph2_imagden,& 771 & iband,ibgq,icgq,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,& 772 & nspinor,timcount,usepaw) 773 774 !Arguments ------------------------------------ 775 !scalars 776 integer,intent(in) :: iband,ibgq,icgq,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw 777 real(dp),intent(in) :: fermie1, eig0nk 778 real(dp),intent(in) :: elph2_imagden 779 type(MPI_type),intent(in) :: mpi_enreg 780 !arrays 781 logical,intent(in) :: cycle_bands(nband) 782 real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor) 783 real(dp),intent(in) :: eig0_kq(nband) 784 real(dp),intent(in) :: eig1(2*nband**2) 785 real(dp),intent(out) :: cwave1(2,npw1*nspinor) 786 type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw) 787 type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw) !vz_i 788 789 !Local variables------------------------------- 790 !scalars 791 integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii 792 integer :: ibandkq_me, ierr 793 real(dp) :: facti,factr,eta,delta_E,inv_delta_E,gkkr 794 !arrays 795 real(dp) :: tsec(2) 796 797 ! ********************************************************************* 798 799 DBG_ENTER("COLL") 800 ABI_UNUSED(mpi_enreg%comm_cell) 801 802 call timab(214+timcount,1,tsec) 803 804 !At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf). 805 !Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 . 806 807 ! New logic 11/11/2019: accumulate correction in cwave1 and cwaveprj1, 808 ! then add it to cwavef at the end with a modified blas call 809 cwave1 = zero 810 811 if (usepaw==1) then 812 call pawcprj_set_zero(cwaveprj1) 813 end if 814 815 eta = elph2_imagden 816 817 !Loop over WF at k+q subspace 818 ibandkq_me = 0 819 do ibandkq=1,nband 820 821 !TODO MJV: here we have an issue - the cgq are no longer present for all bands! 822 ! we only have diagonal terms for iband iband1 and ibandq in same set of bands 823 ! 1) filter with distrb 824 if(cycle_bands(ibandkq)) cycle 825 ibandkq_me = ibandkq_me + 1 826 827 ! 2) get contributions for correction factors of cgq from bands present on this cpu 828 829 delta_E = eig0nk - eig0_kq(ibandkq) 830 inv_delta_E = delta_E / ( delta_E ** 2 + eta ** 2) 831 832 index_eig1=2*ibandkq-1+(iband-1)*2*nband 833 index_cgq=npw1*nspinor*(ibandkq_me-1)+icgq 834 835 if(ibandkq==iband) then 836 gkkr = eig1(index_eig1) - fermie1 837 else 838 gkkr = eig1(index_eig1) 839 end if 840 factr = inv_delta_E * gkkr 841 facti = inv_delta_E * eig1(index_eig1+1) 842 843 ! Apply correction to 1st-order WF 844 !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwave1,facti,factr,index_cgq,npw1,nspinor) 845 do ii=1,npw1*nspinor 846 cwave1(1,ii)=cwave1(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq)) 847 cwave1(2,ii)=cwave1(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq)) 848 end do 849 850 ! In the PAW case, also apply correction to projected WF 851 if (usepaw==1) then 852 index_cprjq=nspinor*(ibandkq_me-1)+ibgq 853 call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1) 854 end if 855 856 end do ! Loop over k+q subspace 857 858 ! 3) reduce over bands to get all contributions to correction 859 ! need MPI reduce over band communicator only 860 call xmpi_sum(cwave1,mpi_enreg%comm_band,ierr) 861 if (usepaw==1) then 862 call pawcprj_mpi_sum(cwaveprj1, mpi_enreg%comm_band, ierr) 863 end if 864 865 ! 4) add correction to the cwave1 866 !Now add on input WF into output WF 867 call cg_zaxpy(npw1*nspinor,(/one,zero/),cwavef,cwave1) 868 869 !Idem for cprj 870 if (usepaw==1) then 871 call pawcprj_zaxpby((/one,zero/),(/one,zero/),cwaveprj,cwaveprj1) 872 end if 873 874 call timab(214+timcount,2,tsec) 875 876 DBG_EXIT("COLL") 877 878 end subroutine full_active_wf1
ABINIT/m_dfpt_vtowfk [ Modules ]
NAME
m_dfpt_vtowfk
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XG, AR, DRH, MB, MVer,XW, MT, GKA) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_dfpt_vtowfk 22 23 use defs_basis 24 use m_abicore 25 use m_errors 26 use m_xmpi 27 use m_mpinfo 28 use m_cgtools 29 use m_wfk 30 use m_rf2 31 use m_dtset 32 use m_dtfil 33 34 use defs_datatypes, only : pseudopotential_type 35 use defs_abitypes, only : MPI_type 36 use m_rf2_init, only : rf2_init 37 use m_time, only : timab 38 use m_pawrhoij, only : pawrhoij_type 39 use m_pawcprj 40 use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type, KPRIME_H_KPRIME 41 use m_spacepar, only : meanvalue_g 42 use m_dfpt_mkrho, only : dfpt_accrho 43 use m_dfpt_cgwf, only : dfpt_cgwf 44 use m_getghc, only : getgsc, getghc_nucdip, getghc_mGGA 45 use m_getgh1c, only : getgh1ndc, getgh1c_mGGA 46 47 implicit none 48 49 private