TABLE OF CONTENTS
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).
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, AR, DRH, MB, MVer,XW, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
cg1(2,mpw1*nspinor*mband*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=number of the k-point
ipert=type of the perturbation
isppol=1 index of current spin component
mband=maximum number of bands
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*mk1mem*nsppol)=pw coefficients of RF
wavefunctions at k,q. They are orthogonalized to the occupied states.
cg1_active(2,mpw1*nspinor*mband*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.
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.
gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
gh0c1_set(2,mpw1*nspinor*mband*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)
PARENTS
dfpt_vtorho
CHILDREN
cg_zcopy,corrmetalwf1,dfpt_accrho,dfpt_cgwf,dotprod_g,getgsc matrixelmt_g,meanvalue_g,pawcprj_alloc,pawcprj_copy,pawcprj_free pawcprj_get,pawcprj_put,rf2_destroy,rf2_init,sqnorm_g,status,timab wfk_read_bks,wrtout
SOURCE
118 #if defined HAVE_CONFIG_H 119 #include "config.h" 120 #endif 121 122 #include "abi_common.h" 123 124 subroutine dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,& 125 & dim_eig2rf,dtfil,dtset,& 126 & edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,& 127 & ek0_k,ek1_k,eloc0_k,enl0_k,enl1_k,& 128 & fermie1,gh0c1_set,gh1c_set,grad_berry,gs_hamkq,& 129 & ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,& 130 & isppol,mband,mcgq,mcprjq,mkmem,mk1mem,& 131 & mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,& 132 & nnsclo_now,npw_k,npw1_k,nspinor,nsppol,& 133 & n4,n5,n6,occ_k,pawrhoij1,prtvol,psps,resid_k,rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,& 134 & ddk_f,wtk_k,nlines_done,cg1_out) 135 136 use defs_basis 137 use defs_datatypes 138 use defs_abitypes 139 use m_profiling_abi 140 use m_errors 141 use m_xmpi 142 use m_cgtools 143 use m_wfk 144 use m_rf2 145 146 use m_pawrhoij, only : pawrhoij_type 147 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_free, pawcprj_get,pawcprj_copy 148 use m_hamiltonian, only : gs_hamiltonian_type,rf_hamiltonian_type,KPRIME_H_KPRIME 149 150 !This section has been created automatically by the script Abilint (TD). 151 !Do not modify the following lines by hand. 152 #undef ABI_FUNC 153 #define ABI_FUNC 'dfpt_vtowfk' 154 use interfaces_14_hidewrite 155 use interfaces_18_timing 156 use interfaces_32_util 157 use interfaces_53_spacepar 158 use interfaces_66_wfs 159 use interfaces_72_response, except_this_one => dfpt_vtowfk 160 !End of the abilint section 161 162 implicit none 163 164 !Arguments ------------------------------------ 165 !scalars 166 integer,intent(in) :: cplex,dim_eig2rf,ibg 167 integer,intent(in) :: ibg1,ibgq,icg,icg1,icgq,idir,ikpt,ipert,isppol 168 integer,intent(in) :: mband,mcgq,mcprjq,mk1mem,mkmem 169 integer,intent(in) :: mpw,mpw1,n4,n5,n6,natom,ncpgr 170 integer,intent(in) :: nnsclo_now,nspinor,nsppol,prtvol 171 integer,optional,intent(in) :: cg1_out 172 integer,intent(in) :: nband_k,npw1_k,npw_k 173 integer,intent(inout) :: nlines_done 174 real(dp),intent(in) :: fermie1,wtk_k 175 type(MPI_type),intent(in) :: mpi_enreg 176 type(datafiles_type),intent(in) :: dtfil 177 type(dataset_type),intent(in) :: dtset 178 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 179 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq,rf_hamk_dir2 180 type(pseudopotential_type),intent(in) :: psps 181 !arrays 182 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),cgq(2,mcgq) 183 real(dp),intent(in) :: eig0_k(nband_k),eig0_kq(nband_k) 184 real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband_k) 185 real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k) 186 real(dp),intent(inout) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) 187 real(dp),intent(inout) :: rhoaug1(cplex*n4,n5,n6,gs_hamkq%nvloc) 188 real(dp),intent(inout) :: cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf) 189 real(dp),intent(inout) :: gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf) 190 real(dp),intent(inout) :: gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf) 191 real(dp),intent(inout) :: edocc_k(nband_k),eeig0_k(nband_k),eig1_k(2*nband_k**2) 192 real(dp),intent(out) :: ek0_k(nband_k),eloc0_k(nband_k) 193 real(dp),intent(inout) :: ek1_k(nband_k) 194 real(dp),intent(out) :: enl0_k(nband_k),enl1_k(nband_k) 195 real(dp),intent(out) :: resid_k(nband_k) 196 type(pawcprj_type),intent(in) :: cprj(natom,nspinor*mband*mkmem*nsppol*gs_hamkq%usecprj) 197 type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq) 198 type(pawcprj_type),intent(inout) :: cprj1(natom,nspinor*mband*mk1mem*nsppol*gs_hamkq%usecprj) 199 type(pawrhoij_type),intent(inout) :: pawrhoij1(natom*gs_hamkq%usepaw) 200 type(wfk_t),intent(inout) :: ddk_f(4) 201 202 !Local variables------------------------------- 203 !scalars 204 integer,parameter :: level=14,tim_fourwf=5 205 integer,save :: nskip=0 206 integer :: counter,iband,idir0,ierr,iexit,igs,igscq,ii,dim_dcwf,inonsc 207 integer :: iorder_cprj,iorder_cprj1,ipw,iscf_mod,ispinor,me,mgscq,nkpt_max 208 integer :: option,opt_gvnl1,quit,test_ddk 209 integer :: tocceig,usedcwavef,ptr,shift_band 210 real(dp) :: aa,ai,ar,eig0nk,resid,residk,scprod,energy_factor 211 character(len=500) :: message 212 type(rf2_t) :: rf2 213 !arrays 214 real(dp) :: tsec(2) 215 real(dp),allocatable :: cwave0(:,:),cwave1(:,:),cwavef(:,:) 216 real(dp),allocatable :: dcwavef(:,:),gh1c_n(:,:),gh0c1(:,:) 217 real(dp),allocatable :: gsc(:,:),gscq(:,:),gvnl1(:,:),gvnlc(:,:) 218 real(dp),pointer :: kinpw1(:) 219 type(pawcprj_type),allocatable :: cwaveprj(:,:),cwaveprj0(:,:),cwaveprj1(:,:) 220 221 ! ********************************************************************* 222 223 DBG_ENTER('COLL') 224 225 !Keep track of total time spent in dfpt_vtowfk 226 call timab(128,1,tsec) 227 228 nkpt_max=50; if (xmpi_paral==1) nkpt_max=-1 229 230 if(prtvol>2 .or. ikpt<=nkpt_max)then 231 write(message,'(2a,i5,2x,a,3f9.5,2x,a)')ch10,' Non-SCF iterations; k pt #',ikpt,'k=',& 232 & gs_hamkq%kpt_k(:),'band residuals:' 233 call wrtout(std_out,message,'PERS') 234 end if 235 236 !Initializations and allocations 237 me=mpi_enreg%me_kpt 238 quit=0 239 240 !The value of iscf must be modified if ddk perturbation 241 iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3 242 243 kinpw1 => gs_hamkq%kinpw_kp 244 ABI_ALLOCATE(gh0c1,(2,npw1_k*nspinor)) 245 ABI_ALLOCATE(gvnlc,(2,npw1_k*nspinor)) 246 ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor)) 247 ABI_ALLOCATE(cwave0,(2,npw_k*nspinor)) 248 ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor)) 249 ABI_ALLOCATE(cwave1,(2,npw1_k*nspinor)) 250 ABI_ALLOCATE(gh1c_n,(2,npw1_k*nspinor)) 251 if (gs_hamkq%usepaw==1) then 252 ABI_ALLOCATE(gsc,(2,npw1_k*nspinor)) 253 else 254 ABI_ALLOCATE(gsc,(0,0)) 255 end if 256 257 !Read the npw and kg records of wf files 258 call status(0,dtfil%filstat,iexit,level,'before WffRead') 259 test_ddk=0 260 if ((ipert==natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and.& 261 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 262 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.& 263 & ipert==natom+10.or.ipert==natom+11) then 264 test_ddk=1 265 if(ipert==natom+10.or.ipert==natom+11) test_ddk=0 266 end if 267 268 !Additional stuff for PAW 269 ABI_DATATYPE_ALLOCATE(cwaveprj0,(0,0)) 270 if (gs_hamkq%usepaw==1) then 271 ! 1-Compute all <g|S|Cnk+q> 272 igscq=0 273 mgscq=mpw1*nspinor*mband 274 ABI_STAT_ALLOCATE(gscq,(2,mgscq), ierr) 275 ABI_CHECK(ierr==0, "out of memory in gscq") 276 277 call getgsc(cgq,cprjq,gs_hamkq,gscq,ibgq,icgq,igscq,ikpt,isppol,mcgq,mcprjq,& 278 & mgscq,mpi_enreg,natom,nband_k,npw1_k,dtset%nspinor,select_k=KPRIME_H_KPRIME) 279 ! 2-Initialize additional scalars/arrays 280 iorder_cprj=0;iorder_cprj1=0 281 dim_dcwf=npw1_k*nspinor;if (ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) dim_dcwf=0 282 ABI_ALLOCATE(dcwavef,(2,dim_dcwf)) 283 if (gs_hamkq%usecprj==1) then 284 ABI_DATATYPE_DEALLOCATE(cwaveprj0) 285 ABI_DATATYPE_ALLOCATE(cwaveprj0,(natom,nspinor)) 286 call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj) 287 end if 288 ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,nspinor)) 289 ABI_DATATYPE_ALLOCATE(cwaveprj1,(natom,nspinor)) 290 call pawcprj_alloc(cwaveprj ,0,gs_hamkq%dimcprj) 291 call pawcprj_alloc(cwaveprj1,0,gs_hamkq%dimcprj) 292 else 293 igscq=0;mgscq=0;dim_dcwf=0 294 ABI_ALLOCATE(gscq,(0,0)) 295 ABI_ALLOCATE(dcwavef,(0,0)) 296 ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0)) 297 ABI_DATATYPE_ALLOCATE(cwaveprj1,(0,0)) 298 end if 299 300 energy_factor=two 301 if(ipert==natom+10.or.ipert==natom+11) energy_factor=six 302 303 !For rf2 perturbation : 304 if(ipert==natom+10.or.ipert==natom+11) then 305 call rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,& 306 mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f) 307 end if 308 309 call timab(139,1,tsec) 310 311 !====================================================================== 312 !================== LOOP OVER BANDS ================================== 313 !====================================================================== 314 315 do iband=1,nband_k 316 317 ! Skip bands not treated by current proc 318 if( (mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me)) cycle 319 320 ! Get ground-state wavefunctions 321 ptr = 1+(iband-1)*npw_k*nspinor+icg 322 call cg_zcopy(npw_k*nspinor,cg(1,ptr),cwave0) 323 324 ! Get PAW ground state projected WF (cprj) 325 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.ipert/=natom+10.and.ipert/=natom+11) then 326 idir0 = idir 327 if(ipert==natom+3.or.ipert==natom+4) idir0 =1 328 call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,natom,iband,ibg,ikpt,iorder_cprj,& 329 & isppol,mband,mkmem,natom,1,nband_k,nspinor,nsppol,dtfil%unpaw,& 330 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,& 331 & icpgr=idir0,ncpgr=ncpgr) 332 end if 333 334 ! Get first-order wavefunctions 335 ptr = 1+(iband-1)*npw1_k*nspinor+icg1 336 call cg_zcopy(npw1_k*nspinor,cg1(1,ptr),cwavef) 337 338 ! Read PAW projected 1st-order WF (cprj) 339 ! Unuseful for the time being (will be recomputed in dfpt_cgwf) 340 ! if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 341 ! call pawcprj_get(gs_hamkq%atindx1,cwaveprj,cprj1,natom,iband,ibg1,ikpt,iorder_cprj1,& 342 ! & isppol,mband,mk1mem,natom,1,nband_k,nspinor,nsppol,dtfil%unpaw1, 343 ! & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 344 ! end if 345 346 ! Filter the wavefunctions for large modified kinetic energy 347 ! The GS wavefunctions should already be non-zero 348 do ispinor=1,nspinor 349 igs=(ispinor-1)*npw1_k 350 do ipw=1+igs,npw1_k+igs 351 if(kinpw1(ipw-igs)>huge(zero)*1.d-11)then 352 cwavef(1,ipw)=zero 353 cwavef(2,ipw)=zero 354 end if 355 end do 356 end do 357 358 359 ! If electric field, the derivative of the wf should be read, and multiplied by i. 360 if(test_ddk==1) then 361 ii = wfk_findk(ddk_f(1), gs_hamkq%kpt_k) 362 ABI_CHECK(ii == ikpt, "ii != ikpt") 363 call wfk_read_bks(ddk_f(1), iband, ikpt, isppol, xmpio_single, cg_bks=gvnl1) 364 365 ! Multiplication by -i 366 ! MVeithen 021212 : use + i instead, 367 ! See X. Gonze, Phys. Rev. B 55, 10337 (1997) Eq. (79) 368 ! the operator used to compute the first-order derivative 369 ! of the wavefunctions with respect to an electric field 370 ! is $+i \frac{d}{dk}$ 371 ! This change will affect the computation of the 2dtes from non 372 ! stationary expressions, see dfpt_nstdy.f and dfpt_nstwf.f 373 do ipw=1,npw1_k*nspinor 374 ! aa=gvnl1(1,ipw) 375 ! gvnl1(1,ipw)=gvnl1(2,ipw) 376 ! gvnl1(2,ipw)=-aa 377 aa=gvnl1(1,ipw) 378 gvnl1(1,ipw)=-gvnl1(2,ipw) 379 gvnl1(2,ipw)=aa 380 end do 381 end if 382 383 ! Unlike in GS calculations, the inonsc loop is inside the band loop 384 ! nnsclo_now=number of non-self-consistent loops for the current vtrial 385 ! (often 1 for SCF calculation, =nstep for non-SCF calculations) 386 do inonsc=1,nnsclo_now 387 388 counter=100*iband+inonsc 389 390 ! Note that the following translation occurs in the called routine : 391 ! iband->band, nband_k->nband, npw_k->npw, npw1_k->npw1 392 eig0nk=eig0_k(iband) 393 usedcwavef=gs_hamkq%usepaw;if (dim_dcwf==0) usedcwavef=0 394 if (inonsc==1) usedcwavef=2*usedcwavef 395 opt_gvnl1=0;if (ipert==natom+2) opt_gvnl1=1 396 if (ipert==natom+2.and.gs_hamkq%usepaw==1.and.inonsc==1) opt_gvnl1=2 397 398 if ( (ipert/=natom+10 .and. ipert/=natom+11) .or. abs(occ_k(iband))>tol8 ) then 399 call dfpt_cgwf(iband,dtset%berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,rf2,dcwavef,& 400 & eig0nk,eig0_kq,eig1_k,gh0c1,gh1c_n,grad_berry,gsc,gscq,gs_hamkq,gvnlc,gvnl1,icgq,& 401 & idir,ipert,igscq,mcgq,mgscq,mpi_enreg,mpw1,natom,nband_k,dtset%nbdbuf,dtset%nline,& 402 & npw_k,npw1_k,nspinor,opt_gvnl1,prtvol,quit,resid,rf_hamkq,dtset%dfpt_sciss,dtset%tolrde,& 403 & dtset%tolwfr,usedcwavef,dtset%wfoptalg,nlines_done) 404 resid_k(iband)=resid 405 else 406 resid_k(iband)=zero 407 end if 408 409 if (ipert/=natom+10 .and. ipert/= natom+11) then 410 ! At this stage, the 1st order function cwavef is orthogonal to cgq (unlike 411 ! when it is input to dfpt_cgwf). Here, restore the "active space" content 412 ! of the first-order wavefunction, to give cwave1. 413 ! PAW: note that dcwavef (1st-order change of WF due to overlap change) 414 ! remains in the subspace orthogonal to cgq 415 call corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,edocc_k,eig1_k,fermie1,gh0c1,& 416 & iband,ibgq,icgq,gs_hamkq%istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,& 417 & occ_k,rocceig,0,gs_hamkq%usepaw,tocceig) 418 else 419 tocceig=0 420 call cg_zcopy(npw1_k*nspinor,cwavef,cwave1) 421 if (gs_hamkq%usepaw==1) then 422 call pawcprj_copy(cwaveprj,cwaveprj1) 423 end if 424 end if 425 426 if (abs(occ_k(iband))<= tol8) then 427 ek0_k(iband)=zero 428 ek1_k(iband)=zero 429 eeig0_k(iband)=zero 430 enl0_k(iband)=zero 431 enl1_k(iband)=zero 432 eloc0_k(iband)=zero 433 nskip=nskip+1 434 else 435 ! Compute the 0-order kinetic operator contribution (with cwavef) 436 call meanvalue_g(ar,kinpw1,0,gs_hamkq%istwf_k,mpi_enreg,npw1_k,nspinor,cwavef,cwavef,0) 437 ! There is an additional factor of 2 with respect to the bare matrix element 438 ek0_k(iband)=energy_factor*ar 439 ! Compute the 1-order kinetic operator contribution (with cwave1 and cwave0), if needed. 440 ! Note that this is called only for ddk or strain, so that npw1_k=npw_k 441 if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4)then 442 call matrixelmt_g(ai,ar,rf_hamkq%dkinpw_k,gs_hamkq%istwf_k,0,npw_k,nspinor,cwave1,cwave0,& 443 & mpi_enreg%me_g0, mpi_enreg%comm_fft) 444 ! There is an additional factor of 4 with respect to the bare matrix element 445 ek1_k(iband)=two*energy_factor*ar 446 end if 447 448 ! Compute eigenvalue part of total energy (with cwavef) 449 if (gs_hamkq%usepaw==1) then 450 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gsc,mpi_enreg%me_g0,& 451 & mpi_enreg%comm_spinorfft) 452 else 453 call sqnorm_g(scprod,gs_hamkq%istwf_k,npw1_k*nspinor,cwavef,mpi_enreg%me_g0,& 454 & mpi_enreg%comm_fft) 455 end if 456 eeig0_k(iband)=-energy_factor*(eig0_k(iband)- (dtset%dfpt_sciss) )*scprod 457 458 ! Compute nonlocal psp contributions to nonlocal energy: 459 ! <G|Vnl|C1nk(perp)> is contained in gvnlc (with cwavef) 460 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gvnlc,mpi_enreg%me_g0,& 461 & mpi_enreg%comm_spinorfft) 462 enl0_k(iband)=energy_factor*scprod 463 464 if(ipert/=natom+10.and.ipert/=natom+11) then 465 ! <G|Vnl1|Cnk> is contained in gvnl1 (with cwave1) 466 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,gvnl1,mpi_enreg%me_g0,& 467 & mpi_enreg%comm_spinorfft) 468 enl1_k(iband)=two*energy_factor*scprod 469 end if 470 471 ! Removal of the 1st-order kinetic energy from the 1st-order non-local part. 472 if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4) then 473 enl1_k(iband)=enl1_k(iband)-ek1_k(iband) 474 end if 475 476 ! Accumulate 1st-order density (only at the last inonsc) 477 ! Accumulate zero-order potential part of the 2nd-order total energy 478 ! BUGFIX from Max Stengel: need to initialize eloc at each inonsc iteration, in case nnonsc > 1 479 eloc0_k(iband) = zero 480 option=2;if (iscf_mod>0.and.inonsc==nnsclo_now) option=3 481 call dfpt_accrho(counter,cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,eloc0_k(iband),& 482 & dtfil%filstat,gs_hamkq,iband,idir,ipert,isppol,dtset%kptopt,mpi_enreg,natom,nband_k,ncpgr,& 483 & npw_k,npw1_k,nspinor,occ_k,option,pawrhoij1,prtvol,rhoaug1,tim_fourwf,tocceig,wtk_k) 484 if(ipert==natom+10.or.ipert==natom+11) eloc0_k(iband)=energy_factor*eloc0_k(iband)/two 485 486 if(ipert==natom+10.or.ipert==natom+11) then 487 shift_band=(iband-1)*npw1_k*nspinor 488 call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,& 489 & rf2%RHS_Stern(:,1+shift_band:npw1_k*nspinor+shift_band),mpi_enreg%me_g0, mpi_enreg%comm_spinorfft) 490 ek1_k(iband)=two*energy_factor*scprod 491 end if 492 493 end if ! End of non-zero occupation 494 495 ! Exit loop over inonsc if converged and if non-self-consistent 496 if (iscf_mod<0 .and. resid<dtset%tolwfr) exit 497 498 end do ! End loop over inonsc 499 500 ! Get first-order eigenvalues and wavefunctions 501 ptr = 1+(iband-1)*npw1_k*nspinor+icg1 502 if (.not. present(cg1_out)) then 503 call cg_zcopy(npw1_k*nspinor,cwave1,cg1(1,ptr)) 504 end if 505 if(dim_eig2rf > 0) then 506 if (.not. present(cg1_out)) then 507 cg1_active(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)=cwavef(:,:) 508 end if 509 gh1c_set(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)=gh1c_n(:,:) 510 gh0c1_set(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)=gh0c1(:,:) 511 end if 512 513 ! PAW: write first-order projected wavefunctions 514 if (psps%usepaw==1.and.gs_hamkq%usecprj==1) then 515 call pawcprj_put(gs_hamkq%atindx,cwaveprj,cprj1,natom,iband,ibg1,ikpt,iorder_cprj1,isppol,& 516 & mband,mk1mem,natom,1,nband_k,gs_hamkq%dimcprj,nspinor,nsppol,dtfil%unpaw1,& 517 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,to_be_gathered=.true.) 518 end if 519 520 end do 521 522 !====================================================================== 523 !================== END LOOP OVER BANDS ============================== 524 !====================================================================== 525 526 !For rf2 perturbation 527 if(ipert==natom+10.or.ipert==natom+11) call rf2_destroy(rf2) 528 529 !Find largest resid over bands at this k point 530 residk=maxval(resid_k(:)) 531 if (prtvol>2 .or. ikpt<=nkpt_max) then 532 do ii=0,(nband_k-1)/8 533 write(message,'(1p,8e10.2)')(resid_k(iband),iband=1+ii*8,min(nband_k,8+ii*8)) 534 call wrtout(std_out,message,'PERS') 535 end do 536 end if 537 538 call timab(139,2,tsec) 539 call timab(130,1,tsec) 540 541 ABI_DEALLOCATE(cwave0) 542 ABI_DEALLOCATE(cwavef) 543 ABI_DEALLOCATE(cwave1) 544 ABI_DEALLOCATE(gh0c1) 545 ABI_DEALLOCATE(gvnlc) 546 ABI_DEALLOCATE(gvnl1) 547 ABI_DEALLOCATE(gh1c_n) 548 549 if (gs_hamkq%usepaw==1) then 550 call pawcprj_free(cwaveprj) 551 call pawcprj_free(cwaveprj1) 552 if (gs_hamkq%usecprj==1) then 553 call pawcprj_free(cwaveprj0) 554 end if 555 end if 556 ABI_DEALLOCATE(dcwavef) 557 ABI_DEALLOCATE(gscq) 558 ABI_DEALLOCATE(gsc) 559 ABI_DATATYPE_DEALLOCATE(cwaveprj0) 560 ABI_DATATYPE_DEALLOCATE(cwaveprj) 561 ABI_DATATYPE_DEALLOCATE(cwaveprj1) 562 563 !################################################################### 564 565 !Write the number of one-way 3D ffts skipped until now (in case of fixed occupation numbers) 566 if(iscf_mod>0 .and. (prtvol>2 .or. ikpt<=nkpt_max))then 567 write(message,'(a,i0)')' dfpt_vtowfk : number of one-way 3D ffts skipped in vtowfk3 until now =',nskip 568 call wrtout(std_out,message,'PERS') 569 end if 570 571 if(prtvol<=2 .and. ikpt==nkpt_max+1)then 572 write(message,'(3a)') ch10,' dfpt_vtowfk : prtvol=0, 1 or 2, do not print more k-points.',ch10 573 call wrtout(std_out,message,'PERS') 574 end if 575 576 if (residk>dtset%tolwfr .and. iscf_mod<=0 .and. iscf_mod/=-3) then 577 write(message,'(a,2i0,a,es13.5)')'Wavefunctions not converged for nnsclo,ikpt=',nnsclo_now,& 578 & ikpt,' max resid=',residk 579 MSG_WARNING(message) 580 end if 581 582 call timab(130,2,tsec) 583 call timab(128,2,tsec) 584 585 DBG_EXIT('COLL') 586 587 end subroutine dfpt_vtowfk