TABLE OF CONTENTS
ABINIT/afterscfloop [ Functions ]
NAME
afterscfloop
FUNCTION
Perform all calculations needed after the SCF loop, independent of the call to scfcv (with or without atomic displacements), and exclusive of print or write purposes, or deallocations.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) 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> cpus= cpu time limit in seconds deltae=change in energy between the previous and present SCF cycle dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables in this dataset | mband=maximum number of bands | mgfft=maximum size of 1D FFTs (see NOTES at beginning of scfcv) | mkmem=maximum number of k points in core memory | mpw = maximum number of plane waves | natom=number of atoms in cell | nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv) | nkpt=number of k points in Brillouin zone | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | nsym=number of symmetries in space group 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_constrained_DFT)/d(xred) (hartree) grewtn(3,natom)=d(Ewald)/d(xred) (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree) gsqcut=cutoff value on G**2 for (large) sphere inside FFT box. gsqcut=(boxcut**2)*dtset%ecut/(2._dp*(Pi**2) hdr <type(hdr_type)>=the header of wf, den and pot files indsym(4,nsym,natom)=index showing transformation of atom labels under symmetry operations (computed in symatm) intgres(nspden,natom)=integrated residuals from constrained DFT. They are also Lagrange parameters, or gradients with respect to constraints. irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istep=number of the SCF iteration istep_fock_outer=number of outer SCF iteration in the double loop approach istep_mix=number of inner SCF iteration in the double loop approach kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere kxc(nfftf,nkxc)=XC kernel 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 (see NOTES at beginning of scfcv) ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc moved_atm_inside: if==1, the atoms are allowed to move. mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor n3xccc=dimension of the xccc3d array (0 or nfft). nattyp(dtset%ntypat)=number of atoms of each type nfftf= - PAW only - number of FFT grid points for the "fine" grid (see NOTES at beginning of scfcv) ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv) ngfftf(18)= - PAW only - contain all needed information about 3D FFT for the "fine" grid nhat(nfftf,nspden*psps%usepaw)= -PAW only- compensation density nkxc=dimension of kxc npwarr(nkpt)=number of planewaves in basis and on boundary for each k nvresid(nfftf,nspden)=array for the residual of the density/potential occ(mband*nkpt*nsppol)=occupancies of bands at various k points optres=0: the potential residual has been computed in scfcv 1: the density residual has been computed in scfcv paw_an(my_natom*usepaw) <type(paw_an_type)>=paw arrays given on angular mesh 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 pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data pel(3)=reduced coordinates of the electronic polarization (a. u.) pel_cg(3) = reduced coordinates of the electronic polarization (a. u.) computed in the SCF loop ph1df(2,3*(2*mgfftf+1)*natom)= - PAW only - 1-dim structure factor phases for the "fine" grid (see NOTES at beginning of scfcv) phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases pion(3)=reduced coordinates of the ionic polarization (a. u.) prtfor=1 only if forces have to be printed (0 otherwise) prtxml=1 if XML file has to be output psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f) pwind_alloc = first dimension of pwind pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations res2=density/potential residual (squared) resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins residm=maximum value from resid array (except for nbdbuf highest bands) rhog(2,nfftf)=Fourier transform of total electron density (including compensation density in PAW) rhor(nfftf,nspden)=total electron density (including compensation density in PAW) rprimd(3,3)=dimensional primitive translations in real space (bohr) stress_needed=1 if stresses are needed, 0 otherwise strscondft(6)=cDFT correction to stress strsxc(6)=xc correction to stress symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates tollist(12)=list of tolerances usecprj=1 if cprj datastructure has been allocated vhartr(nfftf)=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) vxcavg=vxc average 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 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
conv_retcode= Non-zero if convergence is not achieved. elfr(nfftf,nspden)=electron localization function grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space results_gs <type(results_gs_type)>=results (energy and its components, forces and its components, the stress tensor) of a ground-state computation (should be made a pure output quantity) taug(2,nfftf)=Fourier transform of total kinetic energy density taur(nfftf,nspden)=total kinetic energy density in real space ==== if forces are required ==== 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. gred(3,natom)=symmetrized grtn = d(etotal)/d(xred) gresid(3,natom)=forces due to the residual of the 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 tensor is required ==== 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
computed_forces=1 if forces have been computed, 0 otherwise dtefield <type(efield_type)> = variables related to Berry phase and electric field calculations (see initberry.f). In case dtset%berryopt = 4/6/7/14/16/17, the overlap matrices computed in this routine are stored in dtefield%smat in order to be used in the electric field calculation. electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation energies <type(energies_type)>=all part of total energy. | entropy(IN)=entropy due to the occupation number smearing (if metal) | e_localpsp(IN)=local psp energy (hartree) | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree) | e_chempot(IN)=energy from spatially varying chemical potential (hartree) | e_ewald(IN)=Ewald energy (hartree) | e_vdw_dftd(IN)=VdW DFT-D energy | e_hartree(IN)=Hartree part of total energy (hartree units) | e_corepsp(IN)=psp core-core energy | e_hybcomp_E0=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent den | e_kinetic(IN)=kinetic energy part of total energy. | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy. | e_xc(IN)=exchange-correlation energy (hartree) | e_xcdc(IN)=exchange-correlation double-counting energy (hartree) | e_paw(IN)=PAW spherical part energy | e_pawdc(IN)=PAW spherical part double-counting energy | e_elecfield(OUT)=the term of the energy functional that depends explicitely !!HONG | on the electric field: | enefield = -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j for fixed E/ebar | = Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j for fixed D/d etotal=total energy, might be correct by improved polarization computation forold(3,natom)=old forces xred(3,natom)=reduced dimensionless atomic coordinates ===== if dtset%densfor_pred==3 .and. moved_atm_inside==1 ===== ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse grid) ph1df(2,3*(2*mgfftf+1)*natom)=1-dim structure factor phases (fine PAW grid) wvl <type(wvl_data)>=all wavelets data.
NOTES
SOURCE
276 subroutine afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,& 277 & deltae,diffor,dtefield,dtfil,dtset,eigen,electronpositron,elfr,& 278 & energies,etotal,favg,fcart,fock,forold,grchempottn,grcondft,& 279 & gred,gresid,grewtn,grhf,grhor,grvdw,& 280 & grxc,gsqcut,hdr,extfpmd,indsym,intgres,irrzon,istep,istep_fock_outer,istep_mix,& 281 & kg,kxc,lrhor,maxfor,mcg,mcprj,mgfftf,& 282 & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,& 283 & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,& 284 & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,prtxml,& 285 & psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,& 286 & rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,taug,& 287 & taur,tollist,usecprj,vhartr,vpsp,vtrial,vxc,vxctau,vxcavg,wvl,& 288 & xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero,conv_retcode) 289 290 !Arguments ------------------------------------ 291 !scalars 292 integer,intent(in) :: istep,istep_fock_outer,istep_mix 293 integer,intent(in) :: mcg,mcprj,mgfftf,moved_atm_inside,my_natom,n3xccc,nfftf,ngrvdw,nkxc 294 integer,intent(in) :: optres,prtfor,prtxml,pwind_alloc,stress_needed,usecprj 295 integer,intent(inout) :: computed_forces 296 real(dp),intent(in) :: cpus,deltae,gsqcut,res2,residm 297 real(dp),intent(in) :: qvpotzero 298 real(dp),intent(inout) :: diffor,etotal,maxfor,vxcavg 299 type(MPI_type),intent(inout) :: mpi_enreg 300 type(datafiles_type),intent(in) :: dtfil 301 type(dataset_type),intent(inout) :: dtset 302 type(efield_type),intent(inout) :: dtefield 303 type(electronpositron_type),pointer :: electronpositron 304 type(energies_type),intent(inout) :: energies 305 type(hdr_type),intent(inout) :: hdr 306 type(extfpmd_type),pointer,intent(inout) :: extfpmd 307 type(pawang_type),intent(in) :: pawang 308 type(pawfgr_type),intent(in) :: pawfgr 309 type(pseudopotential_type),intent(in) :: psps 310 type(results_gs_type),intent(inout) :: results_gs 311 type(wvl_data),intent(inout) :: wvl 312 type(fock_type),pointer, intent(inout) :: fock 313 !arrays 314 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom) 315 integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom) 316 integer,intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 317 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%ntypat) 318 integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt) 319 integer,intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym) 320 integer,intent(out) :: conv_retcode 321 real(dp),intent(in) :: grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw) 322 real(dp),intent(in) :: grcondft(:,:) ! (3,natom) if constrainedDFT otherwise (3,0) 323 real(dp),intent(in) :: intgres(:,:) ! (nspden,natom) if constrainedDFT otherwise (nspden,0) 324 real(dp),intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 325 real(dp),intent(in) :: pwnsfac(2,pwind_alloc) 326 real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 327 real(dp),intent(in) :: strscondft(6) 328 real(dp),intent(in) :: tollist(12),vpsp(nfftf) 329 real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden) 330 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 331 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 332 real(dp),intent(inout) :: cg(2,mcg) 333 real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 334 real(dp),intent(inout) :: forold(3,dtset%natom) 335 real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw) 336 real(dp),intent(inout) :: nvresid(nfftf,dtset%nspden),pel(3) 337 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),pel_cg(3) 338 real(dp),intent(inout) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 339 real(dp),intent(inout) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),pion(3) 340 real(dp),intent(inout) :: rprimd(3,3) 341 real(dp),intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),strsxc(6) 342 real(dp),intent(inout) :: vhartr(nfftf),vxc(nfftf,dtset%nspden),vxctau(nfftf,dtset%nspden,4*dtset%usekden) 343 real(dp),intent(inout) :: xccc3d(n3xccc),xcctau3d(n3xccc*dtset%usekden),xred(3,dtset%natom) 344 real(dp),intent(inout) :: favg(3),fcart(3,dtset%natom),gred(3,dtset%natom) 345 real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom) 346 real(dp),intent(inout) :: grxc(3,dtset%natom),kxc(nfftf,nkxc),strten(6) 347 real(dp),intent(inout) :: synlgr(3,dtset%natom) 348 real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taug(:,:),taur(:,:) 349 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj) 350 type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw) 351 type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw) 352 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 353 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 354 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw) 355 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 356 357 !Local variables------------------------------- 358 !scalars 359 integer,parameter :: response=0 360 integer :: bantot,bufsz,calc_pol_ddk,choice,cplex,ierr,ifft,igrad,ishift,ispden 361 integer :: mcg1_3,nfftotf,ngrad,optcut,optfor,optgr0,optgr1,optgr2,optrad,quit,shft 362 integer :: spaceComm_fft,tim_mkrho 363 logical :: save_cg1_3,test_gylmgr,test_nfgd,test_rfgd 364 logical :: wvlbigdft=.false. 365 real(dp) :: c_fermi,dtaur,dtaurzero,ucvol 366 character(len=500) :: message 367 type(paw_dmft_type) :: paw_dmft 368 #if defined HAVE_BIGDFT 369 integer :: ia,ii,mband_cprj 370 logical :: do_last_ortho,compute_wvl_tail=.false. 371 real(dp) :: dum,eexctx,eh,ekin,eloc,enl,eproj,esicdc,evxc,exc,ucvol_local 372 #endif 373 !arrays 374 real(dp) :: gmet(3,3),gprimd(3,3),pelev(3),ptot(3),red_ptot(3),rmet(3,3),tsec(2) 375 real(dp) :: dmatdum(0,0,0,0) 376 real(dp),allocatable :: cg1_3(:,:,:),mpibuf(:,:),qphon(:),rhonow(:,:,:),sqnormgrhor(:,:) 377 real(dp),allocatable :: tauwfg(:,:),tauwfr(:,:),vtrial_local(:,:) 378 #if defined HAVE_BIGDFT 379 integer,allocatable :: dimcprj_srt(:) 380 real(dp),allocatable :: hpsi_tmp(:),xcart(:,:) 381 #endif 382 383 ! ************************************************************************* 384 385 DBG_ENTER("COLL") 386 387 call timab(250,1,tsec) 388 call timab(251,1,tsec) 389 390 !Compute different geometric tensor, as well as ucvol, from rprimd 391 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 392 nfftotf=product(ngfftf(1:3)) 393 394 !MPI FFT communicator 395 spaceComm_fft=mpi_enreg%comm_fft 396 397 !Recompute structure factor phases if atomic positions have changed 398 if (moved_atm_inside==1) then 399 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 400 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 401 else 402 ph1d(:,:)=ph1df(:,:) 403 end if 404 end if 405 406 !---------------------------------------------------------------------- 407 !Wavelet case: transform psi to KS orbitals 408 !---------------------------------------------------------------------- 409 if (dtset%usewvl == 1) then 410 411 ! wvlbigdft indicates that the BigDFT workflow will be followed 412 wvlbigdft=(dtset%wvl_bigdft_comp==1) 413 414 #if defined HAVE_BIGDFT 415 ! Transform to KS orbitals 416 417 ! Need xcart 418 ABI_MALLOC(xcart,(3, dtset%natom)) 419 call xred2xcart(dtset%natom, rprimd, xcart, xred) 420 ucvol_local=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp) 421 422 ! do_last_ortho in case of direct minimization, since 423 ! We never diagonalized the hamiltonian and the eigenvalues are unknown. 424 if ( wvlbigdft) do_last_ortho=(dtset%iscf==0) 425 if (.not.wvlbigdft) do_last_ortho=(.false.) 426 if (do_last_ortho) then 427 call total_energies(wvl%e%energs, istep, mpi_enreg%me_wvl) 428 call write_energies(istep,0,wvl%e%energs,zero,zero,"FINAL") 429 if(.not.wvlbigdft) then 430 ! If ISCF>10, we exit scfcv at a place where bigdft objects 431 ! do not contain the KS potential. Hence, we copy vtrial to wvl%den 432 if(dtset%iscf>=10) call wvl_vtrial_abi2big(1,vtrial,wvl%den) 433 ! hpsi is lost in hpsitopsi, so we recalculate it (needed for last_orthon). 434 call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,& 435 & istep,1,-1,mpi_enreg%me_wvl,dtset%natom,& 436 & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,& 437 & dum,.false.,evxc,wvl,wvlbigdft,xcart,strsxc) 438 if (dtset%iscf==0) then 439 energies%e_kinetic=ekin ; energies%e_hartree=eh 440 energies%e_xc=exc ; energies%e_localpsp=eloc 441 energies%e_nlpsp_vfock=enl ; energies%e_exactX=eexctx 442 energies%e_sicdc=esicdc ; energies%e_xcdc=evxc 443 energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp & 444 & + energies%e_xcdc + two*energies%e_hartree +energies%e_nlpsp_vfock 445 end if 446 end if 447 call last_orthon(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,istep,wvl%wfs%ks,wvl%e%energs%evsum,.true.) 448 if (mpi_enreg%nproc_wvl == 1) nullify(wvl%wfs%ks%psit) 449 call eigensystem_info(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,0.d0,& 450 wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,& 451 wvl%wfs%ks%orbs,wvl%wfs%ks%psi) 452 ! Copy eigenvalues from BigDFT object to "eigen" 453 call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs) 454 ! Copy occupations from BigDFT objects to ABINIT 455 call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs) 456 end if 457 458 ! Tail corrections, pending for wvlbigdft==.false. 459 ! TODO put it at the end of gstate. 460 ! WVL - maybe compute the tail corrections to energy 461 compute_wvl_tail=(dtset%tl_radius>tol12.and.wvlbigdft) 462 if (compute_wvl_tail) then 463 ! Use the tails to improve energy precision. 464 call wvl_tail_corrections(dtset, energies, etotal, mpi_enreg, psps, wvl, xcart) 465 end if 466 467 ! Clean KSwfn parts only needed in the SCF loop. 468 call kswfn_free_scf_data(wvl%wfs%ks, (mpi_enreg%nproc_wvl > 1)) 469 ! Clean denspot parts only needed in the SCF loop. 470 call denspot_free_history(wvl%den%denspot) 471 472 ! If WF have been modified, change the density according to the KS projection. 473 if ( do_last_ortho ) then 474 475 ! Density from new orthogonalized WFs 476 call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den) 477 478 ! PAW: has to update cprj, rhoij and compensation charge density 479 if (psps%usepaw==1) then 480 ! 1-Compute cprj 481 ABI_MALLOC(hpsi_tmp,(size(wvl%wfs%ks%hpsi))) 482 call applyprojectorsonthefly(mpi_enreg%me_wvl,wvl%wfs%ks%orbs,wvl%descr%atoms,wvl%descr%Glr,& 483 & xcart,wvl%descr%h(1),wvl%descr%h(2),wvl%descr%h(3),wvl%wfs%ks%lzd%Glr%wfd,& 484 & wvl%projectors%nlpsp,wvl%wfs%ks%psi,hpsi_tmp,eproj,& 485 & proj_G=wvl%projectors%G,paw=wvl%descr%paw) 486 ABI_FREE(hpsi_tmp) 487 do ii=1,mcprj 488 do ia=1,dtset%natom 489 !Note that cprj should be allocated (i.e. usepcrj=1 imposed in scfcv) 490 cprj(ia,ii)%cp(:,:)= wvl%descr%paw%cprj(ia,ii)%cp(:,:) 491 end do 492 end do 493 ! 2-Compute rhoij 494 ABI_MALLOC(dimcprj_srt,(dtset%natom)) 495 call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 496 mband_cprj=mcprj/(dtset%nspinor*dtset%mkmem*dtset%nsppol) 497 paw_dmft%use_sc_dmft=0 ; paw_dmft%use_dmft=0 ! dmft not used here 498 call pawmkrhoij(atindx,atindx1,cprj,dimcprj_srt,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,& 499 & mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspden,dtset%nspinor,& 500 & dtset%nsppol,occ,mpi_enreg%paral_kgb,paw_dmft,pawrhoij,dtfil%unpaw,dtset%usewvl,dtset%wtk) 501 ABI_FREE(dimcprj_srt) 502 ! 3-Symetrize rhoij, compute nhat and add it to rhor 503 call pawmkrho(1,dum,1,gprimd,0,indsym,0,mpi_enreg,my_natom,dtset%natom,dtset%nspden,dtset%nsym,& 504 & dtset%ntypat,mpi_enreg%paral_kgb,pawang,pawfgr,pawfgrtab,dtset%pawprtvol,pawrhoij,pawrhoij,& 505 & pawtab,(/zero,zero,zero/),rhog,rhor,rhor,rprimd,dtset%symafm,symrec,dtset%typat,ucvol_local,& 506 & dtset%usewvl,xred,pawnhat=nhat) 507 call wvl_rho_abi2big(1,rhor,wvl%den) 508 end if 509 end if 510 511 ABI_FREE(xcart) 512 513 #else 514 BIGDFT_NOTENABLED_ERROR() 515 #endif 516 end if 517 518 call timab(251,2,tsec) 519 call timab(252,1,tsec) 520 521 !---------------------------------------------------------------------- 522 !Polarization Calculation, but not orbital magnetism 523 !---------------------------------------------------------------------- 524 525 if(dtset%berryopt/=0 .AND. dtset%orbmag == 0)then 526 call elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,energies%e_elecfield,gprimd,hdr,& 527 & kg,dtset%mband,mcg,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,dtset%nkpt,& 528 & npwarr,dtset%nsppol,psps%ntypat,pawrhoij,pawtab,pel,pel_cg,pelev,pion,& 529 & psps,pwind,pwind_alloc,pwnsfac,rprimd,ucvol,usecprj,xred) 530 end if 531 532 call timab(252,2,tsec) 533 call timab(253,1,tsec) 534 535 !---------------------------------------------------------------------- 536 !Orbital magnetization calculation using PEAD DDK wavefunctions 537 !---------------------------------------------------------------------- 538 539 if (dtset%berryopt == -2 .AND. dtset%orbmag /= 0) then 540 save_cg1_3 = .TRUE. 541 mcg1_3 = mcg 542 ABI_MALLOC(cg1_3,(2,mcg1_3,3)) 543 544 calc_pol_ddk = 2 545 call berryphase_new(atindx1,cg,cg1_3,cprj,dtefield,dtfil,dtset,psps,& 546 & gprimd,hdr,psps%indlmn,kg,psps%lmnmax,dtset%mband,mcg,mcg1_3,mcprj,& 547 & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,npwarr,dtset%nsppol,psps%ntypat,& 548 & dtset%nkpt,calc_pol_ddk,pawrhoij,pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC 549 & pwind_alloc,pwnsfac,rprimd,save_cg1_3,dtset%typat,ucvol,ab_out,& 550 & usecprj,psps%usepaw,xred,psps%ziontypat) 551 552 if ( .NOT. ALLOCATED(vtrial_local)) then 553 ABI_MALLOC(vtrial_local,(nfftf,dtset%nspden)) 554 end if 555 vtrial_local = vtrial 556 call orbmag(cg,cg1_3,cprj,dtset,eigen,gsqcut,kg,mcg,mcg1_3,mcprj,dtset%mkmem,& 557 & mpi_enreg,dtset%mpw,nfftf,ngfftf,npwarr,occ,paw_an,paw_ij,pawang,pawfgr,& 558 & pawrad,pawtab,psps,rprimd,vtrial_local,vxctau,xred,ylm,ylmgr) 559 560 ABI_FREE(vtrial_local) 561 ABI_FREE(cg1_3) 562 563 end if 564 565 !---------------------------------------------------------------------- 566 !Gradient and Laplacian of the Density Calculation 567 !---------------------------------------------------------------------- 568 569 !We use routine xcden which get gradient of rhor (grhor), and eventually laplacian of rhor (lrhor). 570 if(dtset%prtgden/=0 .or. dtset%prtlden/=0)then 571 572 ! Compute gradient of the electron density 573 ngrad=2 574 cplex=1 575 ishift=0 576 ABI_MALLOC(rhonow,(nfftf,dtset%nspden,ngrad*ngrad)) 577 if(dtset%prtlden/=0)then 578 nullify(lrhor) 579 ABI_MALLOC(lrhor,(nfftf,dtset%nspden)) 580 end if 581 write(message,'(a,a)') ch10, " Compute gradient of the electron density" 582 call wrtout(ab_out,message) 583 if(dtset%prtlden/=0) then 584 write(message,'(a)') " and also Compute Laplacian of the electron density" 585 call wrtout(ab_out,message) 586 end if 587 write(message,'(a)') "--------------------------------------------------------------------------------" 588 call wrtout(ab_out,message) 589 590 ABI_MALLOC(qphon,(3)) 591 qphon(:)=zero 592 if(dtset%prtlden/=0) then 593 call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,qphon,rhor,rhonow,lrhonow=lrhor) 594 else 595 call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,qphon,rhor,rhonow) 596 end if 597 ABI_FREE(qphon) 598 599 ! Copy gradient which has been output in rhonow to grhor (and free rhonow) 600 nullify(grhor) 601 ABI_MALLOC(grhor,(nfftf,dtset%nspden,3)) 602 do ispden=1,dtset%nspden 603 do ifft=1,nfftf 604 grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4) 605 end do 606 end do 607 ABI_FREE(rhonow) 608 609 if(dtset%prtgden/=0) then 610 ! Print result for grhor 611 write(message,'(a,a)') ch10, " Result for gradient of the electron density for each direction (1,2,3):" 612 call wrtout(ab_out,message) 613 write(message,'(a,a,a,a)') ch10," 1rst direction:",ch10,& 614 & "--------------------------------------------------------------------------------" 615 call wrtout(ab_out,message) 616 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,1),optrhor=2,ucvol=ucvol) 617 write(message,'(a)') "--------------------------------------------------------------------------------" 618 call wrtout(ab_out,message) 619 write(message,'(a,a,a,a)') ch10," 2nd direction:",ch10,& 620 & "--------------------------------------------------------------------------------" 621 call wrtout(ab_out,message) 622 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,2),optrhor=2,ucvol=ucvol) 623 write(message,'(a)') "--------------------------------------------------------------------------------" 624 call wrtout(ab_out,message) 625 write(message,'(a,a,a,a)') ch10," 3rd direction:",ch10,& 626 & "--------------------------------------------------------------------------------" 627 call wrtout(ab_out,message) 628 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,3),optrhor=2,ucvol=ucvol) 629 write(message,'(a)') "--------------------------------------------------------------------------------" 630 call wrtout(ab_out,message) 631 end if 632 633 if(dtset%prtlden/=0) then 634 ! Print result for lrhor 635 write(message,'(a,a)') ch10, " Result for Laplacian of the electron density :" 636 call wrtout(ab_out,message) 637 write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------" 638 call wrtout(ab_out,message) 639 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,lrhor,optrhor=3,ucvol=ucvol) 640 write(message,'(a)') "--------------------------------------------------------------------------------" 641 call wrtout(ab_out,message) 642 end if 643 644 write(message,'(a)') "--------------------------------------------------------------------------------" 645 call wrtout(ab_out,message) 646 end if 647 648 !---------------------------------------------------------------------- 649 !Kinetic Energy Density Calculation 650 !---------------------------------------------------------------------- 651 652 call timab(253,2,tsec) 653 call timab(254,1,tsec) 654 655 !We use routine mkrho with option=1 to compute kinetic energy density taur (and taug) 656 if(dtset%usekden==0 .and. dtset%prtelf/=0)then 657 ! tauX are reused in outscfcv for output 658 ! should be deallocated there 659 nullify(taug,taur) 660 ABI_MALLOC(taug,(2,nfftf)) 661 ABI_MALLOC(taur,(nfftf,dtset%nspden)) 662 tim_mkrho=5 663 if(dtset%prtelf/=0) then 664 write(message,'(a,a)') ch10, " Compute ELF" 665 call wrtout(ab_out,message) 666 write(message,'(a)') "--------------------------------------------------------------------------------" 667 call wrtout(ab_out,message) 668 end if 669 write(message,'(a,a)') ch10, " Compute kinetic energy density" 670 call wrtout(ab_out,message) 671 paw_dmft%use_sc_dmft=0 ! dmft not used here 672 paw_dmft%use_dmft=0 ! dmft not used here 673 if (psps%usepaw==0) then 674 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,& 675 & npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,& 676 & extfpmd=extfpmd,option=1) 677 else 678 ABI_MALLOC(tauwfg,(2,dtset%nfft)) 679 ABI_MALLOC(tauwfr,(dtset%nfft,dtset%nspden)) 680 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,& 681 & npwarr,occ,paw_dmft,phnons,tauwfg,tauwfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,& 682 & extfpmd=extfpmd,option=1) 683 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,tauwfg,taug,tauwfr,taur) 684 ABI_FREE(tauwfg) 685 ABI_FREE(tauwfr) 686 end if 687 ABI_FREE(taug) 688 end if 689 !Print result 690 if(dtset%prtkden/=0) then 691 write(message,'(a,a)') ch10, "Result for kinetic energy density :" 692 call wrtout(ab_out,message) 693 write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------" 694 call wrtout(ab_out,message) 695 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,taur,optrhor=1,ucvol=ucvol) 696 write(message,'(a)') "--------------------------------------------------------------------------------" 697 call wrtout(ab_out,message) 698 end if 699 700 !---------------------------------------------------------------------- 701 !Electron Localization Function (ELF) Calculation 702 !---------------------------------------------------------------------- 703 704 call timab(254,2,tsec) 705 call timab(255,1,tsec) 706 707 !We use routine xcden to compute gradient of electron density (grhor), 708 !NOTE: If GGA is used, gradient of electron density is already computed 709 !and it is stored in exchange correlation kernel kxc(:,5:7) (nspden=1) or kxc(:,14:19) (nspden=2). 710 !In order to save memory and do not have the same quantity twice 711 !in memory we should use kxc. 712 !But unfortunately only spin up ans spin down component are stored in kxc. 713 !So far we thus use grhor which contains all component (nspden=4) 714 !just like rhonow in xcden instead of kxc. 715 716 if((dtset%prtelf/=0))then 717 if(dtset%nspden<=2) then 718 719 ngrad=2 720 cplex=1 721 if((cplex*dtset%nfft)/=nfftf)then 722 write(message, '(a,a,a,a)' ) ch10,& 723 & ' afterscfloop: ERROR -', ch10, & 724 & ' The density is complex, ELF analysis cannot be performed.' 725 call wrtout(std_out,message) 726 ! ABI_ERROR(message) 727 end if 728 729 if((dtset%prtgden==0) .and. (dtset%prtlden==0)) then 730 ! Compute gradient of the electron density 731 ishift=0 732 ABI_MALLOC(rhonow,(nfftf,dtset%nspden,ngrad*ngrad)) 733 write(message,'(a,a)') ch10, " Compute gradient of the electron density" 734 call wrtout(ab_out,message) 735 ABI_MALLOC(qphon,(3)) 736 qphon(:)=zero 737 call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,qphon,rhor,rhonow) 738 ABI_FREE(qphon) 739 ! Copy gradient which has been output in rhonow to grhor (and free rhonow) 740 ABI_MALLOC(grhor,(nfftf,dtset%nspden,3)) 741 do ispden=1,dtset%nspden 742 do ifft=1,nfftf 743 grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4) 744 end do 745 end do 746 ABI_FREE(rhonow) 747 end if 748 ! Compute square norm of gradient of the electron density (|grhor|**2) 749 if(dtset%nspden==1)then 750 ABI_MALLOC(sqnormgrhor,(nfftf,dtset%nspden)) 751 do ifft=1,nfftf 752 sqnormgrhor(ifft,1) = zero 753 end do 754 elseif(dtset%nspden==2)then 755 ABI_MALLOC(sqnormgrhor,(nfftf,dtset%nspden+1)) 756 ! because we not only want (total and up quantities, but also down) 757 ! Indeed after having token the square norm we can not recover the 758 ! down quantity by substracting total and up quantities (as we do for usual densities) 759 do ispden=1,dtset%nspden+1 760 do ifft=1,nfftf 761 sqnormgrhor(ifft,ispden) = zero 762 end do 763 end do 764 end if 765 766 do igrad=1,3 767 do ispden=1,dtset%nspden !total (and up) 768 do ifft=1,nfftf 769 sqnormgrhor(ifft,ispden) = sqnormgrhor(ifft,ispden) + grhor(ifft,ispden,igrad)**2 770 end do 771 end do 772 if(dtset%nspden==2)then 773 do ifft=1,nfftf !down 774 sqnormgrhor(ifft,3) = sqnormgrhor(ifft,3) + (grhor(ifft,1,igrad)-grhor(ifft,2,igrad))**2 775 end do 776 end if 777 end do 778 779 ! Compute electron localization function (ELF) (here it is elfr) 780 781 nullify(elfr) 782 if(dtset%nspden==1)then 783 ABI_MALLOC(elfr,(nfftf,dtset%nspden)) 784 elseif(dtset%nspden==2)then 785 ABI_MALLOC(elfr,(nfftf,dtset%nspden+1)) 786 ! 1rst is total elf, 2nd is spin-up elf, and 3rd is spin-down elf. (elf_tot /= elf_up + elf_down) 787 end if 788 c_fermi = 3.0d0/10.0d0*((3.0d0*pi**2)**(2.0d0/3.0d0)) 789 790 ! First compute total elf 791 ispden=1 792 do ifft=1,nfftf 793 dtaurzero = c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0) 794 dtaur = taur(ifft,ispden) 795 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden)) 796 ! Ensure that dtaur is always positive or zero, as it should be. 797 if(dtaur<0.0d0)dtaur=0.0d0 798 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 799 if(dtaurzero<(1.0d-20*dtaur)) then 800 elfr(ifft,ispden) = 0.0d0 801 else 802 elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 803 ! For atomic shell studies we could also add the condition that when (dtaur/dtaurzero) 804 ! is very close to zero we actually set it exactly to zero (or --> elfr = 1.0d0 805 ! which is the upper limit of elfr.) 806 end if 807 end do 808 809 ! If spin-dependent densities are avalaible, compute spin-dependent elf 810 ! and offer the possibility to compute total elf in an alternative approach 811 ! (see doc/theory/ELF) 812 if(dtset%nspden==2)then 813 814 ! alternative approach to the total elf 815 if(dtset%prtelf==2)then 816 ispden=1 817 do ifft=1,nfftf 818 dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*( rhor(ifft,ispden+1)**(5.0d0/3.0d0) + & 819 & (rhor(ifft,ispden) - rhor(ifft,ispden+1))**(5.0d0/3.0d0) ) 820 dtaur = taur(ifft,ispden) 821 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+1)/rhor(ifft,ispden+1)) & 822 & - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1))) 823 if(dtaur<0.0d0)dtaur=0.0d0 824 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 825 if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then 826 elfr(ifft,ispden) = 0.0d0 827 else 828 elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 829 end if 830 end do 831 end if 832 833 ! elf_up 834 ispden=2 835 do ifft=1,nfftf 836 dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0) 837 dtaur = taur(ifft,ispden) 838 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden)) 839 if(dtaur<0.0d0)dtaur=0.0d0 840 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 841 if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then 842 elfr(ifft,ispden) = 0.0d0 843 else 844 elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 845 end if 846 end do 847 848 ! elf_down 849 ispden=1 850 do ifft=1,nfftf 851 dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*(rhor(ifft,ispden)-rhor(ifft,ispden+1))**(5.0d0/3.0d0) 852 dtaur = taur(ifft,ispden)-taur(ifft,ispden+1) 853 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1))) 854 if(dtaur<0.0d0)dtaur=0.0d0 855 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 856 if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then 857 elfr(ifft,ispden+2) = 0.0d0 858 else 859 elfr(ifft,ispden+2) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 860 end if 861 end do 862 863 end if !endif dtset%nspden==2 864 865 ! Print result for elfr 866 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,elfr,optrhor=4,ucvol=ucvol) 867 868 ABI_FREE(grhor) 869 ABI_FREE(sqnormgrhor) 870 871 else 872 message ='ELF is not yet implemented for non collinear spin cases.' 873 ABI_WARNING(message) 874 875 ABI_MALLOC(elfr,(nfftf,dtset%nspden)) 876 do ispden=1,dtset%nspden 877 do ifft=1,nfftf 878 elfr(ifft,ispden) = -2.0d0 879 end do 880 end do 881 ! even if elf is not computed we want to finish the abinit run. 882 ! To ensure that users won't use the _ELF file which will be produced 883 ! we set elf to -2.0 (a meaningless value) 884 885 end if ! endif dtset%nspden<=2 886 887 write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------" 888 call wrtout(ab_out,message) 889 write(message,'(a)') " End of ELF section" 890 call wrtout(ab_out,message) 891 892 if (dtset%usekden==0) then 893 ABI_FREE(taur) 894 end if 895 896 end if !endif prtelf/=0 897 898 !###################################################################### 899 !Compute forces (if they were not computed during the elec. iterations) 900 !and stresses (if requested by user) 901 !---------------------------------------------------------------------- 902 903 call timab(255,2,tsec) 904 call timab(256,1,tsec) 905 906 optfor=0 907 908 if (computed_forces==0.and.dtset%optforces>0.and.dtset%iscf>=0) then 909 if (dtset%nstep>0.or.dtfil%ireadwf==1) optfor=1 910 end if 911 912 if (optfor>0.or.stress_needed>0) then 913 914 ! PAW: eventually, compute g_l(r).Y_lm(r) gradients (if not already done) 915 if (psps%usepaw==1) then 916 test_nfgd =any(pawfgrtab(:)%nfgd==0) 917 test_rfgd =any(pawfgrtab(:)%rfgd_allocated==0) 918 test_gylmgr=any(pawfgrtab(:)%gylmgr_allocated==0) 919 if (test_nfgd.or.& 920 & (test_gylmgr.and.dtset%pawstgylm==1).or.& 921 & (test_rfgd.and.stress_needed==1.and.dtset%pawstgylm==1).or.& 922 (test_rfgd.and.dtset%pawstgylm==0)) then 923 optcut=0;optgr0=0;optgr1=dtset%pawstgylm;optgr2=0 924 optrad=1-dtset%pawstgylm;if (stress_needed==1) optrad=1 925 if (dtset%usewvl==0) then 926 call nhatgrid(atindx1,gmet,my_natom,dtset%natom,nattyp,ngfftf,dtset%ntypat,& 927 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,& 928 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 929 & comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft) 930 else 931 shft=0 932 #if defined HAVE_BIGDFT 933 shft=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,4) 934 call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,& 935 & wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,& 936 & nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,& 937 & wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,& 938 & wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,& 939 & pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shft,xred) 940 #endif 941 end if 942 end if 943 end if 944 945 call forstr(atindx1,cg,cprj,diffor,dtefield,dtset,& 946 & eigen,electronpositron,energies,favg,fcart,fock,& 947 & forold,gred,grchempottn,grcondft,gresid,grewtn,& 948 & grhf,grvdw,grxc,gsqcut,extfpmd,indsym,& 949 & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,& 950 & n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,& 951 & npwarr,dtset%ntypat,nvresid,occ,optfor,optres,& 952 & paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,& 953 & psps,rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,& 954 & ucvol,usecprj,vhartr,vpsp,vxc,vxctau,wvl,xccc3d,xcctau3d,xred,ylm,ylmgr,qvpotzero) 955 end if 956 957 ! Init values with MAGIC_UNDEF if not computed. 958 if (optfor==1) computed_forces=1 959 if (optfor==1) diffor = MAGIC_UNDEF 960 if (stress_needed==0) strten = MAGIC_UNDEF 961 if (computed_forces==0) fcart = MAGIC_UNDEF 962 if (dtset%prtstm/=0) strten(:)=zero 963 964 call timab(256,2,tsec) 965 call timab(257,1,tsec) 966 967 !If SCF convergence was not reached (for dtset%nstep>0), 968 !print a warning to the output file (non-dummy arguments: dtset%nstep, 969 !residm, diffor - infos from tollist have been saved inside ) 970 choice=3 971 call scprqt(choice,cpus,deltae,diffor,dtset,& 972 & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,& 973 & dtfil%fnameabo_app_eig,dtfil%filnam_ds(1),& 974 & 1,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,maxfor,& 975 & moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,& 976 & dtset%nstep,occ,optres,prtfor,prtxml,quit,& 977 & res2,resid,residm,response,tollist,psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,& 978 & electronpositron=electronpositron, fock=fock) 979 980 !output POSCAR and FORCES files, VASP style, for PHON code and friends. 981 if (dtset%prtposcar == 1) then 982 call prtposcar(fcart, dtfil%filnam_ds(4), dtset%natom, dtset%ntypat, rprimd, dtset%typat, ucvol, xred, dtset%znucl) 983 end if ! prtposcar 984 985 if(allocated(qphon)) then 986 ABI_FREE(qphon) 987 end if 988 989 !get current operator on wavefunctions 990 if (dtset%prtspcur == 1) then 991 call spin_current(cg,dtfil,dtset,gprimd,hdr,kg,mcg,mpi_enreg,psps) 992 end if 993 994 !Electron-positron stuff: if last calculation was a positron minimization, 995 !exchange electron and positron data in order to 996 !get electronic quantities in global variables 997 if (dtset%positron/=0) then 998 electronpositron%scf_converged=.false. 999 if (dtset%positron<0.and.electronpositron_calctype(electronpositron)==1) then 1000 call exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,gred,mcg,mcprj,& 1001 & mpi_enreg,my_natom,nfftf,ngfftf,nhat,npwarr,occ,paw_an,pawrhoij,rhog,rhor,strten,usecprj,vhartr) 1002 end if 1003 end if 1004 1005 !If PAW+U and density mixing, has to update nocc_mmp 1006 if (psps%usepaw==1.and.dtset%usepawu/=0.and.(dtset%iscf>0.or.dtset%iscf==-3)) then 1007 call setnoccmmp(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,& 1008 & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,& 1009 & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,& 1010 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1011 end if 1012 1013 !Update the content of the header (evolving variables) 1014 bantot=hdr%bantot 1015 if (dtset%positron==0) then 1016 call hdr%update(bantot,etotal,energies%e_fermie,energies%e_fermih,residm,rprimd,occ,& 1017 pawrhoij,xred,dtset%amu_orig(:,1),& 1018 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1019 else 1020 call hdr%update(bantot,electronpositron%e0,energies%e_fermie,energies%e_fermih,residm,rprimd,occ,& 1021 pawrhoij,xred,dtset%amu_orig(:,1),& 1022 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1023 end if 1024 1025 #ifdef HAVE_LOTF 1026 if(dtset%ionmov==23 .and. mpi_enreg%nproc_band>1) then 1027 bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom 1028 ABI_MALLOC(mpibuf,(3,bufsz)) 1029 mpibuf(:,1:dtset%natom)=gred(:,1:dtset%natom) 1030 mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom) 1031 if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom) 1032 mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/)) 1033 call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr) 1034 gred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_band 1035 fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_band 1036 if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_band 1037 strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_band 1038 ABI_FREE(mpibuf) 1039 end if 1040 #endif 1041 1042 !In case of FFT parallelisation, has to synchronize positions and forces 1043 !to avoid numerical noise 1044 if (mpi_enreg%nproc_fft>1) then 1045 bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom 1046 ABI_MALLOC(mpibuf,(3,bufsz)) 1047 mpibuf(:,1:dtset%natom)=gred(:,1:dtset%natom) 1048 mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom) 1049 if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom) 1050 mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/)) 1051 call xmpi_sum(mpibuf,mpi_enreg%comm_fft,ierr) 1052 gred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_fft 1053 fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_fft 1054 if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_fft 1055 strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_fft 1056 ABI_FREE(mpibuf) 1057 end if 1058 1059 !results_gs%energies = energies 1060 call energies_copy(energies,results_gs%energies) 1061 results_gs%etotal =etotal 1062 results_gs%deltae =deltae 1063 results_gs%diffor =diffor 1064 results_gs%residm =residm 1065 results_gs%res2 =res2 1066 results_gs%fcart(:,:) =fcart(:,:) 1067 results_gs%gred(:,:) =gred(:,:) 1068 results_gs%grchempottn(:,:)=grchempottn(:,:) 1069 results_gs%gresid(:,:)=gresid(:,:) 1070 results_gs%grewtn(:,:)=grewtn(:,:) 1071 results_gs%grxc(:,:) =grxc(:,:) 1072 results_gs%berryopt =dtset%berryopt 1073 results_gs%pel(1:3) =pel(1:3) 1074 results_gs%pion(1:3) =pion(1:3) 1075 results_gs%strten(1:6)=strten(1:6) 1076 results_gs%synlgr(:,:)=synlgr(:,:) 1077 results_gs%vxcavg =vxcavg 1078 if (ngrvdw>0) results_gs%grvdw(1:3,1:ngrvdw)=grvdw(1:3,1:ngrvdw) 1079 if (associated(extfpmd)) then 1080 results_gs%entropy_extfpmd=extfpmd%entropy 1081 results_gs%nelect_extfpmd=extfpmd%nelect 1082 results_gs%shiftfactor_extfpmd=extfpmd%shiftfactor 1083 end if 1084 1085 results_gs%intgres(:,:)=zero 1086 results_gs%grcondft(:,:)=zero 1087 if(any(dtset%constraint_kind(:)/=0))then 1088 results_gs%intgres(1:dtset%nspden,:) =intgres(1:dtset%nspden,:) 1089 results_gs%grcondft(:,:) =grcondft(:,:) 1090 endif 1091 1092 if (dtset%nstep == 0 .and. dtset%occopt>=3.and.dtset%occopt<=8) then 1093 results_gs%etotal = results_gs%etotal - dtset%tsmear * results_gs%entropy 1094 end if 1095 1096 !This call is only for testing purpose: 1097 !test of the nonlop routine (DFPT vs Finite Differences) 1098 if (dtset%useria==112233) then 1099 call nonlop_test(cg,eigen,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,dtset%mgfft,dtset%mkmem,& 1100 & mpi_enreg,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,dtset%nkpt,& 1101 & dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,paw_ij,& 1102 & pawtab,ph1d,psps,rprimd,dtset%typat,xred) 1103 end if 1104 1105 call timab(257,2,tsec) 1106 call timab(250,2,tsec) 1107 1108 DBG_EXIT("COLL") 1109 1110 #if !defined HAVE_BIGDFT 1111 if (.false.) write(std_out,*) vtrial(1,1) 1112 #endif 1113 1114 end subroutine afterscfloop
ABINIT/m_afterscfloop [ Modules ]
NAME
m_afterscfloop
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (XG) 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_afterscfloop 23 24 use defs_basis 25 use defs_wvltypes 26 use m_energies 27 use m_errors 28 use m_abicore 29 use m_efield 30 use m_ab7_mixing 31 use m_hdr 32 use m_dtset 33 use m_dtfil 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_xmpi, only : xmpi_sum, xmpi_comm_rank,xmpi_comm_size 40 use m_berryphase_new, only : berryphase_new 41 use m_geometry, only : xred2xcart, metric 42 use m_crystal, only : prtposcar 43 use m_results_gs , only : results_gs_type 44 use m_electronpositron, only : electronpositron_type, electronpositron_calctype, exchange_electronpositron 45 use m_paw_dmft, only : paw_dmft_type 46 use m_pawang, only : pawang_type 47 use m_pawrad, only : pawrad_type 48 use m_pawtab, only : pawtab_type 49 use m_pawrhoij, only : pawrhoij_type 50 use m_paw_an, only : paw_an_type 51 use m_paw_ij, only : paw_ij_type 52 use m_pawfgrtab, only : pawfgrtab_type 53 use m_pawcprj, only : pawcprj_type,pawcprj_getdim 54 use m_pawfgr, only : pawfgr_type 55 use m_paw_mkrho, only : pawmkrho 56 use m_paw_nhat, only : nhatgrid,wvl_nhatgrid 57 use m_paw_occupancies, only : pawmkrhoij 58 use m_paw_correlations, only : setnoccmmp 59 use m_fock, only : fock_type 60 use m_kg, only : getph 61 use m_spin_current, only : spin_current 62 use m_mkrho, only : mkrho, prtrhomxmn 63 use m_elpolariz, only : elpolariz 64 use m_orbmag, only : orbmag 65 use m_nonlop_test, only : nonlop_test 66 use m_common, only : scprqt 67 use m_xctk, only : xcden 68 use m_forstr, only : forstr 69 use m_wvl_rho, only : wvl_mkrho 70 use m_wvl_psi, only : wvl_psitohpsi, wvl_tail_corrections 71 use m_fourier_interpol, only : transgrid 72 73 #ifdef HAVE_BIGDFT 74 use m_abi2big 75 use BigDFT_API, only : last_orthon, & 76 & kswfn_free_scf_data, denspot_free_history,& 77 & write_energies, total_energies, XC_potential,& 78 & eigensystem_info, applyprojectorsonthefly 79 #endif 80 81 implicit none 82 83 private