TABLE OF CONTENTS
ABINIT/fock_getghc [ Functions ]
NAME
fock_getghc
FUNCTION
Compute the matrix elements <G|Vx|psi> of the Fock operator.
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (CMartins,FJ,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 .
INPUTS
cwavef(2,npw*nspinor*ndat)= planewave coefficients of wavefunctions on which Fock operator is applied. cwaveprj <type(pawcprj_type> = <cwavevf|proj> gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied mpi_enreg= information about MPI parallelization
SIDE EFFECTS
ghc(2,npw*ndat)= matrix elements <G|H|C> or <G|H-lambda.S|C> (if sij_opt>=0 or =-1 in getghc) contains the fock exchange term for cwavef at the end.
NOTES
The current version assumes that : * nspinor = 1 * no "my_nspinor" * no restriction to the value of istwfk_bz (but must be tested in all case) * all the data for the occupied states (cgocc_bz) are the same as those for the current states (cg)
PARENTS
fock2ACE,forstrnps,getghc
CHILDREN
bare_vqg,dotprod_g,fftpac,fourdp,fourwf,hartre,load_kprime_hamiltonian matr3inv,nonlop,pawdijhat,pawmknhat_psipsi,sphereboundary,strfock,timab xmpi_sum
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 subroutine fock_getghc(cwavef,cwaveprj,ghc,gs_ham,mpi_enreg) 49 50 use defs_basis 51 use defs_abitypes 52 use m_errors 53 use m_xmpi 54 use m_cgtools, only :dotprod_g 55 use m_fock 56 use m_hamiltonian, only : gs_hamiltonian_type,load_kprime_hamiltonian,K_H_KPRIME,load_k_hamiltonian 57 58 use m_pawcprj 59 use m_pawdij, only : pawdijhat 60 61 use defs_datatypes, only: pseudopotential_type 62 use m_pawrhoij, only : pawrhoij_type, pawrhoij_free, pawrhoij_alloc 63 64 !This section has been created automatically by the script Abilint (TD). 65 !Do not modify the following lines by hand. 66 #undef ABI_FUNC 67 #define ABI_FUNC 'fock_getghc' 68 use interfaces_18_timing 69 use interfaces_32_util 70 use interfaces_52_fft_mpi_noabirule 71 use interfaces_53_ffts 72 use interfaces_56_xc 73 use interfaces_65_paw 74 use interfaces_66_nonlocal 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 ! Scalars 81 type(MPI_type),intent(in) :: mpi_enreg 82 type(gs_hamiltonian_type),target,intent(inout) :: gs_ham 83 ! Arrays 84 type(pawcprj_type),intent(inout) :: cwaveprj(:,:) 85 real(dp),intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k) 86 real(dp),intent(inout) :: ghc(:,:) 87 88 !Local variables------------------------------- 89 ! Scalars 90 integer,parameter :: tim_fourwf0=0,tim_fourdp0=0,ndat1=1 91 integer :: bdtot_jindex,choice,cplex_fock,cplex_dij,cpopt,i1,i2,i3,ia,iatom 92 integer :: iband_cprj,ider,idir,idir1,ier,ind,ipert,ipw,ifft,itypat,izero,jband,jbg,jcg,jkg 93 integer :: jkpt,my_jsppol,jstwfk,lmn2_size,mgfftf,mpw,n1,n2,n3,n4,n5,n6 94 integer :: n1f,n2f,n3f,n4f,n5f,n6f,natom,nband_k,ndij,nfft,nfftf,nfftotf,nhat12_grdim,nnlout 95 integer :: npw,npwj,nspden_fock,nspinor,paw_opt,signs,tim_nonlop 96 logical :: need_ghc 97 real(dp),parameter :: weight1=one 98 real(dp) :: doti,eigen,imcwf,imcwocc,imvloc,invucvol,recwf,recwocc,revloc,occ,wtk 99 type(fock_common_type),pointer :: fockcommon 100 type(fock_BZ_type),pointer :: fockbz 101 ! Arrays 102 integer :: ngfft(18),ngfftf(18) 103 integer,pointer :: gboundf(:,:),kg_occ(:,:),gbound_kp(:,:) 104 real(dp) :: enlout_dum(1),dotr(6),fockstr(6),for1(3),qphon(3),qvec_j(3),tsec(2),gsc_dum(2,0),rhodum(2,1) 105 real(dp) :: rhodum0(0,1,1),str(3,3) 106 real(dp), allocatable :: dummytab(:,:),dijhat(:,:,:),dijhat_tmp(:,:),ffnl_kp_dum(:,:,:,:) 107 real(dp), allocatable :: gvnlc(:,:),ghc1(:,:),ghc2(:,:),grnhat12(:,:,:,:),grnhat_12(:,:,:,:,:),forikpt(:,:) 108 real(dp), allocatable :: rho12(:,:,:),rhog_munu(:,:),rhor_munu(:,:),vlocpsi_r(:) 109 real(dp), allocatable :: vfock(:),psilocal(:,:,:),vectin_dum(:,:),vqg(:),forout(:,:),strout(:,:) 110 real(dp), allocatable,target ::cwavef_r(:,:,:,:) 111 real(dp), ABI_CONTIGUOUS pointer :: cwaveocc_r(:,:,:,:) 112 type(pawcprj_type),pointer :: cwaveocc_prj(:,:) 113 114 real(dp) :: rprimd(3,3),for12(3) 115 116 ! ************************************************************************* 117 !return 118 call timab(1504,1,tsec) 119 call timab(1505,1,tsec) 120 121 ABI_CHECK(associated(gs_ham%fockcommon),"fock_common must be associated!") 122 fockcommon => gs_ham%fockcommon 123 ABI_CHECK(associated(gs_ham%fockbz),"fock_bz must be associated!") 124 fockbz => gs_ham%fockbz 125 126 ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!") 127 ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!") 128 if (fockcommon%usepaw==1) then 129 ABI_CHECK((size(cwaveprj,1)==gs_ham%natom.and.size(cwaveprj,2)==gs_ham%nspinor),"error on cwaveprj dims") 130 end if 131 need_ghc=(size(ghc,2)>0) 132 133 !Some constants 134 invucvol=1.d0/sqrt(gs_ham%ucvol) 135 call matr3inv(gs_ham%gprimd,rprimd) 136 cplex_fock=2;nspden_fock=1 137 natom=fockcommon%natom 138 nspinor=gs_ham%nspinor 139 mpw=maxval(fockbz%npwarr) 140 npw=gs_ham%npw_k 141 ider=0;izero=0 142 if (fockcommon%usepaw==1) then 143 nfft =fockcommon%pawfgr%nfftc ; ngfft =fockcommon%pawfgr%ngfftc 144 nfftf=fockcommon%pawfgr%nfft ; ngfftf=fockcommon%pawfgr%ngfft 145 mgfftf=fockcommon%pawfgr%mgfft 146 else 147 nfft =gs_ham%nfft ; nfftf =nfft 148 ngfft=gs_ham%ngfft ; ngfftf=ngfft 149 mgfftf=gs_ham%mgfft 150 end if 151 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3) 152 n4=ngfft(4);n5=ngfft(5);n6=ngfft(6) 153 n1f=ngfftf(1);n2f=ngfftf(2);n3f=ngfftf(3) 154 n4f=ngfftf(4);n5f=ngfftf(5);n6f=ngfftf(6) 155 156 ! =========================== 157 ! === Initialize arrays === 158 ! =========================== 159 ! transient optfor and optstress 160 ! fockcommon%optfor=.false. 161 ! fockcommon%optstr=.false. 162 !*Initialization of local pointers 163 !*Initialization of the array cwavef_r 164 !*cwavef_r = current wavefunction in r-space 165 ABI_ALLOCATE(cwavef_r,(2,n4f,n5f,n6f)) 166 !*rhormunu = overlap matrix between cwavef and (jkpt,mu) in R-space 167 ABI_ALLOCATE(rhor_munu,(cplex_fock,nfftf)) 168 !*rhogmunu = overlap matrix between cwavef and (jkpt,mu) in G-space 169 ABI_ALLOCATE(rhog_munu,(2,nfftf)) 170 !*dummytab = variables for fourwf 171 ABI_ALLOCATE(dummytab,(2,nfft)) 172 !*vfock = Fock potential 173 ABI_ALLOCATE(vfock,(cplex_fock*nfftf)) 174 !*vqg = 4pi/(G+q)**2 175 ABI_ALLOCATE(vqg,(nfftf)) 176 177 !*Initialization of the array ghc1 178 !*ghc1 will contain the exact exchange contribution to the Hamiltonian 179 ABI_ALLOCATE(ghc1,(2,npw)) 180 ABI_ALLOCATE(ghc2,(2,npw)) 181 ghc1=zero 182 ghc2=zero 183 !*Initialization of the array vlocpsi_r 184 !*vlocpsi_r = partial local Fock operator applied to cwavef in r-space and summed over all occupied (jkpt,mu) 185 ABI_ALLOCATE(vlocpsi_r,(cplex_fock*nfftf)) 186 vlocpsi_r=zero 187 188 !*Additional arrays in case of paw 189 if (fockcommon%usepaw==1) then 190 nhat12_grdim=0 191 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 192 ider=3 193 ABI_ALLOCATE(strout,(2,npw*nspinor)) 194 end if 195 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 196 ider=3 197 ABI_ALLOCATE(forout,(2,npw*nspinor)) 198 ABI_ALLOCATE(forikpt,(3,natom)) 199 forikpt=zero 200 end if 201 ABI_ALLOCATE(grnhat_12,(2,nfftf,nspinor**2,3,natom*(ider/3))) 202 ABI_ALLOCATE(gvnlc,(2,npw*nspinor)) 203 ABI_ALLOCATE(grnhat12,(2,nfftf,nspinor**2,3*nhat12_grdim)) 204 end if 205 206 if (fockcommon%usepaw==1.or.fockcommon%optstr) then 207 ABI_ALLOCATE(gboundf,(2*mgfftf+8,2)) 208 call sphereboundary(gboundf,gs_ham%istwf_k,gs_ham%kg_k,mgfftf,npw) 209 else 210 gboundf=>gs_ham%gbound_k 211 end if 212 213 ! ========================================== 214 ! === Get cwavef in real space using FFT === 215 ! ========================================== 216 cwavef_r=zero 217 call fourwf(0,rhodum0,cwavef,rhodum,cwavef_r,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,& 218 & mgfftf,mpi_enreg,ndat1,ngfftf,npw,1,n4f,n5f,n6f,0,mpi_enreg%paral_kgb,tim_fourwf0,weight1,weight1,& 219 & use_gpu_cuda=gs_ham%use_gpu_cuda) 220 cwavef_r=cwavef_r*invucvol 221 222 ! ===================================================== 223 ! === Select the states in cgocc_bz with the same spin === 224 ! ===================================================== 225 !* Initialization of the indices/shifts, according to the value of isppol 226 !* bdtot_jindex = shift to be applied on the location of data in the array occ_bz ? 227 bdtot_jindex=0 228 !* jbg = shift to be applied on the location of data in the array cprj/occ 229 jbg=0;jcg=0 230 my_jsppol=fockcommon%isppol 231 if((fockcommon%isppol==2).and.(mpi_enreg%nproc_kpt/=1)) my_jsppol=1 232 233 call timab(1505,2,tsec) 234 call timab(1506,1,tsec) 235 236 !=================================== 237 !=== Loop on the k-points in IBZ === 238 !=================================== 239 jkg=0 240 241 if (associated(gs_ham%ph3d_kp)) then 242 nullify (gs_ham%ph3d_kp) 243 end if 244 245 do jkpt=1,fockbz%mkpt 246 !* nband_k = number of bands at point k_j 247 nband_k=fockbz%nbandocc_bz(jkpt,my_jsppol) 248 !* wtk = weight in BZ of this k point 249 wtk=fockbz%wtk_bz(jkpt) !*sqrt(gs_ham%ucvol) 250 !* jstwfk= how is stored the wavefunction 251 jstwfk=fockbz%istwfk_bz(jkpt) 252 !* npwj= number of plane wave in basis for the wavefunction 253 npwj=fockbz%npwarr(jkpt) 254 !* Basis sphere of G vectors 255 if (allocated(fockbz%cgocc)) then 256 gbound_kp => fockbz%gbound_bz(:,:,jkpt) 257 kg_occ => fockbz%kg_bz(:,1+jkg:npwj+jkg) 258 end if 259 !* Load k^prime hamiltonian in the gs_ham datastructure 260 ! Note: ffnl_kp / ph3d_kp / gbound_kp are not used 261 262 if (associated(gs_ham%ph3d_kp)) then 263 ABI_ALLOCATE(gs_ham%ph3d_kp,(2,npwj,gs_ham%matblk)) 264 end if 265 266 call load_kprime_hamiltonian(gs_ham,kpt_kp=fockbz%kptns_bz(:,jkpt),& 267 & istwf_kp=jstwfk,npw_kp=npwj,kg_kp=fockbz%kg_bz(:,1+jkg:npwj+jkg)) 268 !* Some temporary allocations needed for PAW 269 if (fockcommon%usepaw==1) then 270 ABI_ALLOCATE(vectin_dum,(2,npwj*nspinor)) 271 vectin_dum=zero 272 ABI_ALLOCATE(ffnl_kp_dum,(npwj,0,gs_ham%lmnmax,gs_ham%ntypat)) 273 call load_kprime_hamiltonian(gs_ham,ffnl_kp=ffnl_kp_dum) 274 end if 275 276 ! ====================================== 277 ! === Calculate the vector q=k_i-k_j === 278 ! ====================================== 279 !* Evaluation of kpoint_j, the considered k-point in reduced coordinates 280 ! kpoint_j(:)=fockbz%kptns_bz(:,jkpt) 281 !* the vector qvec is expressed in reduced coordinates. 282 ! qvec(:)=kpoint_i(:)-kpoint_j(:) 283 qvec_j(:)=gs_ham%kpt_k(:)-fockbz%kptns_bz(:,jkpt) 284 call bare_vqg(qvec_j,fockcommon%gsqcut,gs_ham%gmet,fockcommon%usepaw,fockcommon%hyb_mixing,& 285 & fockcommon%hyb_mixing_sr,fockcommon%hyb_range_fock,nfftf,fockbz%nkpt_bz,ngfftf,gs_ham%ucvol,vqg) 286 287 288 289 ! ================================================= 290 ! === Loop on the band indices jband of cgocc_k === 291 ! ================================================= 292 293 do jband=1,nband_k 294 295 !* occ = occupancy of jband at this k point 296 occ=fockbz%occ_bz(jband+bdtot_jindex,my_jsppol) 297 if(occ<tol8) cycle 298 299 ! ============================================== 300 ! === Get cwaveocc_r in real space using FFT === 301 ! ============================================== 302 if (allocated(fockbz%cwaveocc_bz)) then 303 cwaveocc_r => fockbz%cwaveocc_bz(:,:,:,:,jband+jbg,my_jsppol) 304 else 305 ABI_ALLOCATE(cwaveocc_r,(2,n4f,n5f,n6f)) 306 cwaveocc_r=zero 307 call fourwf(1,rhodum0,fockbz%cgocc(:,1+jcg+npwj*(jband-1):jcg+jband*npwj,my_jsppol),rhodum,cwaveocc_r, & 308 & gbound_kp,gbound_kp,jstwfk,kg_occ,kg_occ,mgfftf,mpi_enreg,ndat1,ngfftf,& 309 & npwj,1,n4f,n5f,n6f,tim_fourwf0,mpi_enreg%paral_kgb,0,weight1,weight1,use_gpu_cuda=gs_ham%use_gpu_cuda) 310 cwaveocc_r=cwaveocc_r*invucvol 311 end if 312 313 ! ================================================ 314 ! === Get the overlap density matrix rhor_munu === 315 ! ================================================ 316 !* Calculate the overlap density matrix in real space = conj(cwaveocc_r)*cwavef_r 317 !* rhor_munu will contain the overlap density matrix. 318 ! vfock=-int{conj(cwaveocc_r)*cwavef_r*dr'/|r-r'|} 319 320 call timab(1508,1,tsec) 321 ind=0 322 do i3=1,n3f 323 do i2=1,n2f 324 do i1=1,n1f 325 ind=ind+1 326 recwf =cwavef_r(1,i1,i2,i3) ; imcwf =cwavef_r(2,i1,i2,i3) 327 recwocc=cwaveocc_r(1,i1,i2,i3) ; imcwocc=cwaveocc_r(2,i1,i2,i3) 328 rhor_munu(1,ind)= recwocc*recwf+imcwocc*imcwf 329 rhor_munu(2,ind)= recwocc*imcwf-imcwocc*recwf 330 end do ! i1 331 end do ! i2 332 end do ! i3 333 call timab(1508,2,tsec) 334 335 ! ======================================================= 336 ! === Add compensation charge density in the PAW case === 337 ! === Get the overlap density matrix rhor_munu === 338 ! ======================================================= 339 call timab(1509,1,tsec) 340 if (fockcommon%usepaw==1) then 341 iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+jband 342 343 ABI_ALLOCATE(rho12,(2,nfftf,nspinor**2)) 344 345 cwaveocc_prj=>fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+nspinor-1) 346 347 call pawmknhat_psipsi(cwaveprj,cwaveocc_prj,ider,izero,natom,natom,nfftf,ngfftf,& 348 & nhat12_grdim,nspinor,fockcommon%ntypat,fockbz%pawang,fockcommon%pawfgrtab,grnhat12,rho12,& 349 & fockcommon%pawtab,gprimd=gs_ham%gprimd,grnhat_12=grnhat_12,qphon=qvec_j,xred=gs_ham%xred,atindx=gs_ham%atindx) 350 351 rhor_munu(1,:)=rhor_munu(1,:)+rho12(1,:,nspinor) 352 rhor_munu(2,:)=rhor_munu(2,:)-rho12(2,:,nspinor) 353 end if 354 355 !Perform an FFT using fourwf to get rhog_munu = FFT^-1(rhor_munu) 356 call fourdp(cplex_fock,rhog_munu,rhor_munu,-1,mpi_enreg,nfftf,ngfftf,mpi_enreg%paral_kgb,tim_fourdp0) 357 call timab(1509,2,tsec) 358 359 if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 360 call strfock(gs_ham%gprimd,fockcommon%gsqcut,fockstr,fockcommon%hyb_mixing,fockcommon%hyb_mixing_sr,& 361 & fockcommon%hyb_range_fock,mpi_enreg,nfftf,ngfftf,fockbz%nkpt_bz,rhog_munu,gs_ham%ucvol,qvec_j) 362 fockcommon%stress_ikpt(:,fockcommon%ieigen)=fockcommon%stress_ikpt(:,fockcommon%ieigen)+fockstr(:)*occ*wtk 363 if (fockcommon%usepaw==0.and.(.not.need_ghc)) then 364 if (allocated(fockbz%cgocc)) then 365 ABI_DEALLOCATE(cwaveocc_r) 366 end if 367 cycle 368 end if 369 end if 370 371 ! =================================================== 372 ! === Calculate the local potential vfockloc_munu === 373 ! =================================================== 374 !* Apply the Poisson solver to "rhog_munu" while taking into account the effect of the vector "qvec" 375 !* This is precisely what is done in the subroutine hartre, with option cplex=2. 376 !* vfock will contain the local Fock potential, the result of hartre routine. 377 !* vfock = FFT( rhog_munu/|g+qvec|^2 ) 378 call timab(1510,1,tsec) 379 #if 0 380 381 call hartre(cplex_fock,fockcommon%gsqcut,fockcommon%usepaw,mpi_enreg,nfftf,ngfftf,& 382 & mpi_enreg%paral_kgb,rhog_munu,rprimd,vfock,divgq0=fock%divgq0,qpt=qvec_j) 383 384 #else 385 do ifft=1,nfftf 386 rhog_munu(1,ifft) = rhog_munu(1,ifft) * vqg(ifft) 387 rhog_munu(2,ifft) = rhog_munu(2,ifft) * vqg(ifft) 388 end do 389 390 391 call fourdp(cplex_fock,rhog_munu,vfock,+1,mpi_enreg,nfftf,ngfftf,mpi_enreg%paral_kgb,tim_fourdp0) 392 393 #endif 394 call timab(1510,2,tsec) 395 396 !=============================================================== 397 !======== Calculate Dij_Fock_hat contribution in case of PAW === 398 !=============================================================== 399 400 if (fockcommon%usepaw==1) then 401 qphon=qvec_j;nfftotf=product(ngfftf(1:3)) 402 cplex_dij=cplex_fock;ndij=nspden_fock 403 ABI_ALLOCATE(dijhat,(cplex_dij*gs_ham%dimekb1,natom,ndij)) 404 dijhat=zero 405 do iatom=1,natom 406 ipert=iatom 407 itypat=gs_ham%typat(iatom) 408 lmn2_size=fockcommon%pawtab(itypat)%lmn2_size 409 ABI_ALLOCATE(dijhat_tmp,(cplex_dij*lmn2_size,ndij)) 410 dijhat_tmp=zero 411 call pawdijhat(cplex_fock,cplex_dij,dijhat_tmp,gs_ham%gprimd,iatom,ipert,& 412 !& natom,ndij,nfftf,nfftotf,nspden_fock,my_jsppol,fockbz%pawang,fockcommon%pawfgrtab(iatom),& 413 & natom,ndij,nfftf,nfftotf,nspden_fock,nspden_fock,fockbz%pawang,fockcommon%pawfgrtab(iatom),& 414 & fockcommon%pawtab(itypat),vfock,qphon,gs_ham%ucvol,gs_ham%xred) 415 dijhat(1:cplex_dij*lmn2_size,iatom,:)=dijhat_tmp(1:cplex_dij*lmn2_size,:) 416 ABI_DEALLOCATE(dijhat_tmp) 417 end do 418 signs=2; cpopt=2;idir=0; paw_opt=1;nnlout=1;tim_nonlop=1 419 if(need_ghc) then 420 choice=1 421 call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,& 422 & ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,gvnlc,enl=dijhat,& 423 & select_k=K_H_KPRIME) 424 ghc2=ghc2-gvnlc*occ*wtk 425 end if 426 427 ! Forces calculation 428 429 if (fockcommon%optfor.and.(fockcommon%ieigen/=0)) then 430 choice=2; dotr=zero;doti=zero;cpopt=4 431 432 do iatom=1,natom 433 do idir=1,3 434 call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,& 435 & ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,& 436 & forout,enl=dijhat,iatom_only=iatom,& 437 & select_k=K_H_KPRIME) 438 call dotprod_g(dotr(idir),doti,gs_ham%istwf_k,npw,2,cwavef,forout,mpi_enreg%me_g0,mpi_enreg%comm_fft) 439 for1(idir)=zero 440 do ifft=1,fockcommon%pawfgrtab(iatom)%nfgd 441 ind=fockcommon%pawfgrtab(iatom)%ifftsph(ifft) 442 for1(idir)=for1(idir)+vfock(2*ind-1)*grnhat_12(1,ind,1,idir,iatom)-& 443 & vfock(2*ind)*grnhat_12(2,ind,1,idir,iatom) 444 end do 445 end do 446 do idir=1,3 447 for12(idir)=rprimd(1,idir)*for1(1)+rprimd(2,idir)*for1(2)+rprimd(3,idir)*for1(3) 448 forikpt(idir,iatom)=forikpt(idir,iatom)-(for12(idir)*gs_ham%ucvol/nfftf+dotr(idir))*occ*wtk 449 end do 450 end do 451 end if 452 453 ! Stresses calculation 454 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 455 signs=2;choice=3;cpopt=4 456 457 ! first contribution 458 dotr=zero 459 do idir=1,6 460 call nonlop(choice,cpopt,cwaveocc_prj,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,& 461 & ndat1,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,& 462 & strout,enl=dijhat,select_k=K_H_KPRIME) 463 call dotprod_g(dotr(idir),doti,gs_ham%istwf_k,npw,2,cwavef,strout,mpi_enreg%me_g0,mpi_enreg%comm_fft) 464 fockcommon%stress_ikpt(idir,fockcommon%ieigen)=fockcommon%stress_ikpt(idir,fockcommon%ieigen)-& 465 & dotr(idir)*occ*wtk/gs_ham%ucvol 466 end do 467 ! second contribution 468 str=zero 469 do iatom=1,natom 470 do idir=1,3 471 do idir1=1,3 472 do ifft=1,fockcommon%pawfgrtab(iatom)%nfgd 473 ind=fockcommon%pawfgrtab(iatom)%ifftsph(ifft) 474 str(idir,idir1)=str(idir,idir1)+(vfock(2*ind-1)*grnhat_12(1,ind,1,idir,iatom)-& 475 & vfock(2*ind)*grnhat_12(2,ind,1,idir,iatom))*fockcommon%pawfgrtab(iatom)%rfgd(idir1,ifft) 476 477 end do 478 end do 479 end do 480 end do 481 do idir=1,3 482 fockstr(idir)=str(idir,idir) 483 end do 484 fockstr(4)=(str(3,2)+str(2,3))*half 485 fockstr(5)=(str(3,1)+str(1,3))*half 486 fockstr(6)=(str(1,2)+str(2,1))*half 487 do idir=1,6 488 fockcommon%stress_ikpt(idir,fockcommon%ieigen)=fockcommon%stress_ikpt(idir,fockcommon%ieigen)+& 489 & fockstr(idir)/nfftf*occ*wtk 490 end do 491 492 ! third contribution 493 doti=zero 494 do ifft=1,nfftf 495 doti=doti+vfock(2*ifft-1)*rho12(1,ifft,nspinor)-vfock(2*ifft)*rho12(2,ifft,nspinor) 496 end do 497 fockcommon%stress_ikpt(1:3,fockcommon%ieigen)=fockcommon%stress_ikpt(1:3,fockcommon%ieigen)-doti/nfftf*occ*wtk 498 ! doti=zero 499 ! do ifft=1,nfftf 500 ! doti=doti+vfock(2*ifft-1)*rhor_munu(1,ifft)-vfock(2*ifft)*rhor_munu(2,ifft) 501 ! end do 502 ! fockcommon%stress_ikpt(1:3,fockcommon%ieigen)=fockcommon%stress_ikpt(1:3,fockcommon%ieigen)+doti/nfftf*occ*wtk*half 503 end if ! end stresses 504 505 ABI_DEALLOCATE(dijhat) 506 ABI_DEALLOCATE(rho12) 507 end if !end PAW 508 509 ! ============================================================= 510 ! === Apply the local potential vfockloc_munu to cwaveocc_r === 511 ! ============================================================= 512 call timab(1507,1,tsec) 513 ind=0 514 do i3=1,ngfftf(3) 515 do i2=1,ngfftf(2) 516 do i1=1,ngfftf(1) 517 ind=ind+1 518 ! ind=i1+ngfftf(1)*(i2-1+ngfftf(2)*(i3-1)) 519 revloc=vfock(2*ind-1) ; imvloc=vfock(2*ind) 520 recwocc=cwaveocc_r(1,i1,i2,i3) ; imcwocc=cwaveocc_r(2,i1,i2,i3) 521 vlocpsi_r(2*ind-1)=vlocpsi_r(2*ind-1)-(revloc*recwocc-imvloc*imcwocc)*occ*wtk 522 vlocpsi_r(2*ind )=vlocpsi_r(2*ind )-(revloc*imcwocc+imvloc*recwocc)*occ*wtk 523 end do 524 end do 525 end do 526 call timab(1507,2,tsec) 527 if (allocated(fockbz%cgocc)) then 528 ABI_DEALLOCATE(cwaveocc_r) 529 end if 530 531 end do ! jband 532 533 ! ================================================ 534 ! === End : update of shifts and deallocations === 535 ! ================================================ 536 !* Update of the shifts to be applied (reminder : mkmem is not 0, nspinor=1) 537 jcg=jcg+npwj*nband_k 538 jbg=jbg+nband_k 539 bdtot_jindex=bdtot_jindex+nband_k 540 jkg=jkg+npwj 541 if (fockcommon%usepaw==1) then 542 ABI_DEALLOCATE(vectin_dum) 543 ABI_DEALLOCATE(ffnl_kp_dum) 544 end if 545 if (associated(gs_ham%ph3d_kp)) then 546 ABI_DEALLOCATE(gs_ham%ph3d_kp) 547 end if 548 end do ! jkpt 549 550 if (fockcommon%usepaw==1) then 551 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 552 call xmpi_sum(forikpt,mpi_enreg%comm_hf,ier) 553 do iatom=1,natom !Loop over atom 554 ia=gs_ham%atindx(iatom) 555 fockcommon%forces_ikpt(:,ia,fockcommon%ieigen)=forikpt(:,iatom) 556 end do 557 end if 558 end if 559 if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 560 call xmpi_sum(fockcommon%stress_ikpt,mpi_enreg%comm_hf,ier) 561 end if 562 563 if (.not.need_ghc) then 564 565 ! =============================== 566 ! === Deallocate local arrays === 567 ! =============================== 568 ABI_DEALLOCATE(cwavef_r) 569 ABI_DEALLOCATE(ghc1) 570 ABI_DEALLOCATE(ghc2) 571 ABI_DEALLOCATE(rhor_munu) 572 ABI_DEALLOCATE(rhog_munu) 573 ABI_DEALLOCATE(vlocpsi_r) 574 ABI_DEALLOCATE(dummytab) 575 ABI_DEALLOCATE(vfock) 576 ABI_DEALLOCATE(vqg) 577 if (fockcommon%usepaw==1) then 578 ABI_DEALLOCATE(gvnlc) 579 ABI_DEALLOCATE(grnhat12) 580 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 581 ABI_DEALLOCATE(forikpt) 582 ABI_DEALLOCATE(forout) 583 end if 584 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 585 ABI_DEALLOCATE(strout) 586 end if 587 ABI_DEALLOCATE(grnhat_12) 588 end if 589 if(fockcommon%usepaw==1.or.fockcommon%optstr) then 590 ABI_DEALLOCATE(gboundf) 591 end if 592 !*Restore gs_ham datastructure 593 594 if (associated(gs_ham%ph3d_kp)) then 595 ABI_ALLOCATE(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk)) 596 end if 597 call load_kprime_hamiltonian(gs_ham,kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,& 598 & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k) 599 600 ! if (fockcommon%ieigen/=0) fockcommon%ieigen=0 601 return 602 end if 603 604 605 call timab(1506,2,tsec) 606 call timab(1511,1,tsec) 607 608 !*Restore gs_ham datastructure 609 610 if (associated(gs_ham%ph3d_kp)) then 611 ABI_ALLOCATE(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk)) 612 end if 613 call load_kprime_hamiltonian(gs_ham,kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,& 614 & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k) 615 616 !* Perform an FFT using fourwf to get ghc1 = FFT^-1(vlocpsi_r) 617 ABI_ALLOCATE(psilocal,(cplex_fock*n4f,n5f,n6f)) 618 call fftpac(1,mpi_enreg,nspden_fock,cplex_fock*n1f,n2f,n3f,cplex_fock*n4f,n5f,n6f,ngfft,vlocpsi_r,psilocal,2) 619 620 call fourwf(0,rhodum0,rhodum,ghc1,psilocal,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,& 621 & mgfftf,mpi_enreg,ndat1,ngfftf,1,npw,n4f,n5f,n6f,3,mpi_enreg%paral_kgb,tim_fourwf0,weight1,weight1,& 622 & use_gpu_cuda=gs_ham%use_gpu_cuda) 623 ABI_DEALLOCATE(psilocal) 624 625 ghc1=ghc1*sqrt(gs_ham%ucvol)+ghc2 626 627 !* If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc 628 ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_hf + ghc1(:,:) 629 630 call xmpi_sum(ghc,mpi_enreg%comm_hf,ier) 631 632 call timab(1511,2,tsec) 633 634 635 ! =============================== 636 ! === Deallocate local PAW arrays === 637 ! =============================== 638 639 if (fockcommon%usepaw==1) then 640 ABI_DEALLOCATE(gvnlc) 641 ABI_DEALLOCATE(grnhat12) 642 if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then 643 ABI_DEALLOCATE(forikpt) 644 ABI_DEALLOCATE(forout) 645 end if 646 if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then 647 ABI_DEALLOCATE(strout) 648 end if 649 ABI_DEALLOCATE(grnhat_12) 650 end if 651 if(fockcommon%usepaw==1.or.fockcommon%optstr) then 652 ABI_DEALLOCATE(gboundf) 653 end if 654 655 ! ============================================ 656 ! === Calculate the contribution to energy === 657 ! ============================================ 658 !* Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction 659 !* at each convergence step. (consistent definition with the defintion of hartree energy) 660 if (fockcommon%ieigen/=0) then 661 eigen=zero 662 !* Dot product of cwavef and ghc 663 !* inspired from the routine 53_spacepar/meanvalue_g but without the reference to parallelism and filtering 664 if(gs_ham%istwf_k==2) then 665 eigen=half*cwavef(1,1)*ghc1(1,1) 666 else 667 eigen=cwavef(1,1)*ghc1(1,1)+cwavef(2,1)*ghc1(2,1) 668 end if 669 do ipw=2,npw 670 eigen=eigen+cwavef(1,ipw)*ghc1(1,ipw)+cwavef(2,ipw)*ghc1(2,ipw) 671 end do 672 if(gs_ham%istwf_k>=2) eigen=two*eigen 673 call xmpi_sum(eigen,mpi_enreg%comm_hf,ier) 674 fockcommon%eigen_ikpt(fockcommon%ieigen)= eigen 675 if(fockcommon%use_ACE==0) fockcommon%ieigen = 0 676 end if 677 678 ! =============================== 679 ! === Deallocate local arrays === 680 ! =============================== 681 ABI_DEALLOCATE(cwavef_r) 682 ABI_DEALLOCATE(ghc1) 683 ABI_DEALLOCATE(ghc2) 684 ABI_DEALLOCATE(rhor_munu) 685 ABI_DEALLOCATE(rhog_munu) 686 ABI_DEALLOCATE(vlocpsi_r) 687 ABI_DEALLOCATE(dummytab) 688 ABI_DEALLOCATE(vfock) 689 ABI_DEALLOCATE(vqg) 690 691 call timab(1504,2,tsec) 692 693 end subroutine fock_getghc