TABLE OF CONTENTS
ABINIT/forstr [ Functions ]
NAME
forstr
FUNCTION
Drives the computation of forces and/or stress tensor
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx cg(2,mcg)=wavefunctions (may be read from disk instead of input) cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn> dtefield <type(efield_type)> = variables related to Berry phase dtset <type(dataset_type)>=all input variables in this dataset | berryopt = 4: electric field is on -> add the contribution of the | - \Omega E.P term to the total energy | /= 4: electric field is off | from Etot(npw) data (at fixed geometry), used for making | Pulay correction to stress tensor (hartree). Should be <=0. | ecut=cut-off energy for plane wave basis sphere (Ha) | ecutsm=smearing energy for plane wave kinetic energy (Ha) | effmass_free=effective mass for electrons (1. in common case) | efield = cartesian coordinates of the electric field in atomic units | ionmov=governs the movement of atoms (see help file) | densfor_pred=governs the mixed electronic-atomic part of the preconditioner | istwfk(nkpt)=input option parameter that describes the storage of wfs | kptns(3,nkpt)=reduced coordinates of k points in Brillouin zone | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem=maximum number of k points in core memory | mpw = maximum number of plane waves | natom=number of atoms in cell | nband(nkpt*nsppol)=number of bands to be included in summation at each k point | nfft=(effective) number of FFT grid points (for this processor) | ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft | nkpt=number of k points in Brillouin zone | nloalg(3)=governs the choice of the algorithm for non-local operator. | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | pawprtvol=control print volume and debugging output for PAW | prtvol=integer controlling volume of printed output | symafm(nsym)=(anti)ferromagnetic part of symmetry operations | tfkinfunc=1 if use of Thomas-Fermi kinetic functional | =2 if use of recursion method | typat(natom)=type integer for each atom in cell | wtk(nkpt)=weights associated with various k points | nsym=number of symmetries in space group energies <type(energies_type)>=all part of total energy. | e_localpsp(IN)=local psp energy (hartree) | e_hartree(IN)=Hartree part of total energy (hartree units) | e_corepsp(IN)=psp core-core energy | e_kinetic(IN)=kinetic energy part of total energy. eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) fock <type(fock_type)>= quantities to calculate Fock exact exchange grchempottn(3,natom)=d(E_chemical potential)/d(xred) (hartree) grcondft(3,natom)=d(E_constrainedDFT)/d(xred) (hartree) grewtn(3,natom)=d(Ewald)/d(xred) (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree) gsqcut=cutoff value on G**2 for (large) sphere inside FFT box. gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2) indsym(4,nsym,natom)=index showing transformation of atom labels under symmetry operations (computed in symatm) kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfftf= -PAW ONLY- maximum size of 1D FFTs for the fine grid (mgfftf=mgfft for norm-conserving potential runs) mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor n3xccc=dimension of the xccc3d array (0 or nfftf). nattyp(ntypat)=number of atoms of each type nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs) ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid (ngfftf=ngfft for norm-conserving potential runs) ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwarr(nkpt)=number of planewaves in basis and on boundary for each k ntypat=number of types of atoms nvresid(nfftf,nspden)=array for the residual of the density/potential occ(mband*nkpt*nsppol)=occupancies of bands at various k points optfor=1 if computation of forces is required optres=0 if the potential residual has to be used for forces corrections =1 if the density residual has to be used for forces corrections paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases ph1df(2,3*(2*mgfftf+1)*natom)=-PAW only- 1-dim structure factor phases for the fine grid psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum rhog(2,nfftf)=Fourier transform of charge density (bohr^-3) rhor(nfftf,nspden)=array for electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) strscondft(6)=cDFT correction to stress strsxc(6)=xc correction to stress stress_needed=1 if computation of stress tensor is required symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates ucvol=unit cell volume in bohr**3 usecprj=1 if cprj datastructure is stored in memory vhartr(nfftf)=array for holding Hartree potential vpsp(nfftf)=array for holding local psp vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space vxctau(nfft,nspden,4*usekden)=(only for meta-GGA): derivative of XC energy density wrt kinetic energy density (depsxcdtau) xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xcctau3d(n3xccc*usekden)=(only for meta-GGA): 3D core electron kinetic energy density for XC core correction xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
OUTPUT
==== if (optfor==1) ==== diffor=maximal absolute value of changes in the components of force between the input and the output. favg(3)=mean of the forces before correction for translational symmetry fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr) at input, previous value of forces, at output, new value. Note : unlike gred, this array has been corrected by enforcing the translational symmetry, namely that the sum of force on all atoms is zero. forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr) gred(3,natom)=symmetrized grtn = d(etotal)/d(xred) gresid(3,natom)=forces due to the residual of the density/potential grhf(3,natom)=Hellman-Feynman derivatives of the total energy grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used maxfor=maximal absolute value of the output array force. synlgr(3,natom)=symmetrized gradients of energy due to nonlocal contributions ==== if (stress_needed==1) ==== strten(6)=components of the stress tensor (hartree/bohr^3) for the 6 unique components of this symmetric 3x3 tensor: Given in order (1,1), (2,2), (3,3), (3,2), (3,1), (2,1).
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr) ===== if psps%usepaw==1 pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data (gradients of rhoij for each atom with respect to atomic positions are computed here) wvl <type(wvl_data)>=all wavelets data.
NOTES
Be careful to the meaning of nfft (size of FFT grids): - In case of norm-conserving calculations the FFT grid is the usual FFT grid. - In case of PAW calculations: Two FFT grids are used; one with nfft points (coarse grid) for the computation of wave functions ; one with nfftf points (fine grid) for the computation of total density.
SOURCE
245 subroutine forstr(atindx1,cg,cprj,diffor,dtefield,dtset,eigen,electronpositron,energies,favg,fcart,fock,& 246 & forold,gred,grchempottn,grcondft,gresid,grewtn,grhf,grvdw,grxc,gsqcut,extfpmd,indsym,& 247 & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,& 248 & nfftf,ngfftf,ngrvdw,nhat,nkxc,npwarr,& 249 & ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,& 250 & pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,psps,rhog,rhor,rprimd,stress_needed,& 251 & strscondft,strsxc,strten,symrec,synlgr,ucvol,usecprj,vhartr,vpsp,& 252 & vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero) 253 254 !Arguments ------------------------------------ 255 !scalars 256 integer,intent(in) :: mcg,mcprj,mgfftf,my_natom,n3xccc,nfftf,ngrvdw,nkxc,ntypat,optfor,optres 257 integer,intent(in) :: stress_needed,usecprj 258 real(dp),intent(in) :: gsqcut,qvpotzero,ucvol 259 real(dp),intent(inout) :: diffor,maxfor 260 type(electronpositron_type),pointer :: electronpositron 261 type(MPI_type),intent(inout) :: mpi_enreg 262 type(efield_type),intent(in) :: dtefield 263 type(dataset_type),intent(in) :: dtset 264 type(energies_type),intent(in) :: energies 265 type(extfpmd_type),pointer,intent(inout) :: extfpmd 266 type(pawang_type),intent(in) :: pawang 267 type(pawfgr_type),intent(in) :: pawfgr 268 type(pseudopotential_type),intent(in) :: psps 269 type(wvl_data),intent(inout) :: wvl 270 type(fock_type),pointer, intent(inout) :: fock 271 !arrays 272 integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 273 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(ntypat),ngfftf(18) 274 integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym) 275 real(dp),intent(in) :: cg(2,mcg) 276 real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 277 real(dp),intent(in) :: grchempottn(3,dtset%natom),grcondft(3,dtset%natom),grewtn(3,dtset%natom) 278 real(dp),intent(in) :: grvdw(3,ngrvdw),kxc(dtset%nfft,nkxc) 279 real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 280 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 281 real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom) 282 real(dp),intent(in) :: rhog(2,nfftf),strscondft(6),strsxc(6),vhartr(nfftf) 283 real(dp),intent(in) :: vpsp(nfftf),vxc(nfftf,dtset%nspden),vxctau(nfftf,dtset%nspden,4*dtset%usekden) 284 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 285 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 286 real(dp),intent(inout) :: forold(3,dtset%natom) 287 real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw),rhor(nfftf,dtset%nspden),rprimd(3,3) 288 real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom) 289 real(dp),intent(inout),target :: nvresid(nfftf,dtset%nspden) 290 real(dp),intent(out) :: favg(3) 291 real(dp),intent(inout) :: fcart(3,dtset%natom),gred(3,dtset%natom) 292 real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom) 293 real(dp),intent(inout) :: grxc(3,dtset%natom),strten(6),synlgr(3,dtset%natom) 294 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj) 295 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 296 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 297 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 298 type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw) 299 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 300 301 !Local variables------------------------------- 302 !scalars 303 304 integer :: comm_grid,ifft,ispden,ncpgr,occopt_,optgr,optgr2,option,optnc,optstr,optstr2,iorder_cprj,ctocprj_choice 305 integer :: idir,iatom,unpaw,mcgbz 306 integer,allocatable :: dimcprj(:) 307 real(dp) ::dum,dum1,ucvol_ 308 logical :: apply_residual 309 !arrays 310 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 311 real(dp) :: kinstr(6),nlstr(6),tsec(2),strdum(6),gmet(3,3),gprimd(3,3),rmet(3,3) 312 real(dp) :: dummy(0) 313 real(dp),allocatable :: grnl(:),vlocal(:,:),vxc_hf(:,:),xcart(:,:),ylmbz(:,:),ylmgrbz(:,:,:) 314 real(dp), ABI_CONTIGUOUS pointer :: resid(:,:) 315 316 ! ************************************************************************* 317 318 call timab(910,1,tsec) 319 call timab(911,1,tsec) 320 321 !Do nothing if nothing is required 322 if (optfor==0.and.stress_needed==0) return 323 324 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW) 325 if (dtset%usewvl==0) then 326 if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then 327 ABI_BUG(' wrong values for nfft, nfftf !') 328 end if 329 if ((psps%usepaw==1.and.pawfgr%mgfft/=mgfftf).or.(psps%usepaw==0.and.dtset%mgfft/=mgfftf)) then 330 ABI_BUG('wrong values for mgfft, mgfftf!') 331 end if 332 end if 333 334 !========================================================================== 335 !Here compute terms common to forces and stresses 336 !========================================================================== 337 338 !output only if (optfor==1) but we have to allocate it 339 ABI_MALLOC(grnl,(3*dtset%natom*optfor)) 340 grnl(:)=zero 341 342 !Compute nonlocal psp + potential Fock ACE parts of forces and stress tensor 343 !-involves summations over wavefunctions at all k points 344 if (dtset%tfkinfunc>0.and.stress_needed==1) then 345 kinstr(1:3)=-two/three*energies%e_kinetic/ucvol ; kinstr(4:6)=zero 346 nlstr(1:6)=zero 347 else if (dtset%usewvl==0) then 348 occopt_=0 ! This means that occ are now fixed 349 if(dtset%usefock==1 .and. associated(fock)) then 350 ! if((dtset%optstress/=0).and.(psps%usepaw==1)) then 351 if((psps%usepaw==1).and.((dtset%optstress/=0).or.(dtset%optforces==2))) then 352 if(dtset%optstress==0) then 353 ctocprj_choice=2 354 ncpgr=3 355 end if 356 if(dtset%optstress/=0) then 357 ctocprj_choice=20*optfor+3*dtset%optstress 358 ncpgr=6*dtset%optstress+3*optfor 359 end if 360 if (allocated(fock%fock_BZ%cwaveocc_prj)) then 361 call pawcprj_free(fock%fock_BZ%cwaveocc_prj) 362 ABI_FREE(fock%fock_BZ%cwaveocc_prj) 363 ABI_MALLOC(fock%fock_BZ%cwaveocc_prj,(dtset%natom,fock%fock_BZ%mcprj)) 364 ABI_MALLOC(dimcprj,(dtset%natom)) 365 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 366 call pawcprj_alloc(fock%fock_BZ%cwaveocc_prj,ncpgr,dimcprj) 367 ABI_FREE(dimcprj) 368 end if 369 iatom=-1;idir=0;iorder_cprj=0;unpaw=26 370 call metric(gmet,gprimd,-1,rmet,rprimd,dum) 371 if (fock%fock_BZ%mkpt/=dtset%mkmem.or.(fock%fock_BZ%mpi_enreg%paral_hf ==1)) then 372 ABI_MALLOC(ylmbz,(dtset%mpw*fock%fock_BZ%mkpt,psps%mpsang*psps%mpsang*psps%useylm)) 373 ABI_MALLOC(ylmgrbz,(dtset%mpw*fock%fock_BZ%mkpt,3,psps%mpsang*psps%mpsang*psps%useylm)) 374 option=1; mcgbz=dtset%mpw*fock%fock_BZ%mkptband*fock%fock_common%my_nsppol 375 call initylmg(gprimd,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,& 376 & psps%mpsang,dtset%mpw,fock%fock_BZ%nbandocc_bz,fock%fock_BZ%mkpt,& 377 & fock%fock_BZ%npwarr,dtset%nsppol,option,rprimd,ylmbz,ylmgrbz) 378 call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,& 379 & iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcgbz,& 380 & fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,fock%fock_BZ%mpi_enreg,psps%mpsang,& 381 & dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,& 382 & dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,& 383 & dtset%nsppol,fock%fock_common%my_nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,& 384 & xred,ylmbz,ylmgrbz) 385 ABI_FREE(ylmbz) 386 ABI_FREE(ylmgrbz) 387 else 388 call ctocprj(fock%fock_common%atindx,fock%fock_BZ%cgocc,ctocprj_choice,fock%fock_BZ%cwaveocc_prj,gmet,gprimd,iatom,idir,& 389 & iorder_cprj,fock%fock_BZ%istwfk_bz,fock%fock_BZ%kg_bz,fock%fock_BZ%kptns_bz,mcg,& 390 & fock%fock_BZ%mcprj,dtset%mgfft,fock%fock_BZ%mkpt,mpi_enreg,psps%mpsang,& 391 & dtset%mpw,dtset%natom,nattyp,fock%fock_BZ%nbandocc_bz,dtset%natom,dtset%ngfft,fock%fock_BZ%mkpt,& 392 & dtset%nloalg,fock%fock_BZ%npwarr,dtset%nspinor,& 393 & dtset%nsppol,fock%fock_common%my_nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,unpaw,& 394 & xred,ylm,ylmgr) 395 end if 396 end if 397 end if 398 call forstrnps(cg,cprj,dtset%ecut,dtset%ecutsm,dtset%effmass_free,eigen,electronpositron,fock,grnl,& 399 & dtset%istwfk,kg,kinstr,nlstr,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,dtset%mkmem,& 400 & mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,& 401 & dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,& 402 & dtset%nucdipmom,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,stress_needed,symrec,dtset%typat,& 403 & usecprj,dtset%usefock,dtset%use_gpu_cuda,dtset%wtk,xred,ylm,ylmgr) 404 !DEBUG 405 ! write(6,*)' after forstrnps, nlstr=',nlstr(1:6) 406 !ENDDEBUG 407 else if (optfor>0) then !WVL 408 ABI_MALLOC(xcart,(3, dtset%natom)) 409 call xred2xcart(dtset%natom, rprimd, xcart, xred) 410 call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart) 411 ABI_FREE(xcart) 412 end if 413 414 call timab(911,2,tsec) 415 call timab(912,1,tsec) 416 417 !PAW: add gradients due to Dij derivatives to non-local term 418 if (psps%usepaw==1) then 419 ABI_MALLOC(vlocal,(nfftf,dtset%nspden)) 420 421 !$OMP PARALLEL DO COLLAPSE(2) 422 do ispden=1,min(dtset%nspden,2) 423 do ifft=1,nfftf 424 vlocal(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft) 425 end do 426 end do 427 if (dtset%nspden==4) then 428 !$OMP PARALLEL DO COLLAPSE(2) 429 do ispden=3,4 430 do ifft=1,nfftf 431 vlocal(ifft,ispden)=vxc(ifft,ispden) 432 end do 433 end do 434 end if 435 ucvol_=ucvol 436 #if defined HAVE_BIGDFT 437 if (dtset%usewvl==1) ucvol_=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp) 438 #endif 439 optgr=optfor;optgr2=0;optstr=stress_needed;optstr2=0 440 comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl 441 call pawgrnl(atindx1,dtset%nspden,dummy,1,dummy,grnl,gsqcut,mgfftf,my_natom,dtset%natom,& 442 & nattyp,nfftf,ngfftf,nhat,nlstr,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,& 443 & pawang,pawfgrtab,pawrhoij,pawtab,ph1df,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred,& 444 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=comm_grid) 445 !DEBUG 446 ! write(6,*)' after pawgrnl, nlstr=',nlstr(1:6) 447 !ENDDEBUG 448 ABI_FREE(vlocal) 449 450 end if 451 call timab(912,2,tsec) 452 call timab(913,1,tsec) 453 454 !========================================================================== 455 !Here compute forces (if required) 456 !========================================================================== 457 if (optfor==1) then 458 apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. & 459 & abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) 460 461 ! If residual is a density residual (and forces from residual asked), 462 ! has to convert it into a potential residual before calling forces routine 463 if (apply_residual) then 464 ABI_MALLOC(resid,(nfftf,dtset%nspden)) 465 option=0; if (dtset%densfor_pred<0) option=1 466 optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2 467 call nres2vres(dtset,gsqcut,psps%usepaw,kxc,mpi_enreg,my_natom,nfftf,ngfftf,nhat,& 468 & nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,& 469 & rhor,rprimd,psps%usepaw,resid,xccc3d,xred,vxc) 470 else 471 resid => nvresid 472 end if 473 474 call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,gred,grchempottn,grcondft,gresid,grewtn,& 475 & grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfftf,& 476 & mpi_enreg,psps%n1xccc,n3xccc,nattyp,& 477 & nfftf,ngfftf,ngrvdw,ntypat,pawrad,pawtab,ph1df,psps,rhog,& 478 & rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,vxctau,wvl%descr,wvl%den,xred,& 479 & electronpositron=electronpositron) 480 481 if (apply_residual) then 482 ABI_FREE(resid) 483 end if 484 end if 485 486 call timab(913,2,tsec) 487 call timab(914,1,tsec) 488 489 !========================================================================== 490 !Here compute stress tensor (if required) 491 !========================================================================== 492 493 if (stress_needed==1.and.dtset%usewvl==0) then 494 ! if (dtset%usefock==1 .and. associated(fock).and.fock%fock_common%optstr.and.psps%usepaw==0) then 495 if (dtset%usefock==1 .and. associated(fock)) then 496 if (fock%fock_common%optstr) then 497 fock%fock_common%stress(1:3)=fock%fock_common%stress(1:3)-(two*energies%e_fock-energies%e_fock0)/ucvol 498 if (n3xccc>0.and.psps%usepaw==0 .and. & 499 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) then 500 ABI_MALLOC(vxc_hf,(nfftf,dtset%nspden)) 501 !compute Vxc^GGA(rho_val) 502 call xchybrid_ncpp_cc(dtset,dum,mpi_enreg,nfftf,ngfftf,n3xccc,rhor,rprimd,strdum,dum1,xccc3d,vxc=vxc_hf,optstr=1) 503 end if 504 end if 505 end if 506 call stress(atindx1,dtset%berryopt,dtefield,energies%e_localpsp,dtset%efield,& 507 & energies%e_hartree,energies%e_corepsp,fock,gsqcut,extfpmd,dtset%ixc,kinstr,mgfftf,& 508 & mpi_enreg,psps%mqgrid_vl,psps%n1xccc,n3xccc,dtset%natom,nattyp,& 509 & nfftf,ngfftf,nlstr,dtset%nspden,dtset%nsym,ntypat,psps,pawrad,pawtab,ph1df,& 510 & dtset%prtvol,psps%qgrid_vl,dtset%red_efieldbar,rhog,rprimd,strten,strscondft,strsxc,symrec,& 511 & dtset%typat,dtset%usefock,dtset%usekden,psps%usepaw,& 512 & dtset%vdw_tol,dtset%vdw_tol_3bt,dtset%vdw_xc,psps%vlspl,vxc,vxctau,vxc_hf,& 513 & psps%xccc1d,xccc3d,xcctau3d,psps%xcccrc,xred,psps%ziontypat,psps%znucltypat,qvpotzero,& 514 & electronpositron=electronpositron) 515 end if 516 517 !Memory deallocation 518 ABI_FREE(grnl) 519 if (allocated(vxc_hf)) then 520 ABI_FREE(vxc_hf) 521 end if 522 523 524 call timab(914,2,tsec) 525 call timab(910,2,tsec) 526 527 end subroutine forstr
ABINIT/forstrnps [ Functions ]
NAME
forstrnps
FUNCTION
Compute nonlocal pseudopotential energy contribution to forces and/or stress tensor as well as kinetic energy contribution to stress tensor.
INPUTS
cg(2,mcg)=wavefunctions (may be read from disk file) cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn> ecut=cut-off energy for plane wave basis sphere (Ha) ecutsm=smearing energy for plane wave kinetic energy (Ha) effmass_free=effective mass for electrons (1. in common case) electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) fock <type(fock_type)>= quantities to calculate Fock exact exchange istwfk(nkpt)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis kpt(3,nkpt)=k points in reduced coordinates mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mkmem=number of k points treated by this node. mpi_enreg=information about MPI parallelization mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mpw= maximum number of plane waves my_natom=number of atoms treated by current processor natom=number of atoms in cell. nband(nkpt)=number of bands at each k point nfft=number of FFT grid points ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpt=number of k points in Brillouin zone nloalg(3)=governs the choice of the algorithm for non-local operator. npwarr(nkpt)=number of planewaves in basis and boundary at each k nspden=Number of spin Density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of elements in symmetry group ntypat=number of types of atoms nucdipmom(3,my_natom)= nuclear dipole moments occ(mband*nkpt*nsppol)=occupation numbers for each band over all k points optfor=1 if computation of forces is required paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations in real space (bohr) stress_needed=1 if computation of stress tensor is required symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless) typat(natom)=type of each atom usecprj=1 if cprj datastructure has been allocated use_gpu_cuda= 0 or 1 to know if we use cuda for nonlop call wtk(nkpt)=weight associated with each k point xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
OUTPUT
if (optfor==1) grnl(3*natom*optfor)=stores grads of nonlocal energy wrt atomic coordinates if (stress_needed==1) kinstr(6)=kinetic energy part of stress tensor (hartree/bohr^3) Store 6 unique components of symmetric 3x3 tensor in the order 11, 22, 33, 32, 31, 21 npsstr(6)=nonlocal pseudopotential energy part of stress tensor (hartree/bohr^3)
SOURCE
601 subroutine forstrnps(cg,cprj,ecut,ecutsm,effmass_free,eigen,electronpositron,fock,& 602 & grnl,istwfk,kg,kinstr,npsstr,kpt,mband,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,& 603 & mpw,my_natom,natom,nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,nsym,& 604 & ntypat,nucdipmom,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,& 605 & stress_needed,symrec,typat,usecprj,usefock,use_gpu_cuda,wtk,xred,ylm,ylmgr) 606 607 !Arguments ------------------------------------ 608 !scalars 609 integer,intent(in) :: mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,nfft,nkpt 610 integer,intent(in) :: nspden,nsppol,nspinor,nsym,ntypat,optfor,stress_needed 611 integer,intent(in) :: usecprj,usefock,use_gpu_cuda 612 real(dp),intent(in) :: ecut,ecutsm,effmass_free 613 type(electronpositron_type),pointer :: electronpositron 614 type(MPI_type),intent(inout) :: mpi_enreg 615 type(pseudopotential_type),intent(in) :: psps 616 !arrays 617 integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol) 618 integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt) 619 integer,intent(in) :: symrec(3,3,nsym),typat(natom) 620 real(dp),intent(in) :: cg(2,mcg) 621 real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kpt(3,nkpt),nucdipmom(3,my_natom) 622 real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom) 623 real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom) 624 real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm) 625 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm) 626 real(dp),intent(out) :: grnl(3*natom*optfor),kinstr(6),npsstr(6) 627 type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj) 628 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 629 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 630 type(fock_type),pointer, intent(inout) :: fock 631 !Local variables------------------------------- 632 !scalars 633 integer,parameter :: tim_rwwf=7 634 integer :: bandpp,bdtot_index,choice,cpopt,dimffnl,dimffnl_str,iband,iband_cprj,iband_last,ibg,icg,ider,ider_str 635 integer :: idir,idir_str,ierr,ii,ikg,ikpt,ilm,ipositron,ipw,ishift,isppol,istwf_k 636 integer :: mband_cprj,me_distrb,my_ikpt,my_nspinor,nband_k,nband_cprj_k,ndat,nkpg 637 integer :: nnlout,npw_k,paw_opt,signs,spaceComm 638 integer :: tim_nonlop,tim_nonlop_prep,usecprj_local,use_ACE_old 639 integer :: blocksize,iblock,iblocksize,ibs,nblockbd 640 real(dp) :: ar,renorm_factor,dfsm,ecutsm_inv,fact_kin,fsm,htpisq,kgc1 641 real(dp) :: kgc2,kgc3,kin,xx 642 type(gs_hamiltonian_type) :: gs_hamk 643 logical :: compute_gbound,usefock_loc 644 character(len=500) :: msg 645 type(fock_common_type),pointer :: fockcommon 646 !arrays 647 integer,allocatable :: kg_k(:,:) 648 real(dp) :: kpoint(3),nonlop_dum(1,1),rmet(3,3),tsec(2) 649 real(dp),allocatable :: cwavef(:,:),enlout(:),ffnl_sav(:,:,:,:),ffnl_str(:,:,:,:) 650 real(dp),allocatable :: ghc_dum(:,:),gprimd(:,:),kpg_k(:,:),kpg_k_sav(:,:) 651 real(dp),allocatable :: kstr1(:),kstr2(:),kstr3(:),kstr4(:),kstr5(:),kstr6(:) 652 real(dp),allocatable :: lambda(:),occblock(:),ph3d(:,:,:),ph3d_sav(:,:,:) 653 real(dp),allocatable :: weight(:),ylm_k(:,:),ylmgr_k(:,:,:) 654 real(dp),allocatable,target :: ffnl(:,:,:,:) 655 type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null() 656 type(pawcprj_type),target,allocatable :: cwaveprj(:,:) 657 type(pawcprj_type),pointer :: cwaveprj_idat(:,:) 658 659 660 !************************************************************************* 661 662 call timab(920,1,tsec) ; call timab(921,-1,tsec) 663 664 !Init mpicomm and me 665 if(mpi_enreg%paral_kgb==1)then 666 spaceComm=mpi_enreg%comm_kpt 667 me_distrb=mpi_enreg%me_kpt 668 else 669 !* In case of HF calculation 670 if (mpi_enreg%paral_hf==1) then 671 spaceComm=mpi_enreg%comm_kpt 672 me_distrb=mpi_enreg%me_kpt 673 else 674 spaceComm=mpi_enreg%comm_cell 675 me_distrb=mpi_enreg%me_cell 676 end if 677 end if 678 679 !Some constants 680 ipositron=abs(electronpositron_calctype(electronpositron)) 681 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 682 !Smearing of plane wave kinetic energy 683 ecutsm_inv=zero;if( ecutsm>1.0d-20) ecutsm_inv=1/ecutsm 684 !htpisq is (1/2) (2 Pi) **2: 685 htpisq=0.5_dp*(two_pi)**2 686 687 !Check that fock is present if want to use fock option 688 compute_gbound=.false. 689 usefock_loc = (usefock==1 .and. associated(fock)) 690 !Arrays initializations 691 grnl(:)=zero 692 if (usefock_loc) then 693 fockcommon =>fock%fock_common 694 fockcommon%optfor=.false. 695 fockcommon%optstr=.false. 696 use_ACE_old=fockcommon%use_ACE 697 fockcommon%use_ACE=0 698 if (fockcommon%optfor) compute_gbound=.true. 699 end if 700 if (stress_needed==1) then 701 kinstr(:)=zero;npsstr(:)=zero 702 if (usefock_loc) then 703 fockcommon%optstr=.TRUE. 704 fockcommon%stress=zero 705 compute_gbound=.true. 706 end if 707 end if 708 709 usecprj_local=usecprj 710 711 if ((usefock_loc).and.(psps%usepaw==1)) then 712 usecprj_local=1 713 if(optfor==1)then 714 fockcommon%optfor=.true. 715 if (.not.allocated(fockcommon%forces_ikpt)) then 716 ABI_MALLOC(fockcommon%forces_ikpt,(3,natom,mband)) 717 end if 718 if (.not.allocated(fockcommon%forces)) then 719 ABI_MALLOC(fockcommon%forces,(3,natom)) 720 end if 721 fockcommon%forces=zero 722 compute_gbound=.true. 723 end if 724 end if 725 726 !Initialize Hamiltonian (k-independent terms) 727 728 call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,nsppol,nspden,natom,& 729 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj_local,& 730 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 731 & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,fock=fock,& 732 & nucdipmom=nucdipmom,use_gpu_cuda=use_gpu_cuda) 733 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 734 735 !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted) 736 if (psps%usepaw==1.and.usecprj_local==1) then 737 call pawcprj_reorder(cprj,gs_hamk%atindx) 738 end if 739 740 !Common data for "nonlop" routine 741 signs=1 ; idir=0 ; ishift=0 ; tim_nonlop=4 ; tim_nonlop_prep=12 742 choice=2*optfor;if (stress_needed==1) choice=10*choice+3 743 if (optfor==1.and.stress_needed==1) ishift=6 744 nnlout=max(1,6*stress_needed+3*natom*optfor) 745 if (psps%usepaw==0) then 746 paw_opt=0 ; cpopt=-1 747 else 748 paw_opt=2 ; cpopt=-1+3*usecprj_local 749 end if 750 751 call timab(921,2,tsec) 752 753 !LOOP OVER SPINS 754 bdtot_index=0;ibg=0;icg=0 755 do isppol=1,nsppol 756 757 ! Continue to initialize the Hamiltonian (PAW DIJ coefficients) 758 call gs_hamk%load_spin(isppol,with_nonlocal=.true.) 759 if (usefock_loc) fockcommon%isppol=isppol 760 761 ! Loop over k points 762 ikg=0 763 do ikpt=1,nkpt 764 if (usefock_loc) fockcommon%ikpt=ikpt 765 nband_k=nband(ikpt+(isppol-1)*nkpt) 766 istwf_k=istwfk(ikpt) 767 npw_k=npwarr(ikpt) 768 kpoint(:)=kpt(:,ikpt) 769 770 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 771 bdtot_index=bdtot_index+nband_k 772 cycle 773 end if 774 775 call timab(922,1,tsec) 776 777 ! Parallelism over FFT and/or bands: define sizes and tabs 778 if (mpi_enreg%paral_kgb==1) then 779 my_ikpt=mpi_enreg%my_kpttab(ikpt) 780 nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp) 781 bandpp=mpi_enreg%bandpp 782 my_bandfft_kpt => bandfft_kpt(my_ikpt) 783 else 784 my_ikpt=ikpt 785 bandpp=mpi_enreg%bandpp 786 nblockbd=nband_k/bandpp 787 end if 788 blocksize=nband_k/nblockbd 789 mband_cprj=mband/mpi_enreg%nproc_band 790 nband_cprj_k=nband_k/mpi_enreg%nproc_band 791 792 ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize)) 793 if (psps%usepaw==1.and.usecprj_local==1) then 794 ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp)) 795 call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj) 796 else 797 ABI_MALLOC(cwaveprj,(0,0)) 798 end if 799 800 if (stress_needed==1) then 801 ABI_MALLOC(kstr1,(npw_k)) 802 ABI_MALLOC(kstr2,(npw_k)) 803 ABI_MALLOC(kstr3,(npw_k)) 804 ABI_MALLOC(kstr4,(npw_k)) 805 ABI_MALLOC(kstr5,(npw_k)) 806 ABI_MALLOC(kstr6,(npw_k)) 807 end if 808 809 ABI_MALLOC(kg_k,(3,mpw)) 810 !$OMP PARALLEL DO 811 do ipw=1,npw_k 812 kg_k(:,ipw)=kg(:,ipw+ikg) 813 end do 814 815 ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 816 if (stress_needed==1) then 817 ABI_MALLOC(ylmgr_k,(npw_k,3,mpsang*mpsang*psps%useylm)) 818 else 819 ABI_MALLOC(ylmgr_k,(0,0,0)) 820 end if 821 if (psps%useylm==1) then 822 !$OMP PARALLEL DO COLLAPSE(2) 823 do ilm=1,mpsang*mpsang 824 do ipw=1,npw_k 825 ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm) 826 end do 827 end do 828 if (stress_needed==1) then 829 !$OMP PARALLEL DO COLLAPSE(2) 830 do ilm=1,mpsang*mpsang 831 do ii=1,3 832 do ipw=1,npw_k 833 ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm) 834 end do 835 end do 836 end do 837 end if 838 end if 839 840 ! Prepare kinetic contribution to stress tensor (Warning : the symmetry 841 ! has not been broken, like in mkkin.f or kpg3.f . It should be, in order to be coherent). 842 if (stress_needed==1) then 843 ABI_MALLOC(gprimd,(3,3)) 844 gprimd=gs_hamk%gprimd 845 !$OMP PARALLEL DO PRIVATE(fact_kin,ipw,kgc1,kgc2,kgc3,kin,xx,fsm,dfsm) & 846 !$OMP&SHARED(ecut,ecutsm,ecutsm_inv,gs_hamk,htpisq,kg_k,kpoint,kstr1,kstr2,kstr3,kstr4,kstr5,kstr6,npw_k) 847 do ipw=1,npw_k 848 ! Compute Cartesian coordinates of (k+G) 849 kgc1=gprimd(1,1)*(kpoint(1)+kg_k(1,ipw))+& 850 & gprimd(1,2)*(kpoint(2)+kg_k(2,ipw))+& 851 & gprimd(1,3)*(kpoint(3)+kg_k(3,ipw)) 852 kgc2=gprimd(2,1)*(kpoint(1)+kg_k(1,ipw))+& 853 & gprimd(2,2)*(kpoint(2)+kg_k(2,ipw))+& 854 & gprimd(2,3)*(kpoint(3)+kg_k(3,ipw)) 855 kgc3=gprimd(3,1)*(kpoint(1)+kg_k(1,ipw))+& 856 & gprimd(3,2)*(kpoint(2)+kg_k(2,ipw))+& 857 & gprimd(3,3)*(kpoint(3)+kg_k(3,ipw)) 858 kin=htpisq* ( kgc1**2 + kgc2**2 + kgc3**2 ) 859 fact_kin=1.0_dp 860 if(kin>ecut-ecutsm)then 861 if(kin>ecut)then 862 fact_kin=0.0_dp 863 else 864 ! See the routine mkkin.f, for the smearing procedure 865 xx=(ecut-kin)*ecutsm_inv 866 ! This kinetic cutoff smoothing function and its xx derivatives 867 ! were produced with Mathematica and the fortran code has been 868 ! numerically checked against Mathematica. 869 fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx)))) 870 dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2 871 ! d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*& 872 ! & (-144+45*xx))))))*fsm**3 873 fact_kin=fsm+kin*(-ecutsm_inv)*dfsm 874 end if 875 end if 876 kstr1(ipw)=fact_kin*kgc1*kgc1 877 kstr2(ipw)=fact_kin*kgc2*kgc2 878 kstr3(ipw)=fact_kin*kgc3*kgc3 879 kstr4(ipw)=fact_kin*kgc3*kgc2 880 kstr5(ipw)=fact_kin*kgc3*kgc1 881 kstr6(ipw)=fact_kin*kgc2*kgc1 882 end do ! ipw 883 ABI_FREE(gprimd) 884 end if 885 886 ! Compute (k+G) vectors 887 nkpg=3*nloalg(3) 888 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 889 if (nkpg>0) then 890 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 891 end if 892 893 ! Compute nonlocal form factors ffnl at all (k+G) 894 ider=0;idir=0;dimffnl=1 895 if (stress_needed==1) then 896 ider=1;dimffnl=2+2*psps%useylm 897 end if 898 ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 899 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 900 & ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,& 901 & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 902 if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then 903 ider_str=1; dimffnl_str=7;idir_str=-7 904 ABI_MALLOC(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,ntypat)) 905 call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 906 & ider_str,idir_str,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,& 907 & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 908 end if 909 910 ! Load k-dependent part in the Hamiltonian datastructure 911 ! - Compute 3D phase factors 912 ! - Prepare various tabs in case of band-FFT parallelism 913 ! - Load k-dependent quantities in the Hamiltonian 914 ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk)) 915 call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,& 916 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.) 917 918 ! Load band-FFT tabs (transposed k-dependent arrays) 919 if (mpi_enreg%paral_kgb==1) then 920 call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav) 921 call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg) 922 call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, & 923 & kg_k =my_bandfft_kpt%kg_k_gather, & 924 & kpg_k =my_bandfft_kpt%kpg_k_gather, & 925 ffnl_k =my_bandfft_kpt%ffnl_gather, & 926 ph3d_k =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound) 927 end if 928 929 ! Setup gemm_nonlop 930 if (gemm_nonlop_use_gemm) then 931 gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt 932 call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax,gs_hamk%ntypat, & 933 & gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, & 934 & gs_hamk%ffnl_k, gs_hamk%ph3d_k, gs_hamk%kpt_k, gs_hamk%kg_k, gs_hamk%kpg_k, & 935 & compute_grad_strain=(stress_needed>0),compute_grad_atom=(optfor>0)) 936 end if 937 938 ! Loop over (blocks of) bands; accumulate forces and/or stresses 939 ! The following is now wrong. In sequential, nblockbd=nband_k/bandpp 940 ! blocksize= bandpp (JB 2016/04/16) 941 ! Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1 942 ABI_MALLOC(lambda,(blocksize)) 943 ABI_MALLOC(occblock,(blocksize)) 944 ABI_MALLOC(weight,(blocksize)) 945 ABI_MALLOC(enlout,(nnlout*blocksize)) 946 occblock=zero;weight=zero;enlout(:)=zero 947 if (usefock_loc) then 948 if (fockcommon%optstr) then 949 ABI_MALLOC(fockcommon%stress_ikpt,(6,nband_k)) 950 fockcommon%stress_ikpt=zero 951 end if 952 end if 953 if ((usefock_loc).and.(psps%usepaw==1)) then 954 if (fockcommon%optfor) then 955 fockcommon%forces_ikpt=zero 956 end if 957 end if 958 959 call timab(922,2,tsec) 960 961 do iblock=1,nblockbd 962 963 iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k) 964 iband_cprj=(iblock-1)*bandpp+1 965 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle 966 967 ! Select occupied bandsddk 968 occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index) 969 if( abs(maxval(occblock))>=tol8 ) then 970 call timab(923,1,tsec) 971 weight(:)=wtk(ikpt)*occblock(:) 972 973 ! gs_hamk%ffnl_k is changed in fock_getghc, so that it is necessary to restore it when stresses are to be calculated. 974 if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then 975 call gs_hamk%load_k(ffnl_k=ffnl) 976 end if 977 978 ! Load contribution from n,k 979 cwavef(:,1:npw_k*my_nspinor*blocksize)=& 980 & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg) 981 if (psps%usepaw==1.and.usecprj_local==1) then 982 call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,& 983 & mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,& 984 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 985 end if 986 987 call timab(923,2,tsec) 988 call timab(924,-1,tsec) 989 990 lambda(1:blocksize)= eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index) 991 if (mpi_enreg%paral_kgb/=1) then 992 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,mpi_enreg,blocksize,nnlout,& 993 & paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef) 994 else 995 call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,lambda,blocksize,& 996 & mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,cwavef) 997 end if 998 if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then 999 call gs_hamk%load_k(ffnl_k=ffnl_str) 1000 end if 1001 1002 ! Accumulate non-local contributions from n,k 1003 if (optfor==1) then 1004 do iblocksize=1,blocksize 1005 ibs=nnlout*(iblocksize-1) 1006 grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout(ibs+1+ishift:ibs+3*natom+ishift) 1007 end do 1008 end if 1009 if (stress_needed==1) then 1010 do iblocksize=1,blocksize 1011 ibs=nnlout*(iblocksize-1) 1012 npsstr(1:6)=npsstr(1:6) + weight(iblocksize)*enlout(ibs+1:ibs+6) 1013 end do 1014 end if 1015 1016 call timab(924,2,tsec) 1017 1018 ! Accumulate stress tensor kinetic contributions 1019 if (stress_needed==1) then 1020 call timab(925,1,tsec) 1021 do iblocksize=1,blocksize 1022 call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,my_nspinor,& 1023 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),& 1024 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0) 1025 kinstr(1)=kinstr(1)+weight(iblocksize)*ar 1026 call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,my_nspinor,& 1027 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),& 1028 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0) 1029 kinstr(2)=kinstr(2)+weight(iblocksize)*ar 1030 call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,my_nspinor,& 1031 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),& 1032 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0) 1033 kinstr(3)=kinstr(3)+weight(iblocksize)*ar 1034 call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,my_nspinor,& 1035 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),& 1036 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0) 1037 kinstr(4)=kinstr(4)+weight(iblocksize)*ar 1038 call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,my_nspinor,& 1039 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),& 1040 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0) 1041 kinstr(5)=kinstr(5)+weight(iblocksize)*ar 1042 call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,my_nspinor,& 1043 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),& 1044 & cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0) 1045 kinstr(6)=kinstr(6)+weight(iblocksize)*ar 1046 end do 1047 call timab(925,2,tsec) 1048 end if 1049 1050 ! Accumulate stress tensor and forces for the Fock part 1051 if (usefock_loc) then 1052 if(fockcommon%optstr.or.fockcommon%optfor) then 1053 call timab(926,1,tsec) 1054 if (mpi_enreg%paral_kgb==1) then 1055 msg='forsrtnps: Paral_kgb is not yet implemented for fock stresses' 1056 ABI_BUG(msg) 1057 end if 1058 ndat=mpi_enreg%bandpp 1059 if (gs_hamk%usepaw==0) cwaveprj_idat => cwaveprj 1060 ABI_MALLOC(ghc_dum,(0,0)) 1061 do iblocksize=1,blocksize 1062 fockcommon%ieigen=(iblock-1)*blocksize+iblocksize 1063 fockcommon%iband=(iblock-1)*blocksize+iblocksize 1064 if (gs_hamk%usepaw==1) then 1065 cwaveprj_idat => cwaveprj(:,(iblocksize-1)*my_nspinor+1:iblocksize*my_nspinor) 1066 end if 1067 call fock_getghc(cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),cwaveprj_idat,& 1068 & ghc_dum,gs_hamk,mpi_enreg) 1069 if (fockcommon%optstr) then 1070 fockcommon%stress(:)=fockcommon%stress(:)+weight(iblocksize)*fockcommon%stress_ikpt(:,fockcommon%ieigen) 1071 end if 1072 if (fockcommon%optfor) then 1073 fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen) 1074 end if 1075 end do 1076 ABI_FREE(ghc_dum) 1077 call timab(926,2,tsec) 1078 end if 1079 end if ! usefock_loc 1080 end if 1081 1082 end do ! End of loop on block of bands 1083 1084 call timab(927,1,tsec) 1085 1086 ! Restore the bandfft tabs 1087 if (mpi_enreg%paral_kgb==1) then 1088 call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav) 1089 end if 1090 1091 ! Increment indexes 1092 bdtot_index=bdtot_index+nband_k 1093 if (mkmem/=0) then 1094 ibg=ibg+my_nspinor*nband_cprj_k 1095 icg=icg+npw_k*my_nspinor*nband_k 1096 ikg=ikg+npw_k 1097 end if 1098 1099 if (usefock_loc) then 1100 if (fockcommon%optstr) then 1101 ABI_FREE(fockcommon%stress_ikpt) 1102 end if 1103 end if 1104 1105 if (psps%usepaw==1) then 1106 call pawcprj_free(cwaveprj) 1107 end if 1108 ABI_FREE(cwaveprj) 1109 ABI_FREE(cwavef) 1110 1111 ABI_FREE(lambda) 1112 ABI_FREE(occblock) 1113 ABI_FREE(weight) 1114 ABI_FREE(enlout) 1115 ABI_FREE(ffnl) 1116 ABI_FREE(kg_k) 1117 ABI_FREE(kpg_k) 1118 ABI_FREE(ylm_k) 1119 ABI_FREE(ylmgr_k) 1120 ABI_FREE(ph3d) 1121 if (stress_needed==1) then 1122 ABI_FREE(kstr1) 1123 ABI_FREE(kstr2) 1124 ABI_FREE(kstr3) 1125 ABI_FREE(kstr4) 1126 ABI_FREE(kstr5) 1127 ABI_FREE(kstr6) 1128 end if 1129 if ((stress_needed==1).and.(usefock_loc).and.(psps%usepaw==1))then 1130 ABI_FREE(ffnl_str) 1131 end if 1132 1133 call timab(927,2,tsec) 1134 1135 end do ! End k point loop 1136 end do ! End loop over spins 1137 1138 call timab(928,1,tsec) 1139 1140 !Stress is equal to dE/d_strain * (1/ucvol) 1141 npsstr(:)=npsstr(:)/gs_hamk%ucvol 1142 1143 !Parallel case: accumulate (n,k) contributions 1144 if (xmpi_paral==1) then 1145 ! Forces 1146 if (optfor==1) then 1147 call timab(65,1,tsec) 1148 call xmpi_sum(grnl,spaceComm,ierr) 1149 call timab(65,2,tsec) 1150 if ((usefock_loc).and.(psps%usepaw==1)) then 1151 call xmpi_sum(fockcommon%forces,spaceComm,ierr) 1152 end if 1153 end if 1154 ! Stresses 1155 if (stress_needed==1) then 1156 call timab(65,1,tsec) 1157 call xmpi_sum(kinstr,spaceComm,ierr) 1158 call xmpi_sum(npsstr,spaceComm,ierr) 1159 if (usefock_loc) then 1160 if (fockcommon%optstr) then 1161 call xmpi_sum(fockcommon%stress,spaceComm,ierr) 1162 end if 1163 end if 1164 call timab(65,2,tsec) 1165 end if 1166 end if 1167 1168 !Do final normalizations and symmetrizations of stress tensor contributions 1169 if (stress_needed==1) then 1170 renorm_factor=-(two_pi**2)/effmass_free/gs_hamk%ucvol 1171 kinstr(:)=kinstr(:)*renorm_factor 1172 if (nsym>1) then 1173 call stresssym(gs_hamk%gprimd,nsym,kinstr,symrec) 1174 call stresssym(gs_hamk%gprimd,nsym,npsstr,symrec) 1175 if (usefock_loc) then 1176 if (fockcommon%optstr) then 1177 call stresssym(gs_hamk%gprimd,nsym,fockcommon%stress,symrec) 1178 end if 1179 end if 1180 end if 1181 end if 1182 1183 !Need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted) 1184 if (psps%usepaw==1.and.usecprj_local==1) then 1185 call pawcprj_reorder(cprj,gs_hamk%atindx1) 1186 end if 1187 1188 !Deallocate temporary space 1189 call gs_hamk%free() 1190 if (usefock_loc) then 1191 fockcommon%use_ACE=use_ACE_old 1192 end if 1193 call timab(928,2,tsec) ; call timab(920,-2,tsec) 1194 1195 end subroutine forstrnps
ABINIT/m_forstr [ Modules ]
NAME
m_forstr
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, AF, AR, MB, 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_forstr 23 24 use defs_basis 25 use defs_wvltypes 26 use m_abicore 27 use m_efield 28 use m_errors 29 use m_xmpi 30 use m_fock 31 use m_cgtools 32 use m_xcdata 33 use m_dtset 34 use m_extfpmd 35 36 use defs_datatypes, only : pseudopotential_type 37 use defs_abitypes, only : MPI_type 38 use m_time, only : timab 39 use m_geometry, only : xred2xcart, metric, stresssym 40 use m_energies, only : energies_type 41 use m_pawang, only : pawang_type 42 use m_pawrad, only : pawrad_type 43 use m_pawtab, only : pawtab_type 44 use m_paw_ij, only : paw_ij_type 45 use m_pawfgrtab, only : pawfgrtab_type 46 use m_pawrhoij, only : pawrhoij_type 47 use m_pawfgr, only : pawfgr_type 48 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_reorder, pawcprj_getdim 49 use m_paw_dfpt, only : pawgrnl 50 use libxc_functionals, only : libxc_functionals_is_hybrid 51 use m_stress, only : stress 52 use m_forces, only : forces 53 use m_initylmg, only : initylmg 54 use m_xchybrid, only : xchybrid_ncpp_cc 55 use m_kg, only : mkkpg 56 use m_hamiltonian, only : init_hamiltonian, gs_hamiltonian_type !,K_H_KPRIME 57 use m_electronpositron, only : electronpositron_type, electronpositron_calctype 58 use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_type, prep_bandfft_tabs, & 59 & bandfft_kpt_savetabs, bandfft_kpt_restoretabs 60 use m_spacepar, only : meanvalue_g, hartre 61 use m_mkffnl, only : mkffnl 62 use m_mpinfo, only : proc_distrb_cycle 63 use m_nonlop, only : nonlop 64 use m_gemm_nonlop, only : make_gemm_nonlop,gemm_nonlop_use_gemm, & 65 & gemm_nonlop_ikpt_this_proc_being_treated 66 use m_fock_getghc, only : fock_getghc 67 use m_prep_kgb, only : prep_nonlop 68 use m_paw_nhat, only : pawmknhat 69 use m_rhotoxc, only : rhotoxc 70 use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll 71 use m_cgprj, only : ctocprj 72 use m_psolver, only : psolver_hartree 73 use m_wvl_psi, only : wvl_nl_gradient 74 use m_fft, only : fourdp 75 use, intrinsic :: iso_c_binding, only : c_loc,c_f_pointer 76 77 implicit none 78 79 private
ABINIT/nres2vres [ Functions ]
NAME
nres2vres
FUNCTION
Convert a density residual into a potential residual using a first order formula: V^res(r)=dV/dn.n^res(r) =V_hartree(n^res)(r) + Kxc.n^res(r)
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver | natom= number of atoms in cell | nspden=number of spin-density components | ntypat=number of atom types | typat(natom)=type (integer) for each atom gsqcut=cutoff value on G**2 for sphere inside fft box izero=if 1, unbalanced components of Vhartree(g) are set to zero kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT nhat(nfft,nspden*usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description nresid(nfft,nspden)= the input density residual n3xccc=dimension of the xccc3d array (0 or nfft). optnc=option for non-collinear magnetism (nspden=4): 1: the whole 2x2 Vres matrix is computed 2: only Vres^{11} and Vres^{22} are computed optxc=0 if LDA part of XC kernel has only to be taken into account (even for GGA) 1 if XC kernel has to be fully taken into -1 if XC kernel does not have to be taken into account pawang <type(pawang_type)>=paw angular mesh and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data rhor(nfft,nspden)=electron density in real space (used only if Kxc was not computed before) rprimd(3,3)=dimensional primitive translation vectors (bohr) usepaw= 0 for non paw calculation; =1 for paw calculation xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) xred(3,natom)=reduced dimensionless atomic coordinates === optional inputs === vxc(cplex*nfft,nspden)=XC GS potential
OUTPUT
vresid(nfft,nspden)= the output potential residual
SOURCE
1251 subroutine nres2vres(dtset,gsqcut,izero,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,& 1252 & nkxc,nresid,n3xccc,optnc,optxc,pawang,pawfgrtab,pawrhoij,pawtab,& 1253 & rhor,rprimd,usepaw,vresid,xccc3d,xred,vxc) 1254 1255 !Arguments ------------------------------------ 1256 !scalars 1257 integer,intent(in) :: izero,my_natom,n3xccc,nfft,nkxc,optnc,optxc,usepaw 1258 real(dp),intent(in) :: gsqcut 1259 type(MPI_type),intent(in) :: mpi_enreg 1260 type(dataset_type),intent(in) :: dtset 1261 type(pawang_type),intent(in) :: pawang 1262 !arrays 1263 integer,intent(in) :: ngfft(18) 1264 real(dp),intent(in) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden) 1265 real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3),xccc3d(n3xccc),xred(3,dtset%natom) 1266 real(dp),intent(inout) :: nhat(nfft,dtset%nspden*usepaw) 1267 real(dp),intent(out) :: vresid(nfft,dtset%nspden) 1268 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*usepaw) 1269 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*usepaw) 1270 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw) 1271 real(dp),intent(in) :: vxc(nfft,dtset%nspden) !FR TODO:cplex? 1272 1273 !Local variables------------------------------- 1274 !scalars 1275 integer :: cplex,ider,idir,ipert,ispden,nhatgrdim,nkxc_cur,option,me,nproc,comm,usexcnhat 1276 logical :: has_nkxc_gga,non_magnetic_xc 1277 real(dp) :: dum,energy,m_norm_min,ucvol,vxcavg 1278 character(len=500) :: message 1279 type(xcdata_type) :: xcdata 1280 !arrays 1281 integer :: nk3xc 1282 real(dp) :: dummy6(6),gmet(3,3),gprimd(3,3),qq(3),rmet(3,3) 1283 real(dp),allocatable :: dummy(:),kxc_cur(:,:),nhatgr(:,:,:) 1284 real(dp),allocatable :: nresg(:,:),rhor0(:,:),vhres(:) 1285 1286 ! ************************************************************************* 1287 1288 !Compatibility tests: 1289 has_nkxc_gga=(nkxc==7.or.nkxc==19) 1290 1291 if(optxc<-1.or.optxc>1)then 1292 write(message,'(a,i0)')' Wrong value for optxc ',optxc 1293 ABI_BUG(message) 1294 end if 1295 1296 if((optnc/=1.and.optnc/=2).or.(dtset%nspden/=4.and.optnc/=1))then 1297 write(message,'(a,i0)')' Wrong value for optnc ',optnc 1298 ABI_BUG(message) 1299 end if 1300 1301 if(dtset%icoulomb==1.and.optxc/=-1)then 1302 write(message,'(a)')' This routine is not compatible with icoulomb==1 and optxc/=-1 !' 1303 ABI_BUG(message) 1304 end if 1305 1306 if(dtset%nspden==4.and.dtset%xclevel==2.and.optxc==1.and.(.not.has_nkxc_gga))then 1307 ABI_ERROR(' Wrong values for optxc and nkxc !') 1308 end if 1309 1310 qq=zero 1311 nkxc_cur=0 1312 m_norm_min=EPSILON(0.0_dp)**2 1313 usexcnhat=0;if (usepaw==1) usexcnhat=maxval(pawtab(1:dtset%ntypat)%usexcnhat) 1314 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 1315 if (dtset%xclevel==1.or.optxc==0) nkxc_cur= 2*min(dtset%nspden,2)-1 ! LDA: nkxc=1,3 1316 if (dtset%xclevel==2.and.optxc==1)nkxc_cur=12*min(dtset%nspden,2)-5 ! GGA: nkxc=7,19 1317 ABI_MALLOC(vhres,(nfft)) 1318 1319 !Compute different geometric tensor, as well as ucvol, from rprimd 1320 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1321 1322 !Compute density residual in reciprocal space 1323 if (dtset%icoulomb==0) then 1324 ABI_MALLOC(nresg,(2,nfft)) 1325 ABI_MALLOC(dummy,(nfft)) 1326 dummy(:)=nresid(:,1) 1327 call fourdp(1,nresg,dummy,-1,mpi_enreg,nfft,1,ngfft,0) 1328 ABI_FREE(dummy) 1329 end if 1330 1331 !For GGA, has to recompute gradients of nhat 1332 nhatgrdim=0 1333 if ((nkxc==nkxc_cur.and.has_nkxc_gga).or.(optxc==-1.and.has_nkxc_gga).or.& 1334 & (optxc/=-1.and.nkxc/=nkxc_cur)) then 1335 if (usepaw==1.and.dtset%xclevel==2.and.usexcnhat>0.and.dtset%pawnhatxc>0) then 1336 nhatgrdim=1 1337 ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3)) 1338 ider=1;cplex=1;ipert=0;idir=0 1339 call pawmknhat(dum,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 1340 & nfft,ngfft,nhatgrdim,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,& 1341 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qq,rprimd,ucvol,dtset%usewvl,xred,& 1342 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1343 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 1344 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 1345 else 1346 ABI_MALLOC(nhatgr,(0,0,0)) 1347 end if 1348 else 1349 ABI_MALLOC(nhatgr,(0,0,0)) 1350 end if 1351 1352 ABI_MALLOC(dummy,(0)) 1353 !First case: Kxc has already been computed 1354 !----------------------------------------- 1355 if (nkxc==nkxc_cur.or.optxc==-1) then 1356 1357 ! Compute VH(n^res)(r) 1358 if (dtset%icoulomb == 0) then 1359 call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,& 1360 &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres) 1361 else 1362 comm=mpi_enreg%comm_cell 1363 nproc=xmpi_comm_size(comm) 1364 me=xmpi_comm_rank(comm) 1365 call psolver_hartree(energy, (/ rprimd(1,1) / dtset%ngfft(1), & 1366 & rprimd(2,2) / dtset%ngfft(2), rprimd(3,3) / dtset%ngfft(3) /), dtset%icoulomb, & 1367 & me, comm, dtset%nfft, dtset%ngfft(1:3), nproc, dtset%nscforder, dtset%nspden, & 1368 & nresid(:,1), vhres, dtset%usewvl) 1369 end if 1370 1371 ! Compute Kxc(r).n^res(r) 1372 if (optxc/=-1) then 1373 1374 ! Collinear magnetism or non-polarized 1375 if (dtset%nspden/=4) then 1376 !Note: imposing usexcnhat=1 avoid nhat to be substracted 1377 call dfpt_mkvxc(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,& 1378 & nkxc,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy) 1379 else 1380 !FR call routine for Non-collinear magnetism 1381 ABI_MALLOC(rhor0,(nfft,dtset%nspden)) 1382 rhor0(:,:)=rhor(:,:)-nresid(:,:) 1383 !Note: imposing usexcnhat=1 avoid nhat to be substracted 1384 call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,& 1385 & nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d) 1386 ABI_FREE(rhor0) 1387 end if 1388 1389 else 1390 vresid=zero 1391 end if 1392 1393 end if 1394 1395 !2nd case: Kxc has to be computed 1396 !-------------------------------- 1397 if (nkxc/=nkxc_cur.and.optxc/=-1) then 1398 1399 ! Has to use the "initial" density to compute Kxc 1400 ABI_MALLOC(rhor0,(nfft,dtset%nspden)) 1401 rhor0(:,:)=rhor(:,:)-nresid(:,:) 1402 1403 ! Compute VH(n^res) and XC kernel (Kxc) together 1404 ABI_MALLOC(kxc_cur,(nfft,nkxc_cur)) 1405 1406 option=2;if (dtset%xclevel==2.and.optxc==0) option=12 1407 1408 call hartre(1,gsqcut,dtset%icutcoul,izero,mpi_enreg,nfft,ngfft,& 1409 &dtset%nkpt,dtset%rcut,nresg,rprimd,dtset%vcutgeo,vhres) 1410 call xcdata_init(xcdata,dtset=dtset) 1411 1412 ! To be adjusted for the call to rhotoxc 1413 nk3xc=1 1414 call rhotoxc(energy,kxc_cur,mpi_enreg,nfft,ngfft,& 1415 & nhat,usepaw,nhatgr,nhatgrdim,nkxc_cur,nk3xc,non_magnetic_xc,n3xccc,option,& 1416 & rhor0,rprimd,dummy6,usexcnhat,vresid,vxcavg,xccc3d,xcdata,vhartr=vhres) !vresid=work space 1417 if (dtset%nspden/=4) then 1418 ABI_FREE(rhor0) 1419 end if 1420 1421 ! Compute Kxc(r).n^res(r) 1422 1423 if (dtset%nspden/=4) then 1424 ! Collinear magnetism or non-polarized 1425 call dfpt_mkvxc(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,& 1426 & nkxc_cur,non_magnetic_xc,dtset%nspden,0,2,qq,nresid,rprimd,1,vresid,dummy) 1427 else 1428 ! Non-collinear magnetism 1429 ABI_MALLOC(rhor0,(nfft,dtset%nspden)) 1430 rhor0(:,:)=rhor(:,:)-nresid(:,:) 1431 call dfpt_mkvxc_noncoll(1,dtset%ixc,kxc_cur,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat,usepaw,nhatgr,nhatgrdim,& 1432 & nkxc,non_magnetic_xc,dtset%nspden,0,2,2,qq,rhor0,nresid,rprimd,1,vxc,vresid,xccc3d) 1433 ABI_FREE(rhor0) 1434 end if 1435 1436 ABI_FREE(kxc_cur) 1437 end if 1438 1439 !if (nhatgrdim>0) then 1440 ABI_FREE(nhatgr) 1441 !end if 1442 1443 !Assemble potential residual: V^res(r)=VH(n^res)(r) + Kxc(r).n^res(r) 1444 !-------------------------------------------------------------------- 1445 do ispden=1,dtset%nspden/optnc 1446 vresid(:,ispden)=vresid(:,ispden)+vhres(:) 1447 end do 1448 1449 if (dtset%icoulomb==0) then 1450 ABI_FREE(nresg) 1451 end if 1452 ABI_FREE(vhres) 1453 ABI_FREE(dummy) 1454 1455 end subroutine nres2vres