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.
COPYRIGHT
Copyright (C) 2005-2018 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 .
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) 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) irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istep=number of the SCF iteration 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=informations 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 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 vxcavg=vxc average xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 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 fred, this array has been corrected by enforcing the translational symmetry, namely that the sum of force on all atoms is zero. fred(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. dtorbmag <type(orbmag_type)> = variables related to orbital magnetization 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_nonlocalpsp(IN)=nonlocal pseudopotential 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
PARENTS
scfcv
CHILDREN
applyprojectorsonthefly,chern_number,denspot_free_history eigensystem_info,elpolariz,energies_copy,exchange_electronpositron forstr,getph,hdr_update,kswfn_free_scf_data,last_orthon,metric,mkrho nhatgrid,nonlop_test,pawcprj_getdim,pawmkrho,pawmkrhoij,prtposcar prtrhomxmn,scprqt,setnoccmmp,spin_current,timab,total_energies write_energies,wrtout,wvl_eigen_abi2big,wvl_mkrho,wvl_nhatgrid wvl_occ_abi2big,wvl_psitohpsi,wvl_rho_abi2big,wvl_tail_corrections wvl_vtrial_abi2big,xcden,xmpi_sum,xred2xcart
SOURCE
198 #if defined HAVE_CONFIG_H 199 #include "config.h" 200 #endif 201 202 #include "abi_common.h" 203 204 subroutine afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,& 205 & deltae,diffor,dtefield,dtfil,dtorbmag,dtset,eigen,electronpositron,elfr,& 206 & energies,etotal,favg,fcart,fock,forold,fred,grchempottn,& 207 & gresid,grewtn,grhf,grhor,grvdw,& 208 & grxc,gsqcut,hdr,indsym,irrzon,istep,kg,kxc,lrhor,maxfor,mcg,mcprj,mgfftf,& 209 & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,& 210 & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,& 211 & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,prtxml,& 212 & psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,& 213 & rhog,rhor,rprimd,stress_needed,strsxc,strten,symrec,synlgr,taug,& 214 & taur,tollist,usecprj,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,& 215 & xccc3d,xred,ylm,ylmgr,qvpotzero,conv_retcode) 216 217 use defs_basis 218 use defs_datatypes 219 use defs_abitypes 220 use defs_parameters 221 use defs_wvltypes 222 use m_energies 223 use m_errors 224 use m_profiling_abi 225 use m_efield 226 use m_orbmag 227 use m_ab7_mixing 228 use m_hdr 229 230 use m_xmpi, only : xmpi_sum, xmpi_comm_rank,xmpi_comm_size 231 use m_geometry, only : xred2xcart, metric 232 use m_crystal, only : prtposcar 233 use m_results_gs , only : results_gs_type 234 use m_electronpositron, only : electronpositron_type,electronpositron_calctype,exchange_electronpositron 235 use m_dtset, only : dtset_copy, dtset_free 236 use m_paw_dmft, only : paw_dmft_type 237 use m_pawang, only : pawang_type 238 use m_pawrad, only : pawrad_type 239 use m_pawtab, only : pawtab_type 240 use m_pawrhoij, only : pawrhoij_type 241 use m_paw_an, only : paw_an_type 242 use m_paw_ij, only : paw_ij_type 243 use m_pawfgrtab, only : pawfgrtab_type 244 use m_pawcprj, only : pawcprj_type,pawcprj_getdim 245 use m_pawfgr, only : pawfgr_type 246 use m_fock, only : fock_type 247 use m_kg, only : getph 248 249 #ifdef HAVE_BIGDFT 250 use m_abi2big 251 use BigDFT_API, only : last_orthon, & 252 & kswfn_free_scf_data, denspot_free_history,& 253 & write_energies, total_energies, XC_potential,& 254 & eigensystem_info, applyprojectorsonthefly 255 #endif 256 257 !This section has been created automatically by the script Abilint (TD). 258 !Do not modify the following lines by hand. 259 #undef ABI_FUNC 260 #define ABI_FUNC 'afterscfloop' 261 use interfaces_14_hidewrite 262 use interfaces_18_timing 263 use interfaces_56_xc 264 use interfaces_62_wvl_wfs 265 use interfaces_65_paw 266 use interfaces_66_nonlocal 267 use interfaces_67_common 268 use interfaces_94_scfcv, except_this_one => afterscfloop 269 !End of the abilint section 270 271 implicit none 272 273 !Arguments ------------------------------------ 274 !scalars 275 integer,intent(in) :: istep,mcg,mcprj,mgfftf,moved_atm_inside,my_natom,n3xccc,nfftf,ngrvdw,nkxc 276 integer,intent(in) :: optres,prtfor,prtxml,pwind_alloc,stress_needed,usecprj 277 integer,intent(inout) :: computed_forces 278 real(dp),intent(in) :: cpus,deltae,gsqcut,res2,residm 279 real(dp),intent(in) :: qvpotzero 280 real(dp),intent(inout) :: diffor,etotal,maxfor,vxcavg 281 type(MPI_type),intent(inout) :: mpi_enreg 282 type(datafiles_type),intent(in) :: dtfil 283 type(dataset_type),intent(inout) :: dtset 284 type(efield_type),intent(inout) :: dtefield 285 type(orbmag_type),intent(inout) :: dtorbmag 286 type(electronpositron_type),pointer :: electronpositron 287 type(energies_type),intent(inout) :: energies 288 type(hdr_type),intent(inout) :: hdr 289 type(pawang_type),intent(in) :: pawang 290 type(pawfgr_type),intent(in) :: pawfgr 291 type(pseudopotential_type),intent(in) :: psps 292 type(results_gs_type),intent(inout) :: results_gs 293 type(wvl_data),intent(inout) :: wvl 294 type(fock_type),pointer, intent(inout) :: fock 295 !arrays 296 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom) 297 integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom) 298 integer,intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 299 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%ntypat) 300 integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt) 301 integer,intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym) 302 integer,intent(out) :: conv_retcode 303 real(dp),intent(in) :: grchempottn(3,dtset%natom),grewtn(3,dtset%natom),grvdw(3,ngrvdw) 304 real(dp),intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 305 real(dp),intent(in) :: pwnsfac(2,pwind_alloc) 306 real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 307 real(dp),intent(in) :: rprimd(3,3),tollist(12),vpsp(nfftf) 308 real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden) 309 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 310 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 311 real(dp),intent(inout) :: cg(2,mcg) 312 real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 313 real(dp),intent(inout) :: forold(3,dtset%natom) 314 real(dp),intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw) 315 real(dp),intent(inout) :: nvresid(nfftf,dtset%nspden),pel(3) 316 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),pel_cg(3) 317 real(dp),intent(inout) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 318 real(dp),intent(inout) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),pion(3) 319 real(dp),intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),strsxc(6) 320 real(dp),intent(inout) :: vhartr(nfftf),vxc(nfftf,dtset%nspden) 321 real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom) 322 real(dp),intent(inout) :: favg(3),fcart(3,dtset%natom),fred(3,dtset%natom) 323 real(dp),intent(inout) :: gresid(3,dtset%natom),grhf(3,dtset%natom) 324 real(dp),intent(inout) :: grxc(3,dtset%natom),kxc(nfftf,nkxc),strten(6) 325 real(dp),intent(inout) :: synlgr(3,dtset%natom) 326 real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taug(:,:),taur(:,:) 327 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj) 328 type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw) 329 type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw) 330 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 331 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 332 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw) 333 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 334 335 !Local variables------------------------------- 336 !scalars 337 integer,parameter :: response=0 338 integer :: bantot,bufsz,choice,cplex,ierr,ifft,igrad,ishift,ispden,nfftotf,ngrad 339 integer :: optcut,optfor,optgr0,optgr1,optgr2,optrad,quit,shft 340 integer :: spaceComm_fft,tim_mkrho 341 logical :: test_gylmgr,test_nfgd,test_rfgd 342 logical :: wvlbigdft=.false. 343 real(dp) :: c_fermi,dtaur,dtaurzero 344 real(dp) :: ucvol 345 character(len=500) :: message 346 type(paw_dmft_type) :: paw_dmft 347 #if defined HAVE_BIGDFT 348 integer :: ia,ii,mband_cprj 349 logical :: do_last_ortho,compute_wvl_tail=.false. 350 real(dp) :: dum,eexctx,eh,ekin,eloc,enl,eproj,esicdc,evxc,exc,ucvol_local 351 #endif 352 !arrays 353 real(dp) :: gmet(3,3),gprimd(3,3),pelev(3),rmet(3,3),tsec(2) 354 real(dp) :: dmatdum(0,0,0,0) 355 real(dp),allocatable :: mpibuf(:,:),qphon(:),rhonow(:,:,:),sqnormgrhor(:,:) 356 #if defined HAVE_BIGDFT 357 integer,allocatable :: dimcprj_srt(:) 358 real(dp),allocatable :: hpsi_tmp(:),xcart(:,:) 359 #endif 360 361 ! ************************************************************************* 362 363 DBG_ENTER("COLL") 364 365 call timab(250,1,tsec) 366 call timab(251,1,tsec) 367 368 !Compute different geometric tensor, as well as ucvol, from rprimd 369 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 370 nfftotf=product(ngfftf(1:3)) 371 372 !MPI FFT communicator 373 spaceComm_fft=mpi_enreg%comm_fft 374 375 !Recompute structure factor phases if atomic positions have changed 376 if (moved_atm_inside==1) then 377 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 378 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 379 else 380 ph1d(:,:)=ph1df(:,:) 381 end if 382 end if 383 384 !---------------------------------------------------------------------- 385 !Wavelet case: transform psi to KS orbitals 386 !---------------------------------------------------------------------- 387 if (dtset%usewvl == 1) then 388 389 ! wvlbigdft indicates that the BigDFT workflow will be followed 390 wvlbigdft=(dtset%wvl_bigdft_comp==1) 391 392 #if defined HAVE_BIGDFT 393 ! Transform to KS orbitals 394 395 ! Need xcart 396 ABI_ALLOCATE(xcart,(3, dtset%natom)) 397 call xred2xcart(dtset%natom, rprimd, xcart, xred) 398 ucvol_local=product(wvl%den%denspot%dpbox%hgrids)*real(product(wvl%den%denspot%dpbox%ndims),dp) 399 400 ! do_last_ortho in case of direct minimization, since 401 ! We never diagonalized the hamiltonian and the eigenvalues are unknown. 402 if ( wvlbigdft) do_last_ortho=(dtset%iscf==0) 403 if (.not.wvlbigdft) do_last_ortho=(.false.) 404 if (do_last_ortho) then 405 call total_energies(wvl%e%energs, istep, mpi_enreg%me_wvl) 406 call write_energies(istep,0,wvl%e%energs,zero,zero,"FINAL") 407 if(.not.wvlbigdft) then 408 ! If ISCF>10, we exit scfcv at a place where bigdft objects 409 ! do not contain the KS potential. Hence, we copy vtrial to wvl%den 410 if(dtset%iscf>=10) call wvl_vtrial_abi2big(1,vtrial,wvl%den) 411 ! hpsi is lost in hpsitopsi, so we recalculate it (needed for last_orthon). 412 call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,& 413 & istep,1,-1,mpi_enreg%me_wvl,dtset%natom,& 414 & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,& 415 & dum,.false.,evxc,wvl,wvlbigdft,xcart,strsxc) 416 if (dtset%iscf==0) then 417 energies%e_kinetic=ekin ; energies%e_hartree=eh 418 energies%e_xc=exc ; energies%e_localpsp=eloc 419 energies%e_nonlocalpsp=enl ; energies%e_exactX=eexctx 420 energies%e_sicdc=esicdc ; energies%e_xcdc=evxc 421 energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp & 422 & + energies%e_xcdc + two*energies%e_hartree +energies%e_nonlocalpsp 423 end if 424 end if 425 call last_orthon(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,istep,wvl%wfs%ks,wvl%e%energs%evsum,.true.) 426 if (mpi_enreg%nproc_wvl == 1) nullify(wvl%wfs%ks%psit) 427 call eigensystem_info(mpi_enreg%me_wvl,mpi_enreg%nproc_wvl,0.d0,& 428 wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,& 429 wvl%wfs%ks%orbs,wvl%wfs%ks%psi) 430 ! Copy eigenvalues from BigDFT object to "eigen" 431 call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs) 432 ! Copy occupations from BigDFT objects to ABINIT 433 call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs) 434 end if 435 436 ! Tail corrections, pending for wvlbigdft==.false. 437 ! TODO put it at the end of gstate. 438 ! WVL - maybe compute the tail corrections to energy 439 compute_wvl_tail=(dtset%tl_radius>tol12.and.wvlbigdft) 440 if (compute_wvl_tail) then 441 ! Use the tails to improve energy precision. 442 call wvl_tail_corrections(dtset, energies, etotal, mpi_enreg, psps, wvl, xcart) 443 end if 444 445 ! Clean KSwfn parts only needed in the SCF loop. 446 call kswfn_free_scf_data(wvl%wfs%ks, (mpi_enreg%nproc_wvl > 1)) 447 ! Clean denspot parts only needed in the SCF loop. 448 call denspot_free_history(wvl%den%denspot) 449 450 ! If WF have been modified, change the density according to the KS projection. 451 if ( do_last_ortho ) then 452 453 ! Density from new orthogonalized WFs 454 call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den) 455 456 ! PAW: has to update cprj, rhoij and compensation charge density 457 if (psps%usepaw==1) then 458 ! 1-Compute cprj 459 ABI_ALLOCATE(hpsi_tmp,(size(wvl%wfs%ks%hpsi))) 460 call applyprojectorsonthefly(mpi_enreg%me_wvl,wvl%wfs%ks%orbs,wvl%descr%atoms,wvl%descr%Glr,& 461 & xcart,wvl%descr%h(1),wvl%descr%h(2),wvl%descr%h(3),wvl%wfs%ks%lzd%Glr%wfd,& 462 & wvl%projectors%nlpsp,wvl%wfs%ks%psi,hpsi_tmp,eproj,& 463 & proj_G=wvl%projectors%G,paw=wvl%descr%paw) 464 ABI_DEALLOCATE(hpsi_tmp) 465 do ii=1,mcprj 466 do ia=1,dtset%natom 467 !Note that cprj should be allocated (i.e. usepcrj=1 imposed in scfcv) 468 cprj(ia,ii)%cp(:,:)= wvl%descr%paw%cprj(ia,ii)%cp(:,:) 469 end do 470 end do 471 ! 2-Compute rhoij 472 ABI_ALLOCATE(dimcprj_srt,(dtset%natom)) 473 call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 474 mband_cprj=mcprj/(dtset%nspinor*dtset%mkmem*dtset%nsppol) 475 paw_dmft%use_sc_dmft=0 ; paw_dmft%use_dmft=0 ! dmft not used here 476 call pawmkrhoij(atindx,atindx1,cprj,dimcprj_srt,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,& 477 & mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspinor,dtset%nsppol,& 478 & occ,mpi_enreg%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij,dtfil%unpaw,dtset%usewvl,dtset%wtk) 479 ABI_DEALLOCATE(dimcprj_srt) 480 ! 3-Symetrize rhoij, compute nhat and add it to rhor 481 call pawmkrho(dum,1,gprimd,0,indsym,0,mpi_enreg,my_natom,dtset%natom,dtset%nspden,dtset%nsym,& 482 & dtset%ntypat,mpi_enreg%paral_kgb,pawang,pawfgr,pawfgrtab,dtset%pawprtvol,pawrhoij,pawrhoij,& 483 & pawtab,(/zero,zero,zero/),rhog,rhor,rhor,rprimd,dtset%symafm,symrec,dtset%typat,ucvol_local,& 484 & dtset%usewvl,xred,pawnhat=nhat) 485 call wvl_rho_abi2big(1,rhor,wvl%den) 486 end if 487 end if 488 489 ABI_DEALLOCATE(xcart) 490 491 #else 492 BIGDFT_NOTENABLED_ERROR() 493 #endif 494 end if 495 496 call timab(251,2,tsec) 497 call timab(252,1,tsec) 498 499 !---------------------------------------------------------------------- 500 !Polarization Calculation 501 !---------------------------------------------------------------------- 502 503 if(dtset%berryopt/=0)then 504 call elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,energies%e_elecfield,gprimd,hdr,& 505 & kg,dtset%mband,mcg,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,dtset%nkpt,& 506 & npwarr,dtset%nsppol,psps%ntypat,pawrhoij,pawtab,pel,pel_cg,pelev,pion,& 507 & psps,pwind,pwind_alloc,pwnsfac,rprimd,ucvol,usecprj,xred) 508 end if 509 510 !---------------------------------------------------------------------- 511 ! Orbital magnetization calculations 512 !---------------------------------------------------------------------- 513 if(dtset%orbmag==1) then 514 call chern_number(atindx1,cg,cprj,dtset,dtorbmag,gmet,gprimd,kg,& 515 & mcg,size(cprj,2),mpi_enreg,npwarr,pawang,pawrad,pawtab,pwind,pwind_alloc,& 516 & symrec,usecprj,psps%usepaw,xred) 517 end if 518 519 call timab(252,2,tsec) 520 call timab(253,1,tsec) 521 522 !---------------------------------------------------------------------- 523 !Gradient and Laplacian of the Density Calculation 524 !---------------------------------------------------------------------- 525 526 !We use routine xcden which get gradient of rhor (grhor), and eventually laplacian of rhor (lrhor). 527 if(dtset%prtgden/=0 .or. dtset%prtlden/=0)then 528 529 ! Compute gradient of the electron density 530 ngrad=2 531 cplex=1 532 ishift=0 533 ABI_ALLOCATE(rhonow,(nfftf,dtset%nspden,ngrad*ngrad)) 534 if(dtset%prtlden/=0)then 535 nullify(lrhor) 536 ABI_ALLOCATE(lrhor,(nfftf,dtset%nspden)) 537 end if 538 write(message,'(a,a)') ch10, " Compute gradient of the electron density" 539 call wrtout(ab_out,message,'COLL') 540 if(dtset%prtlden/=0) then 541 write(message,'(a)') " and also Compute Laplacian of the electron density" 542 call wrtout(ab_out,message,'COLL') 543 end if 544 write(message,'(a)') "--------------------------------------------------------------------------------" 545 call wrtout(ab_out,message,'COLL') 546 547 ABI_ALLOCATE(qphon,(3)) 548 qphon(:)=zero 549 if(dtset%prtlden/=0) then 550 call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,dtset%paral_kgb,qphon,rhor,rhonow,lrhonow=lrhor) 551 else 552 call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,dtset%paral_kgb,qphon,rhor,rhonow) 553 end if 554 ABI_DEALLOCATE(qphon) 555 556 ! Copy gradient which has been output in rhonow to grhor (and free rhonow) 557 nullify(grhor) 558 ABI_ALLOCATE(grhor,(nfftf,dtset%nspden,3)) 559 do ispden=1,dtset%nspden 560 do ifft=1,nfftf 561 grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4) 562 end do 563 end do 564 ABI_DEALLOCATE(rhonow) 565 566 if(dtset%prtgden/=0) then 567 ! Print result for grhor 568 write(message,'(a,a)') ch10, " Result for gradient of the electron density for each direction (1,2,3):" 569 call wrtout(ab_out,message,'COLL') 570 write(message,'(a,a,a,a)') ch10," 1rst direction:",ch10,& 571 & "--------------------------------------------------------------------------------" 572 call wrtout(ab_out,message,'COLL') 573 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,1),optrhor=2,ucvol=ucvol) 574 write(message,'(a)') "--------------------------------------------------------------------------------" 575 call wrtout(ab_out,message,'COLL') 576 write(message,'(a,a,a,a)') ch10," 2nd direction:",ch10,& 577 & "--------------------------------------------------------------------------------" 578 call wrtout(ab_out,message,'COLL') 579 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,2),optrhor=2,ucvol=ucvol) 580 write(message,'(a)') "--------------------------------------------------------------------------------" 581 call wrtout(ab_out,message,'COLL') 582 write(message,'(a,a,a,a)') ch10," 3rd direction:",ch10,& 583 & "--------------------------------------------------------------------------------" 584 call wrtout(ab_out,message,'COLL') 585 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,grhor(:,:,3),optrhor=2,ucvol=ucvol) 586 write(message,'(a)') "--------------------------------------------------------------------------------" 587 call wrtout(ab_out,message,'COLL') 588 end if 589 590 if(dtset%prtlden/=0) then 591 ! Print result for lrhor 592 write(message,'(a,a)') ch10, " Result for Laplacian of the electron density :" 593 call wrtout(ab_out,message,'COLL') 594 write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------" 595 call wrtout(ab_out,message,'COLL') 596 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,lrhor,optrhor=3,ucvol=ucvol) 597 write(message,'(a)') "--------------------------------------------------------------------------------" 598 call wrtout(ab_out,message,'COLL') 599 end if 600 601 write(message,'(a)') "--------------------------------------------------------------------------------" 602 call wrtout(ab_out,message,'COLL') 603 end if 604 605 !---------------------------------------------------------------------- 606 !Kinetic Energy Density Calculation 607 !---------------------------------------------------------------------- 608 609 call timab(253,2,tsec) 610 call timab(254,1,tsec) 611 612 !We use routine mkrho with option=1 to compute kinetic energy density taur (and taug) 613 if(dtset%prtkden/=0 .or. dtset%prtelf/=0)then 614 nullify(taug,taur) 615 ! tauX are reused in outscfcv for output 616 ! should be deallocated there 617 ABI_ALLOCATE(taug,(2,nfftf)) 618 ABI_ALLOCATE(taur,(nfftf,dtset%nspden)) 619 tim_mkrho=5 620 if(dtset%prtelf/=0) then 621 write(message,'(a,a)') ch10, " Compute ELF" 622 call wrtout(ab_out,message,'COLL') 623 write(message,'(a)') "--------------------------------------------------------------------------------" 624 call wrtout(ab_out,message,'COLL') 625 end if 626 write(message,'(a,a)') ch10, " Compute kinetic energy density" 627 call wrtout(ab_out,message,'COLL') 628 paw_dmft%use_sc_dmft=0 ! dmft not used here 629 paw_dmft%use_dmft=0 ! dmft not used here 630 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,& 631 & npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 632 ! Print result 633 if(dtset%prtkden/=0) then 634 write(message,'(a,a)') ch10, "Result for kinetic energy density :" 635 call wrtout(ab_out,message,'COLL') 636 write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------" 637 call wrtout(ab_out,message,'COLL') 638 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,taur,optrhor=1,ucvol=ucvol) 639 write(message,'(a)') "--------------------------------------------------------------------------------" 640 call wrtout(ab_out,message,'COLL') 641 end if 642 ABI_DEALLOCATE(taug) 643 end if 644 645 !---------------------------------------------------------------------- 646 !Electron Localization Function (ELF) Calculation 647 !---------------------------------------------------------------------- 648 649 call timab(254,2,tsec) 650 call timab(255,1,tsec) 651 652 !We use routine xcden to compute gradient of electron density (grhor), 653 !NOTE: If GGA is used, gradient of electron density is already computed 654 !and it is stored in exchange correlation kernel kxc(:,5:7) (nspden=1) or kxc(:,14:19) (nspden=2). 655 !In order to save memory and do not have the same quantity twice 656 !in memory we should use kxc. 657 !But unfortunately only spin up ans spin down component are stored in kxc. 658 !So far we thus use grhor which contains all component (nspden=4) 659 !just like rhonow in xcden instead of kxc. 660 661 if((dtset%prtelf/=0))then 662 if(dtset%nspden<=2) then 663 664 ngrad=2 665 cplex=1 666 if((cplex*dtset%nfft)/=nfftf)then 667 write(message, '(a,a,a,a)' ) ch10,& 668 & ' afterscfloop: ERROR -', ch10, & 669 & ' The density is complex, ELF analysis cannot be performed.' 670 call wrtout(std_out,message,'COLL') 671 ! MSG_ERROR(message) 672 end if 673 674 if((dtset%prtgden==0) .and. (dtset%prtlden==0)) then 675 ! Compute gradient of the electron density 676 ishift=0 677 ABI_ALLOCATE(rhonow,(nfftf,dtset%nspden,ngrad*ngrad)) 678 write(message,'(a,a)') ch10, " Compute gradient of the electron density" 679 call wrtout(ab_out,message,'COLL') 680 ABI_ALLOCATE(qphon,(3)) 681 qphon(:)=zero 682 call xcden (cplex,gprimd,ishift,mpi_enreg,nfftf,ngfftf,ngrad,dtset%nspden,dtset%paral_kgb,qphon,rhor,rhonow) 683 ABI_DEALLOCATE(qphon) 684 ! Copy gradient which has been output in rhonow to grhor (and free rhonow) 685 ABI_ALLOCATE(grhor,(nfftf,dtset%nspden,3)) 686 do ispden=1,dtset%nspden 687 do ifft=1,nfftf 688 grhor(ifft,ispden,1:3) = rhonow(ifft,ispden,2:4) 689 end do 690 end do 691 ABI_DEALLOCATE(rhonow) 692 end if 693 ! Compute square norm of gradient of the electron density (|grhor|**2) 694 if(dtset%nspden==1)then 695 ABI_ALLOCATE(sqnormgrhor,(nfftf,dtset%nspden)) 696 do ifft=1,nfftf 697 sqnormgrhor(ifft,1) = zero 698 end do 699 elseif(dtset%nspden==2)then 700 ABI_ALLOCATE(sqnormgrhor,(nfftf,dtset%nspden+1)) 701 ! because we not only want (total and up quantities, but also down) 702 ! Indeed after having token the square norm we can not recover the 703 ! down quantity by substracting total and up quantities (as we do for usual densities) 704 do ispden=1,dtset%nspden+1 705 do ifft=1,nfftf 706 sqnormgrhor(ifft,ispden) = zero 707 end do 708 end do 709 end if 710 711 do igrad=1,3 712 do ispden=1,dtset%nspden !total (and up) 713 do ifft=1,nfftf 714 sqnormgrhor(ifft,ispden) = sqnormgrhor(ifft,ispden) + grhor(ifft,ispden,igrad)**2 715 end do 716 end do 717 if(dtset%nspden==2)then 718 do ifft=1,nfftf !down 719 sqnormgrhor(ifft,3) = sqnormgrhor(ifft,3) + (grhor(ifft,1,igrad)-grhor(ifft,2,igrad))**2 720 end do 721 end if 722 end do 723 724 ! Compute electron localization function (ELF) (here it is elfr) 725 726 nullify(elfr) 727 if(dtset%nspden==1)then 728 ABI_ALLOCATE(elfr,(nfftf,dtset%nspden)) 729 elseif(dtset%nspden==2)then 730 ABI_ALLOCATE(elfr,(nfftf,dtset%nspden+1)) 731 ! 1rst is total elf, 2nd is spin-up elf, and 3rd is spin-down elf. (elf_tot /= elf_up + elf_down) 732 end if 733 c_fermi = 3.0d0/10.0d0*((3.0d0*pi**2)**(2.0d0/3.0d0)) 734 735 ! First compute total elf 736 ispden=1 737 do ifft=1,nfftf 738 dtaurzero = c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0) 739 dtaur = taur(ifft,ispden) 740 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden)) 741 ! Ensure that dtaur is always positive or zero, as it should be. 742 if(dtaur<0.0d0)dtaur=0.0d0 743 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 744 if(dtaurzero<(1.0d-20*dtaur)) then 745 elfr(ifft,ispden) = 0.0d0 746 else 747 elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 748 ! For atomic shell studies we could also add the condition that when (dtaur/dtaurzero) 749 ! is very close to zero we actually set it exactly to zero (or --> elfr = 1.0d0 750 ! which is the upper limit of elfr.) 751 end if 752 end do 753 754 ! If spin-dependent densities are avalaible, compute spin-dependent elf 755 ! and offer the possibility to compute total elf in an alternative approach 756 ! (see doc/theory/ELF) 757 if(dtset%nspden==2)then 758 759 ! alternative approach to the total elf 760 if(dtset%prtelf==2)then 761 ispden=1 762 do ifft=1,nfftf 763 dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*( rhor(ifft,ispden+1)**(5.0d0/3.0d0) + & 764 & (rhor(ifft,ispden) - rhor(ifft,ispden+1))**(5.0d0/3.0d0) ) 765 dtaur = taur(ifft,ispden) 766 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+1)/rhor(ifft,ispden+1)) & 767 & - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1))) 768 if(dtaur<0.0d0)dtaur=0.0d0 769 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 770 if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then 771 elfr(ifft,ispden) = 0.0d0 772 else 773 elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 774 end if 775 end do 776 end if 777 778 ! elf_up 779 ispden=2 780 do ifft=1,nfftf 781 dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*rhor(ifft,ispden)**(5.0d0/3.0d0) 782 dtaur = taur(ifft,ispden) 783 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden)/rhor(ifft,ispden)) 784 if(dtaur<0.0d0)dtaur=0.0d0 785 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 786 if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then 787 elfr(ifft,ispden) = 0.0d0 788 else 789 elfr(ifft,ispden) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 790 end if 791 end do 792 793 ! elf_down 794 ispden=1 795 do ifft=1,nfftf 796 dtaurzero = 2.0d0**(2.0d0/3.0d0)*c_fermi*(rhor(ifft,ispden)-rhor(ifft,ispden+1))**(5.0d0/3.0d0) 797 dtaur = taur(ifft,ispden)-taur(ifft,ispden+1) 798 dtaur = dtaur - (1.0d0/8.0d0)*(sqnormgrhor(ifft,ispden+2)/(rhor(ifft,ispden)-rhor(ifft,ispden+1))) 799 if(dtaur<0.0d0)dtaur=0.0d0 800 ! To avoid NaN values we check that dtaurzero is not to small compare to dtaur 801 if(abs(dtaurzero)<abs(1.0d-20*dtaur)) then 802 elfr(ifft,ispden+2) = 0.0d0 803 else 804 elfr(ifft,ispden+2) = 1.0d0/(1.0d0 + (dtaur/dtaurzero)**2) 805 end if 806 end do 807 808 end if !endif dtset%nspden==2 809 810 ! Print result for elfr 811 call prtrhomxmn(ab_out,mpi_enreg,nfftf,ngfftf,dtset%nspden,1,elfr,optrhor=4,ucvol=ucvol) 812 813 ABI_DEALLOCATE(grhor) 814 ABI_DEALLOCATE(sqnormgrhor) 815 816 else 817 message ='ELF is not yet implemented for non collinear spin cases.' 818 MSG_WARNING(message) 819 820 ABI_ALLOCATE(elfr,(nfftf,dtset%nspden)) 821 do ispden=1,dtset%nspden 822 do ifft=1,nfftf 823 elfr(ifft,ispden) = -2.0d0 824 end do 825 end do 826 ! even if elf is not computed we want to finish the abinit run. 827 ! To ensure that users won't use the _ELF file which will be produced 828 ! we set elf to -2.0 (a meaningless value) 829 830 end if ! endif dtset%nspden<=2 831 832 write(message,'(a,a)') ch10, "--------------------------------------------------------------------------------" 833 call wrtout(ab_out,message,'COLL') 834 write(message,'(a)') " End of ELF section" 835 call wrtout(ab_out,message,'COLL') 836 837 end if !endif prtelf/=0 838 839 !###################################################################### 840 !Compute forces (if they were not computed during the elec. iterations) 841 !and stresses (if requested by user) 842 !---------------------------------------------------------------------- 843 844 call timab(255,2,tsec) 845 call timab(256,1,tsec) 846 847 optfor=0 848 849 if (computed_forces==0.and.dtset%optforces>0.and.dtset%iscf>=0) then 850 if (dtset%nstep>0.or.dtfil%ireadwf==1) optfor=1 851 end if 852 853 if (optfor>0.or.stress_needed>0) then 854 855 ! PAW: eventually, compute g_l(r).Y_lm(r) gradients (if not already done) 856 if (psps%usepaw==1) then 857 test_nfgd =any(pawfgrtab(:)%nfgd==0) 858 test_rfgd =any(pawfgrtab(:)%rfgd_allocated==0) 859 test_gylmgr=any(pawfgrtab(:)%gylmgr_allocated==0) 860 if (test_nfgd.or.& 861 & (test_gylmgr.and.dtset%pawstgylm==1).or.& 862 & (test_rfgd.and.stress_needed==1.and.dtset%pawstgylm==1).or.& 863 (test_rfgd.and.dtset%pawstgylm==0)) then 864 optcut=0;optgr0=0;optgr1=dtset%pawstgylm;optgr2=0 865 optrad=1-dtset%pawstgylm;if (stress_needed==1) optrad=1 866 if (dtset%usewvl==0) then 867 call nhatgrid(atindx1,gmet,my_natom,dtset%natom,nattyp,ngfftf,dtset%ntypat,& 868 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,& 869 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 870 & comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft) 871 else 872 shft=0 873 #if defined HAVE_BIGDFT 874 shft=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,4) 875 call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,& 876 & wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,& 877 & nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,& 878 & wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,& 879 & wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,& 880 & pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shft,xred) 881 #endif 882 end if 883 end if 884 end if 885 886 call forstr(atindx1,cg,cprj,diffor,dtefield,dtset,& 887 & eigen,electronpositron,energies,favg,fcart,fock,& 888 & forold,fred,grchempottn,gresid,grewtn,& 889 & grhf,grvdw,grxc,gsqcut,indsym,& 890 & kg,kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,& 891 & n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,& 892 & npwarr,dtset%ntypat,nvresid,occ,optfor,optres,& 893 & paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1df,& 894 & psps,rhog,rhor,rprimd,stress_needed,strsxc,strten,symrec,synlgr,& 895 & ucvol,usecprj,vhartr,vpsp,vxc,wvl,xccc3d,xred,ylm,ylmgr,qvpotzero) 896 end if 897 898 if (optfor==1) computed_forces=1 899 if (optfor==1) diffor=9.9999999999D99 900 if (stress_needed==0) strten(:)=9.9999999999D99 901 if (computed_forces==0) fcart(:,:)=9.9999999999D99 902 if (dtset%prtstm/=0) strten(:)=zero 903 904 call timab(256,2,tsec) 905 call timab(257,1,tsec) 906 907 !If SCF convergence was not reached (for dtset%nstep>0), 908 !print a warning to the output file (non-dummy arguments: dtset%nstep, 909 !residm, diffor - infos from tollist have been saved inside ) 910 choice=3 911 call scprqt(choice,cpus,deltae,diffor,dtset,& 912 & eigen,etotal,favg,fcart,energies%e_fermie,dtfil%fnameabo_app_eig,dtfil%filnam_ds(1),& 913 & 1,dtset%iscf,istep,dtset%kptns,maxfor,& 914 & moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,& 915 & dtset%nstep,occ,optres,prtfor,prtxml,quit,& 916 & res2,resid,residm,response,tollist,psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,& 917 & electronpositron=electronpositron, fock=fock) 918 919 !output POSCAR and FORCES files, VASP style, for PHON code and friends. 920 if (dtset%prtposcar == 1) then 921 call prtposcar(fcart, dtfil%filnam_ds(4), dtset%natom, dtset%ntypat, rprimd, dtset%typat, ucvol, xred, dtset%znucl) 922 end if ! prtposcar 923 924 if(allocated(qphon)) then 925 ABI_DEALLOCATE(qphon) 926 end if 927 928 !get current operator on wavefunctions 929 if (dtset%prtspcur == 1) then 930 call spin_current(cg,dtfil,dtset,gprimd,hdr,kg,mcg,mpi_enreg,psps) 931 end if 932 933 !Electron-positron stuff: if last calculation was a positron minimization, 934 !exchange electron and positron data in order to 935 !get electronic quantities in global variables 936 if (dtset%positron/=0) then 937 electronpositron%scf_converged=.false. 938 if (dtset%positron<0.and.electronpositron_calctype(electronpositron)==1) then 939 call exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,fred,mcg,mcprj,& 940 & mpi_enreg,my_natom,nfftf,ngfftf,nhat,npwarr,occ,paw_an,pawrhoij,rhog,rhor,strten,usecprj,vhartr) 941 end if 942 end if 943 944 !If PAW+U and density mixing, has to update nocc_mmp 945 if (psps%usepaw==1.and.dtset%usepawu>0.and.(dtset%iscf>0.or.dtset%iscf==-3)) then 946 call setnoccmmp(1,0,dmatdum,0,0,indsym,my_natom,dtset%natom,dtset%natpawu,& 947 & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,& 948 & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,& 949 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 950 end if 951 952 !Update the content of the header (evolving variables) 953 bantot=hdr%bantot 954 if (dtset%positron==0) then 955 call hdr_update(hdr,bantot,etotal,energies%e_fermie,residm,rprimd,occ,& 956 & pawrhoij,xred,dtset%amu_orig(:,1),& 957 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 958 else 959 call hdr_update(hdr,bantot,electronpositron%e0,energies%e_fermie,residm,rprimd,occ,& 960 & pawrhoij,xred,dtset%amu_orig(:,1),& 961 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 962 end if 963 964 #ifdef HAVE_LOTF 965 if(dtset%ionmov==23 .and. mpi_enreg%nproc_band>1) then 966 bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom 967 ABI_ALLOCATE(mpibuf,(3,bufsz)) 968 mpibuf(:,1:dtset%natom)=fred(:,1:dtset%natom) 969 mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom) 970 if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom) 971 mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/)) 972 call xmpi_sum(mpibuf,mpi_enreg%comm_band,ierr) 973 fred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_band 974 fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_band 975 if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_band 976 strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_band 977 ABI_DEALLOCATE(mpibuf) 978 end if 979 #endif 980 981 !In case of FFT parallelisation, has to synchronize positions and forces 982 !to avoid numerical noise 983 if (mpi_enreg%nproc_fft>1) then 984 bufsz=2+2*dtset%natom;if (moved_atm_inside==1) bufsz=bufsz+dtset%natom 985 ABI_ALLOCATE(mpibuf,(3,bufsz)) 986 mpibuf(:,1:dtset%natom)=fred(:,1:dtset%natom) 987 mpibuf(:,dtset%natom+1:2*dtset%natom)=fcart(:,1:dtset%natom) 988 if (moved_atm_inside==1) mpibuf(:,2*dtset%natom+1:3*dtset%natom)=xred(:,1:dtset%natom) 989 mpibuf(1:3,bufsz-1:bufsz)=reshape(strten(1:6),(/3,2/)) 990 call xmpi_sum(mpibuf,mpi_enreg%comm_fft,ierr) 991 fred(:,1:dtset%natom)=mpibuf(:,1:dtset%natom)/mpi_enreg%nproc_fft 992 fcart(:,1:dtset%natom)=mpibuf(:,dtset%natom+1:2*dtset%natom)/mpi_enreg%nproc_fft 993 if (moved_atm_inside==1) xred(:,1:dtset%natom)=mpibuf(:,2*dtset%natom+1:3*dtset%natom)/mpi_enreg%nproc_fft 994 strten(1:6)=reshape(mpibuf(1:3,bufsz-1:bufsz),(/6/))/mpi_enreg%nproc_fft 995 ABI_DEALLOCATE(mpibuf) 996 end if 997 998 !results_gs%energies = energies 999 call energies_copy(energies,results_gs%energies) 1000 results_gs%etotal =etotal 1001 results_gs%deltae =deltae 1002 results_gs%diffor =diffor 1003 results_gs%residm =residm 1004 results_gs%res2 =res2 1005 results_gs%fcart(:,:) =fcart(:,:) 1006 results_gs%fred(:,:) =fred(:,:) 1007 results_gs%grchempottn(:,:)=grchempottn(:,:) 1008 results_gs%gresid(:,:)=gresid(:,:) 1009 results_gs%grewtn(:,:)=grewtn(:,:) 1010 results_gs%grxc(:,:) =grxc(:,:) 1011 results_gs%pel(1:3) =pel(1:3) 1012 results_gs%strten(1:6)=strten(1:6) 1013 results_gs%synlgr(:,:)=synlgr(:,:) 1014 results_gs%vxcavg =vxcavg 1015 if (ngrvdw>0) results_gs%grvdw(1:3,1:ngrvdw)=grvdw(1:3,1:ngrvdw) 1016 if (dtset%nstep == 0 .and. dtset%occopt>=3.and.dtset%occopt<=8) then 1017 results_gs%etotal = results_gs%etotal - dtset%tsmear * results_gs%entropy 1018 end if 1019 1020 !This call is only for testing purpose: 1021 !test of the nonlop routine (DFPT vs Finite Differences) 1022 if (dtset%useria==112233) then 1023 call nonlop_test(cg,eigen,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,dtset%mgfft,dtset%mkmem,& 1024 & mpi_enreg,dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,dtset%ngfft,dtset%nkpt,& 1025 & dtset%nloalg,npwarr,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,paw_ij,pawtab,& 1026 & ph1d,psps,rprimd,dtset%typat,xred) 1027 end if 1028 1029 call timab(257,2,tsec) 1030 call timab(250,2,tsec) 1031 1032 DBG_EXIT("COLL") 1033 1034 #if !defined HAVE_BIGDFT 1035 if (.false.) write(std_out,*) vtrial(1,1) 1036 #endif 1037 1038 end subroutine afterscfloop