TABLE OF CONTENTS
ABINIT/energy [ Functions ]
NAME
energy
FUNCTION
Compute electronic energy terms energies%e_eigenvalues, ek and enl from arbitrary (orthonormal) provided wf, ehart, enxc, and eei from provided density and potential, energies%e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree) ek=kinetic energy, ehart=Hartree electron-electron energy, enxc,enxcdc=exchange-correlation energies, eei=local pseudopotential energy, enl=nonlocal pseudopotential energy Also, compute new density from provided wfs, after the evaluation of ehart, enxc, and eei.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, MB, MT, EB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
[add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy cg(2,mpw*nspinor*mband*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of wavefunction dtset <type(dataset_type)>=all input variables for this dataset | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem=number of k points treated by this node. | mpw=maximum dimension for number of planewaves | natom=number of atoms in unit cell | nfft=(effective) number of FFT grid points (for this processor) | nkpt=number of k points | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for polarized | nspinor=number of spinorial components | nsym=number of symmetry elements in space group (at least 1) | occopt=option for occupancies | tsmear=smearing energy or temperature (if metal) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) gsqcut=G^2 cutoff from gsqcut=ecut/(2 Pi^2) indsym(4,nsym,natom)=indirect indexing array for atom labels irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs) ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid (ngfftf=ngfft for norm-conserving potential runs) nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise npwarr(nkpt)=number of planewaves at each k point, and boundary n3xccc=dimension of the xccc3d array (0 or nfftf). occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2) at each k point optene=option for the computation of total energy (direct scheme or double-counting scheme) paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawfgr(natom) <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases ph1d(2,3*(2*mgfft+1)*natom)=phase information related to structure factor psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials | ntypat=number of types of atoms in cell rprimd(3,3)=dimensional real space primitive translations (bohr) symrec(3,3,nsym)=symmetry operations in reciprocal space usexcnhat= -PAW only- flag controling use of compensation density in Vxc vpsp(nfftf)=local pseudopotential in real space (hartree) wfs <type(wvl_projector_type)>=wavefunctions informations for wavelets. wvl <type(wvl_internal_type)>=wavelets internal data xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) xred(3,natom)=reduced coordinates of atoms (dimensionless) ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
OUTPUT
compch_fft=-PAW only- compensation charge inside spheres computed over fine fft grid etotal=total energy (hartree): - computed by direct scheme if optene=0 or 2 - computed by double-counting scheme if optene=1 or 3 resid(mband*nkpt*nsppol)=residuals for each band over all k points (hartree^2) strsxc(6)=exchange-correlation contribution to stress tensor vhartr(nfftf)=work space to hold Hartree potential in real space (hartree) vtrial(nfftf,nspden)=total local potential (hartree) vxc(nfftf,nspden)=work space to hold Vxc(r) in real space (hartree) [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to kinetic energy density (metaGGA cases) (optional output)
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) energies <type(energies_type)>=all part of total energy. | entropy(IN)=entropy due to the occupation number smearing (if metal) | 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_corepsp(IN)=psp core-core energy | e_paw(IN)=PAW spherical part energy | e_pawdc(IN)=PAW spherical part double-counting energy | e_eigenvalues(OUT)=Sum of the eigenvalues - Band energy (Hartree) | e_hartree(OUT)=Hartree part of total energy (hartree units) | e_kinetic(OUT)=kinetic energy part of total energy. | e_nonlocalpsp(OUT)=nonlocal pseudopotential part of total energy. | e_xc(OUT)=exchange-correlation energy (hartree) ==== if optene==0, 2 or 3 | e_localpsp(OUT)=local psp energy (hartree) ==== if optene==1, 2 or 3 | e_xcdc(OUT)=exchange-correlation double-counting energy (hartree) rhog(2,nfftf)=work space for rho(G); save intact on return (? MT 08-12-2008: is that true now ?) rhor(nfftf,nspden)=work space for rho(r); save intact on return (? MT 08-12-2008: is that true now ?) taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density nspinor should not be modified in the call of rdnpw === if psps%usepaw==1 === nhat(nfftf,nspden*usepaw)= compensation charge density pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
NOTES
Be careful to the meaning of nfft (size of FFT grids): - In case of norm-conserving calculations the FFT grid is the usual FFT grid. - In case of PAW calculations: Two FFT grids are used; one with nfft points (coarse grid) for the computation of wave functions ; one with nfftf points (fine grid) for the computation of total density. There is a large amount of overhead in the way this routine do the computation of the energy ! For example, the density has already been precomputed, so why to compute it again here ??
PARENTS
scfcv
CHILDREN
bandfft_kpt_restoretabs,bandfft_kpt_savetabs,destroy_hamiltonian dotprod_vn,fftpac,fourdp,gpu_finalize_ffnl_ph3d,gpu_update_ffnl_ph3d hartre,init_hamiltonian,load_k_hamiltonian,load_spin_hamiltonian mag_constr,make_gemm_nonlop,meanvalue_g,metric,mkffnl,mkkin,mkresi mkrho,nonlop,pawaccrhoij,pawcprj_alloc,pawcprj_free,pawcprj_gather_spin pawmknhat,pawrhoij_alloc,pawrhoij_free,pawrhoij_free_unpacked pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked,prep_bandfft_tabs prep_nonlop,psolver_rhohxc,rhohxcpositron,rhotoxc,symrhoij,timab transgrid,xcdata_init,xmpi_sum
SOURCE
147 #if defined HAVE_CONFIG_H 148 #include "config.h" 149 #endif 150 151 #include "abi_common.h" 152 153 subroutine energy(cg,compch_fft,dtset,electronpositron,& 154 & energies,eigen,etotal,gsqcut,indsym,irrzon,kg,mcg,mpi_enreg,my_natom,nfftf,ngfftf,nhat,& 155 & nhatgr,nhatgrdim,npwarr,n3xccc,occ,optene,paw_dmft,paw_ij,pawang,pawfgr,& 156 & pawfgrtab,pawrhoij,pawtab,phnons,ph1d,psps,resid,rhog,rhor,rprimd,strsxc,symrec,& 157 & taug,taur,usexcnhat,vhartr,vtrial,vpsp,vxc,vxctau,wfs,wvl,wvl_den,wvl_e,xccc3d,xred,ylm,& 158 & add_tfw) ! optional argument 159 160 161 use defs_basis 162 use defs_datatypes 163 use defs_abitypes 164 use defs_wvltypes 165 use m_profiling_abi 166 use m_hamiltonian 167 use m_errors 168 use m_xmpi 169 use m_gemm_nonlop 170 use m_xcdata 171 172 use m_geometry, only : metric 173 use m_kg, only : mkkin 174 use m_energies, only : energies_type 175 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 176 use m_bandfft_kpt, only : bandfft_kpt,bandfft_kpt_type,& 177 & bandfft_kpt_savetabs,bandfft_kpt_restoretabs 178 use m_pawang, only : pawang_type 179 use m_pawtab, only : pawtab_type 180 use m_paw_ij, only : paw_ij_type 181 use m_pawfgrtab, only : pawfgrtab_type 182 use m_pawrhoij, only : pawrhoij_type,pawrhoij_alloc,pawrhoij_free,pawrhoij_init_unpacked,& 183 & pawrhoij_mpisum_unpacked,pawrhoij_free_unpacked,pawrhoij_get_nspden,symrhoij 184 use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_free,pawcprj_gather_spin 185 use m_pawfgr, only : pawfgr_type 186 use m_paw_dmft, only : paw_dmft_type 187 use m_cgtools, only : dotprod_vn 188 189 !This section has been created automatically by the script Abilint (TD). 190 !Do not modify the following lines by hand. 191 #undef ABI_FUNC 192 #define ABI_FUNC 'energy' 193 use interfaces_18_timing 194 use interfaces_32_util 195 use interfaces_53_ffts 196 use interfaces_53_spacepar 197 use interfaces_56_xc 198 use interfaces_62_poisson 199 use interfaces_65_paw 200 use interfaces_66_nonlocal 201 use interfaces_66_wfs 202 use interfaces_67_common, except_this_one => energy 203 !End of the abilint section 204 205 implicit none 206 207 !Arguments ------------------------------------ 208 !scalars 209 integer,intent(in) :: mcg,my_natom,n3xccc,nfftf,nhatgrdim,optene,usexcnhat 210 logical,intent(in),optional :: add_tfw 211 real(dp),intent(in) :: gsqcut 212 real(dp),intent(out) :: compch_fft,etotal 213 type(MPI_type),intent(inout) :: mpi_enreg 214 type(dataset_type),intent(in) :: dtset 215 type(electronpositron_type),pointer :: electronpositron 216 type(energies_type),intent(inout) :: energies 217 type(paw_dmft_type), intent(inout) :: paw_dmft 218 type(pawang_type),intent(in) :: pawang 219 type(pawfgr_type),intent(in) :: pawfgr 220 type(pseudopotential_type),intent(in) :: psps 221 type(wvl_internal_type), intent(in) :: wvl 222 type(wvl_wf_type),intent(inout) :: wfs 223 type(wvl_denspot_type), intent(inout) :: wvl_den 224 type(wvl_energy_terms),intent(inout) ::wvl_e 225 !arrays 226 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise 227 integer, intent(in) :: indsym(4,dtset%nsym,dtset%natom) 228 integer :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)),kg(3,dtset%mpw*dtset%mkmem) 229 integer, intent(in) :: ngfftf(18),npwarr(dtset%nkpt),symrec(3,3,dtset%nsym) 230 real(dp), intent(in) :: cg(2,mcg),eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 231 real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 232 real(dp), intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw) 233 real(dp),intent(in) :: nhatgr(nfftf,dtset%nspden,3*nhatgrdim) 234 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise 235 real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 236 real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 237 real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden) 238 real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden) 239 real(dp), intent(out) :: strsxc(6) 240 real(dp), intent(in) :: rprimd(3,3),vpsp(nfftf),xccc3d(n3xccc),xred(3,dtset%natom) 241 real(dp), intent(out) :: vhartr(nfftf),vtrial(nfftf,dtset%nspden),vxc(nfftf,dtset%nspden) 242 real(dp),intent(out) :: vxctau(nfftf,dtset%nspden,4*dtset%usekden) 243 real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 244 type(paw_ij_type), intent(in) :: paw_ij(my_natom*psps%usepaw) 245 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 246 type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw) 247 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 248 249 !Local variables------------------------------- 250 !scalars 251 integer :: bdtot_index,blocksize,choice,cplex,cplex_rf,cpopt,dimffnl 252 integer :: iband,iband_last,iblock,iblocksize,icg,ider,idir,ierr,ifft,ikg,ikpt,ilm 253 integer :: ipert,ipositron,iresid,ispden,isppol,istwf_k,izero 254 integer :: me_distrb,mpi_comm_sphgrid,my_ikpt,my_nspinor,n1,n2,n3,n4,n5,n6 255 integer :: nband_k,nblockbd,nfftotf,nkpg,nkxc,nk3xc,nnlout,npw_k,nspden_rhoij,option 256 integer :: option_rhoij,paw_opt,signs,spaceComm,tim_mkrho,tim_nonlop 257 logical :: add_tfw_,paral_atom,usetimerev,with_vxctau 258 logical :: wvlbigdft=.false. 259 logical :: non_magnetic_xc 260 real(dp) :: dotr,doti,eeigk,ekk,enlk,evxc,e_xcdc_vxctau,ucvol,ucvol_local,vxcavg 261 !character(len=500) :: message 262 type(gs_hamiltonian_type) :: gs_hamk 263 type(xcdata_type) :: xcdata 264 !arrays 265 integer,allocatable :: kg_k(:,:) 266 real(dp) :: gmet(3,3),gprimd(3,3),kpg_dum(0,0),kpoint(3),nonlop_out(1,1) 267 real(dp) :: qpt(3),rhodum(1),rmet(3,3),tsec(2),ylmgr_dum(1,1,1) 268 real(dp) :: vzeeman(4) 269 real(dp),allocatable :: buffer(:),cgrvtrial(:,:) 270 real(dp),allocatable :: cwavef(:,:),eig_k(:),enlout(:),ffnl(:,:,:,:),ffnl_sav(:,:,:,:) 271 real(dp),allocatable :: kinpw(:),kinpw_sav(:),kxc(:,:),occ_k(:),occblock(:) 272 real(dp),allocatable :: ph3d(:,:,:),ph3d_sav(:,:,:) 273 real(dp),allocatable :: resid_k(:),rhowfg(:,:),rhowfr(:,:),vlocal(:,:,:,:) 274 real(dp),allocatable :: vlocal_tmp(:,:,:),vxctaulocal(:,:,:,:,:),ylm_k(:,:),Vmagconstr(:,:) 275 type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null() 276 type(pawcprj_type),target,allocatable :: cwaveprj(:,:) 277 type(pawcprj_type),pointer :: cwaveprj_gat(:,:) 278 type(pawrhoij_type),pointer :: pawrhoij_unsym(:) 279 280 ! ************************************************************************* 281 282 DBG_ENTER("COLL") 283 284 ! Create variable for non_magnetic_xc 285 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 286 287 !Check that usekden is not 0 if want to use vxctau 288 with_vxctau = (dtset%usekden/=0) 289 290 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW) 291 nfftotf=PRODUCT(ngfftf(1:3)) 292 if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then 293 MSG_BUG('wrong values for nfft, nfftf!') 294 end if 295 296 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 297 wvlbigdft=(dtset%usewvl==1 .and. dtset%wvl_bigdft_comp==1) 298 299 call timab(59,1,tsec) 300 301 !Data for parallelism 302 spaceComm=mpi_enreg%comm_cell 303 if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 304 if(mpi_enreg%paral_hf==1) spaceComm=mpi_enreg%comm_kpt 305 mpi_comm_sphgrid=mpi_enreg%comm_fft 306 if(dtset%usewvl==1) then 307 spaceComm=mpi_enreg%comm_wvl 308 mpi_comm_sphgrid=mpi_enreg%comm_wvl 309 end if 310 me_distrb=mpi_enreg%me_kpt 311 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 312 paral_atom=(my_natom/=dtset%natom) 313 314 !Compute gmet, gprimd and ucvol from rprimd 315 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 316 if (dtset%usewvl == 0) then 317 ucvol_local = ucvol 318 #if defined HAVE_BIGDFT 319 else 320 ! We need to tune the volume when wavelets are used because, not 321 ! all FFT points are used. 322 ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 323 ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(product(wvl_den%denspot%dpbox%ndims), dp) 324 #endif 325 end if 326 327 !Compute Hxc potential from density 328 option=1;nkxc=0 329 ipositron=electronpositron_calctype(electronpositron) 330 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 331 332 if (ipositron/=1) then 333 334 if (dtset%icoulomb == 0) then 335 ! Use the periodic solver to compute Hxc. 336 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr) 337 call xcdata_init(xcdata,dtset=dtset) 338 ABI_ALLOCATE(kxc,(1,nkxc)) 339 ! to be adjusted for the call to rhotoxc 340 nk3xc=1 341 if (ipositron==0) then 342 call rhotoxc(energies%e_xc,kxc, & 343 & mpi_enreg,nfftf,ngfftf,nhat,psps%usepaw,nhatgr,nhatgrdim, & 344 & nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,rprimd,strsxc, & 345 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taug=taug,taur=taur,vhartr=vhartr, & 346 & vxctau=vxctau,exc_vdw_out=energies%e_xc_vdw,add_tfw=add_tfw_) 347 else 348 call rhotoxc(energies%e_xc,kxc, & 349 & mpi_enreg,nfftf,ngfftf,nhat,psps%usepaw,nhatgr,nhatgrdim, & 350 & nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,rprimd,strsxc, & 351 & usexcnhat,vxc,vxcavg,xccc3d,xcdata, & 352 & electronpositron=electronpositron,taug=taug,taur=taur,vhartr=vhartr, & 353 & vxctau=vxctau,exc_vdw_out=energies%e_xc_vdw,add_tfw=add_tfw_) 354 end if 355 ABI_DEALLOCATE(kxc) 356 else if (dtset%usewvl == 0) then 357 ! Use the free boundary solver. 358 call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, & 359 & dtset%icoulomb, dtset%ixc, mpi_enreg, nfftf, & 360 & ngfftf,nhat,psps%usepaw,& 361 & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, & 362 & usexcnhat,psps%usepaw,dtset%usewvl,& 363 & vhartr, vxc, vxcavg, wvl,wvl_den,wvl_e,& 364 & xccc3d,dtset%xclevel,dtset%xc_denpos) 365 end if 366 else 367 energies%e_xc=zero 368 call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfftf,ngfftf,nhat,nkxc,dtset%nspden,n3xccc,& 369 & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos) 370 end if 371 if (ipositron/=0) then 372 call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,nfftf,nfftotf,1,1,electronpositron%vha_ep,& 373 & ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 374 vhartr=vhartr+electronpositron%vha_ep 375 end if 376 377 !Total local potential (for either spin channel) is 378 !Hartree + local psp + Vxc(spin), minus its mean 379 !(Note : this potential should agree with the input vtrial) 380 do ispden=1,min(dtset%nspden,2) 381 do ifft=1,nfftf 382 vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden) 383 end do 384 end do 385 if (dtset%nspden==4) then 386 do ifft=1,nfftf 387 vtrial(ifft,3:4)=vxc(ifft,3:4) 388 end do 389 end if 390 391 !Add the vzeeman pot in the trial pot 392 !Vzeeman might have to be allocated correctly --> to be checked 393 if (any(abs(dtset%zeemanfield(:))>tol8)) then 394 vzeeman(:) = zero 395 if(dtset%nspden==2)then 396 vzeeman(2) = -half*dtset%zeemanfield(3) ! For collinear ispden=2 is rho_up only 397 end if 398 if(dtset%nspden==4)then 399 vzeeman(1)=-half*dtset%zeemanfield(3) 400 vzeeman(2)= half*dtset%zeemanfield(3) 401 vzeeman(3)=-half*dtset%zeemanfield(1) 402 vzeeman(4)= half*dtset%zeemanfield(2) 403 end if 404 do ispden=1,dtset%nspden 405 do ifft=1,nfftf 406 vtrial(ifft,ispden)=vtrial(ifft,ispden)+vzeeman(ispden) 407 end do 408 end do 409 end if 410 411 !Compute the constrained potential for the magnetic moments 412 !NOTE: here in energy.F90 rhor and vtrial are given on nfftf grid 413 !the values coming from mag_constr may be different from those calculated 414 !calling mag_constr with nfft in setvtr and rhotov 415 if (dtset%magconon==1.or.dtset%magconon==2) then 416 ABI_ALLOCATE(Vmagconstr, (nfftf,dtset%nspden)) 417 Vmagconstr = zero 418 call mag_constr(dtset%natom, dtset%spinat, dtset%nspden, dtset%magconon, dtset%magcon_lambda, rprimd, & 419 & mpi_enreg, nfftf, dtset%ngfft, dtset%ntypat, dtset%ratsph, rhor, & 420 & dtset%typat, Vmagconstr, xred) 421 do ispden=1,dtset%nspden 422 do ifft=1,nfftf 423 vtrial(ifft,ispden)=vtrial(ifft,ispden)+Vmagconstr(ifft,ispden) 424 end do 425 end do 426 ABI_DEALLOCATE(Vmagconstr) 427 end if 428 429 !Compute Hartree energy - use up+down rhor 430 if (ipositron/=1) then 431 call dotprod_vn(1,rhor,energies%e_hartree ,doti,nfftf,nfftotf,1,1,vhartr,& 432 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 433 if (ipositron==0) energies%e_hartree=half*energies%e_hartree 434 if (ipositron==2) energies%e_hartree = half *(energies%e_hartree-electronpositron%e_hartree) 435 else 436 energies%e_hartree=zero 437 end if 438 439 !Compute local psp energy - use up+down rhor 440 if (optene/=1) then 441 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfftf,nfftotf,1,1,vpsp,& 442 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 443 end if 444 445 !Compute DC-xc energy - use up+down rhor 446 if (optene>0) then 447 if (ipositron/=1) then 448 if (psps%usepaw==0.or.usexcnhat/=0) then 449 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfftf,nfftotf,dtset%nspden,1,vxc,& 450 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 451 if (with_vxctau)then 452 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfftf,nfftotf,dtset%nspden,1,vxctau(:,:,1),& 453 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 454 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 455 end if 456 else 457 ABI_ALLOCATE(rhowfr,(nfftf,dtset%nspden)) 458 rhowfr=rhor-nhat 459 call dotprod_vn(1,rhowfr,energies%e_xcdc,doti,nfftf,nfftotf,dtset%nspden,1,vxc,& 460 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 461 ABI_DEALLOCATE(rhowfr) 462 end if 463 if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc 464 else 465 energies%e_xcdc=zero 466 end if 467 end if 468 469 energies%e_eigenvalues=zero 470 energies%e_kinetic=zero 471 energies%e_nonlocalpsp=zero 472 bdtot_index=0 473 icg=0 474 475 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 476 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 477 478 !============================================ 479 !==== Initialize most of the Hamiltonian ==== 480 !============================================ 481 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 482 !2) Perform the setup needed for the non-local factors: 483 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 484 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 485 486 call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,& 487 & dtset%natom,dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 488 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 489 & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,& 490 & nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda) 491 492 ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamk%nvloc)) 493 if (with_vxctau) then 494 ABI_ALLOCATE(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4)) 495 end if 496 497 !PAW: additional initializations 498 if (psps%usepaw==1) then 499 ABI_DATATYPE_ALLOCATE(cwaveprj,(dtset%natom,my_nspinor)) 500 call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj) 501 if (mpi_enreg%paral_spinor==1) then 502 ABI_DATATYPE_ALLOCATE(cwaveprj_gat,(dtset%natom,dtset%nspinor)) 503 call pawcprj_alloc(cwaveprj_gat,0,gs_hamk%dimcprj) 504 else 505 cwaveprj_gat => cwaveprj 506 end if 507 if (paral_atom) then 508 ABI_DATATYPE_ALLOCATE(pawrhoij_unsym,(dtset%natom)) 509 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 510 call pawrhoij_alloc(pawrhoij_unsym,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,& 511 & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1) 512 else 513 pawrhoij_unsym => pawrhoij 514 call pawrhoij_init_unpacked(pawrhoij_unsym) 515 end if 516 option_rhoij=1 517 usetimerev=(dtset%kptopt>0.and.dtset%kptopt<3) 518 else 519 ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0)) 520 end if 521 522 !LOOP OVER SPINS 523 do isppol=1,dtset%nsppol 524 ikg=0 525 526 ! Set up local potential vlocal with proper dimensioning, from vtrial 527 ! Also take into account the spin. 528 if(dtset%nspden/=4)then 529 if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then 530 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal,2) 531 if(with_vxctau) then 532 do ispden=1,4 533 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,& 534 & vxctau(:,:,ispden),vxctaulocal(:,:,:,:,ispden),2) 535 end do 536 end if 537 else 538 ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden)) 539 call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial) 540 call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal,2) 541 ABI_DEALLOCATE(cgrvtrial) 542 end if 543 else 544 ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6)) 545 if (psps%usepaw==0) then 546 do ispden=1,dtset%nspden 547 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal_tmp,2) 548 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 549 end do 550 else 551 ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden)) 552 call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial) 553 do ispden=1,dtset%nspden 554 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal_tmp,2) 555 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 556 end do 557 ABI_DEALLOCATE(cgrvtrial) 558 end if 559 ABI_DEALLOCATE(vlocal_tmp) 560 end if 561 562 ! Continue Hamlitonian initializaton 563 call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.) 564 if (with_vxctau) then 565 call load_spin_hamiltonian(gs_hamk,isppol,vxctaulocal=vxctaulocal) 566 end if 567 568 ! Loop over k points 569 do ikpt=1,dtset%nkpt 570 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 571 istwf_k=dtset%istwfk(ikpt) 572 npw_k=npwarr(ikpt) 573 574 ! Skip this k-point if not the proper processor 575 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 576 resid(1+bdtot_index : nband_k+bdtot_index) = zero 577 bdtot_index=bdtot_index+nband_k 578 cycle 579 end if 580 581 ! Parallelism over FFT and/or bands: define sizes and tabs 582 if (mpi_enreg%paral_kgb==1) then 583 my_ikpt=mpi_enreg%my_kpttab(ikpt) 584 nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp) 585 my_bandfft_kpt => bandfft_kpt(my_ikpt) 586 else 587 my_ikpt=ikpt 588 nblockbd=nband_k/mpi_enreg%bandpp 589 !if (nband_k/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1 590 end if 591 blocksize=nband_k/nblockbd 592 593 ABI_ALLOCATE(eig_k,(nband_k)) 594 ABI_ALLOCATE(occ_k,(nband_k)) 595 ABI_ALLOCATE(resid_k,(nband_k)) 596 ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor*blocksize)) 597 resid_k(:)=zero 598 kpoint(:)=dtset%kptns(:,ikpt) 599 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 600 eig_k(:)=eigen(1+bdtot_index:nband_k+bdtot_index) 601 if (minval(eig_k)>1.d100) eig_k=zero 602 cplex=2 ; if (istwf_k>1) cplex=1 603 eeigk=zero ; ekk=zero ; enlk=zero 604 605 ABI_ALLOCATE(kg_k,(3,npw_k)) 606 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 607 608 ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 609 if (psps%useylm==1) then 610 do ilm=1,psps%mpsang*psps%mpsang 611 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 612 end do 613 end if 614 615 ! Compute kinetic energy 616 ABI_ALLOCATE(kinpw,(npw_k)) 617 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0) 618 619 ! Compute kinetic energy of each band 620 do iblock=1,nblockbd 621 do iblocksize=1,blocksize 622 iband=(iblock-1)*blocksize+iblocksize 623 if (abs(occ_k(iband))>tol8) then 624 cwavef(1:2,1:npw_k*my_nspinor)= & 625 & cg(:,1+(iband-1)*npw_k*my_nspinor+icg:iband*npw_k*my_nspinor+icg) 626 call meanvalue_g(dotr,kinpw,0,istwf_k,mpi_enreg,npw_k,my_nspinor,cwavef,cwavef,0) 627 energies%e_kinetic=energies%e_kinetic+dtset%wtk(ikpt)*occ_k(iband)*dotr 628 end if 629 end do 630 end do 631 632 ! Compute nonlocal form factors ffnl at all (k+G): 633 ider=0;dimffnl=1;nkpg=0 634 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 635 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,& 636 & gmet,gprimd,ider,ider,psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,& 637 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 638 & npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,& 639 & psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 640 641 ! Load k-dependent part in the Hamiltonian datastructure 642 ! - Compute 3D phase factors 643 ! - Prepare various tabs in case of band-FFT parallelism 644 ! - Load k-dependent quantities in the Hamiltonian 645 ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk)) 646 call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,& 647 & kinpw_k=kinpw,kg_k=kg_k,ffnl_k=ffnl,ph3d_k=ph3d,& 648 & compute_ph3d=.true.,compute_gbound=(mpi_enreg%paral_kgb/=1)) 649 650 ! Load band-FFT tabs (transposed k-dependent arrays) 651 if (mpi_enreg%paral_kgb==1) then 652 call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kinpw=kinpw_sav) 653 call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg) 654 call load_k_hamiltonian(gs_hamk,npw_fft_k=my_bandfft_kpt%ndatarecv, & 655 & gbound_k =my_bandfft_kpt%gbound, & 656 & kinpw_k =my_bandfft_kpt%kinpw_gather, & 657 & kg_k =my_bandfft_kpt%kg_k_gather, & 658 & ffnl_k =my_bandfft_kpt%ffnl_gather, & 659 & ph3d_k =my_bandfft_kpt%ph3d_gather) 660 end if 661 662 ! Setup gemm_nonlop 663 if (gemm_nonlop_use_gemm) then 664 gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt 665 call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax, & 666 & gs_hamk%ntypat, gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, gs_hamk%ffnl_k,& 667 & gs_hamk%ph3d_k) 668 end if 669 670 #if defined HAVE_GPU_CUDA 671 if (dtset%use_gpu_cuda==1) then 672 call gpu_update_ffnl_ph3d(ph3d,size(ph3d),ffnl,size(ffnl)) 673 end if 674 #endif 675 676 ! Compute nonlocal psp energy (NCPP) or Rhoij (PAW) 677 ABI_ALLOCATE(enlout,(blocksize)) 678 ABI_ALLOCATE(occblock,(blocksize)) 679 do iblock=1,nblockbd 680 iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k) 681 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle 682 683 ! Select occupied bands 684 occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index) 685 if(abs(maxval(occblock))>=tol8 ) then 686 cwavef(:,1:npw_k*my_nspinor*blocksize)=& 687 & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg) 688 689 choice=1-gs_hamk%usepaw ; signs=1 ; idir=0 ; nnlout=blocksize 690 paw_opt=gs_hamk%usepaw;cpopt=gs_hamk%usepaw-1 691 692 if (mpi_enreg%paral_kgb/=1) then 693 tim_nonlop=3 694 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,(/zero/),mpi_enreg,blocksize,nnlout,& 695 & paw_opt,signs,nonlop_out,tim_nonlop,cwavef,cwavef) 696 else 697 tim_nonlop=14 698 call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,(/zero/),blocksize,& 699 & mpi_enreg,nnlout,paw_opt,signs,nonlop_out,tim_nonlop,cwavef,cwavef) 700 end if 701 702 do iblocksize=1,blocksize 703 iband=(iblock-1)*blocksize+iblocksize 704 energies%e_eigenvalues=energies%e_eigenvalues+dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband) 705 energies%e_nonlocalpsp=energies%e_nonlocalpsp+dtset%wtk(ikpt)*occ_k(iband)*enlout(iblocksize) 706 end do 707 708 ! PAW: accumulate rhoij 709 if (psps%usepaw==1) then 710 if (mpi_enreg%paral_spinor==1) then 711 call pawcprj_gather_spin(cwaveprj,cwaveprj_gat,dtset%natom,1,my_nspinor,dtset%nspinor,& 712 & mpi_enreg%comm_spinor,ierr) 713 call pawaccrhoij(gs_hamk%atindx,cplex,cwaveprj_gat,cwaveprj_gat,0,isppol,dtset%natom,dtset%natom,& 714 & dtset%nspinor,occ_k(iband),option_rhoij,pawrhoij_unsym,usetimerev,dtset%wtk(ikpt)) 715 else 716 call pawaccrhoij(gs_hamk%atindx,cplex,cwaveprj,cwaveprj,0,isppol,dtset%natom,dtset%natom,& 717 & dtset%nspinor,occ_k(iband),option_rhoij,pawrhoij_unsym,usetimerev,dtset%wtk(ikpt)) 718 end if 719 end if 720 721 ! End loop on bands 722 end if 723 end do 724 725 ! Compute residual of each band (for informative purposes) 726 call mkresi(cg,eig_k,gs_hamk,icg,ikpt,isppol,mcg,mpi_enreg,nband_k,dtset%prtvol,resid_k) 727 resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:) 728 729 ! Restore the bandfft tabs 730 if (mpi_enreg%paral_kgb==1) then 731 call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kinpw=kinpw_sav) 732 end if 733 734 ! Incremente indexes 735 bdtot_index=bdtot_index+nband_k 736 if (dtset%mkmem/=0) then 737 icg=icg+npw_k*my_nspinor*nband_k 738 ikg=ikg+npw_k 739 end if 740 741 #if defined HAVE_GPU_CUDA 742 if(dtset%use_gpu_cuda==1) then 743 call gpu_finalize_ffnl_ph3d() 744 end if 745 #endif 746 747 ABI_DEALLOCATE(eig_k) 748 ABI_DEALLOCATE(occ_k) 749 ABI_DEALLOCATE(resid_k) 750 ABI_DEALLOCATE(enlout) 751 ABI_DEALLOCATE(occblock) 752 ABI_DEALLOCATE(ffnl) 753 ABI_DEALLOCATE(kinpw) 754 ABI_DEALLOCATE(ph3d) 755 ABI_DEALLOCATE(cwavef) 756 ABI_DEALLOCATE(kg_k) 757 ABI_DEALLOCATE(ylm_k) 758 759 ! End loops on isppol and ikpt 760 end do 761 end do 762 763 call destroy_hamiltonian(gs_hamk) 764 765 if(xmpi_paral==1)then 766 ! Accumulate enl eeig and ek on all proc. 767 ABI_ALLOCATE(buffer,(3+dtset%mband*dtset%nkpt*dtset%nsppol)) 768 buffer(1)=energies%e_nonlocalpsp ; buffer(2)=energies%e_kinetic ; buffer(3)=energies%e_eigenvalues 769 do iresid=1,dtset%mband*dtset%nkpt*dtset%nsppol 770 buffer(iresid+3)=resid(iresid) 771 end do 772 call timab(48,1,tsec) 773 call xmpi_sum(buffer,spaceComm,ierr) 774 call timab(48,2,tsec) 775 energies%e_nonlocalpsp=buffer(1) ; energies%e_kinetic=buffer(2) ; energies%e_eigenvalues=buffer(3) 776 do iresid=1,dtset%mband*dtset%nkpt*dtset%nsppol 777 resid(iresid)=buffer(iresid+3) 778 end do 779 ABI_DEALLOCATE(buffer) 780 ! Accumulate rhoij_ 781 if (psps%usepaw==1) then 782 call pawrhoij_mpisum_unpacked(pawrhoij_unsym,spaceComm,comm2=mpi_enreg%comm_band) 783 end if 784 end if 785 786 !Compute total (free) energy 787 if (optene==0.or.optene==2) then 788 etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc + & 789 & energies%e_localpsp + energies%e_corepsp 790 if (psps%usepaw==0) etotal=etotal + energies%e_nonlocalpsp 791 if (psps%usepaw==1) etotal=etotal + energies%e_paw 792 else if (optene==1.or.optene==3) then 793 etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - & 794 & energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc 795 if (psps%usepaw==1) etotal=etotal + energies%e_pawdc 796 end if 797 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 798 if(dtset%occopt>=3 .and. dtset%occopt<=8) etotal=etotal-dtset%tsmear*energies%entropy 799 800 !Additional stuff for electron-positron 801 if (dtset%positron/=0) then 802 if (ipositron==0) then 803 energies%e_electronpositron =zero 804 energies%edc_electronpositron=zero 805 else 806 energies%e_electronpositron =electronpositron%e_hartree+electronpositron%e_xc 807 energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc 808 if (psps%usepaw==1) then 809 energies%e_electronpositron =energies%e_electronpositron +electronpositron%e_paw 810 energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc 811 end if 812 end if 813 if (optene==0.or.optene==2) electronpositron%e0=etotal 814 if (optene==1.or.optene==3) electronpositron%e0=etotal-energies%edc_electronpositron 815 etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron 816 end if 817 818 !Compute new charge density based on incoming wf 819 !Keep rhor and rhog intact for later use e.g. in stress. (? MT 08-12-2008: is that true now ?) 820 !=== Norm-conserving psps: simply compute rho from WFs 821 !paw_dmft%use_dmft=0 ! dmft not used here 822 !paw_dmft%use_sc_dmft=0 ! dmft not used here 823 if (psps%usepaw==0) then 824 tim_mkrho=3 825 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,& 826 & npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wfs) 827 if(dtset%usekden==1)then 828 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,& 829 & npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl_den,wfs,option=1) 830 end if 831 else 832 ! === PAW case: symmetrize rhoij and add compensation charge density 833 tim_mkrho=3;option=1;choice=1 834 call symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,0,dtset%natom,dtset%nsym,& 835 & dtset%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,& 836 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 837 call pawrhoij_free_unpacked(pawrhoij_unsym) 838 if (paral_atom) then 839 call pawrhoij_free(pawrhoij_unsym) 840 ABI_DATATYPE_DEALLOCATE(pawrhoij_unsym) 841 end if 842 ider=0;izero=0;cplex_rf=1;ipert=0;idir=0;qpt(:)=zero 843 844 call pawmknhat(compch_fft,cplex_rf,ider,idir,ipert,izero,gprimd,& 845 & my_natom,dtset%natom,nfftf,ngfftf,& 846 & 0,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,rhodum,nhat,pawrhoij,pawrhoij,& 847 & pawtab,qpt,rprimd,ucvol_local,dtset%usewvl,xred,& 848 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 849 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 850 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 851 852 ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden)) 853 ABI_ALLOCATE(rhowfg,(2,dtset%nfft)) 854 rhowfr(:,:)=zero 855 856 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,& 857 & npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol_local,wvl_den,wfs) 858 859 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,0,dtset%paral_kgb,pawfgr,rhowfg,rhodum,rhowfr,rhor) 860 ABI_DEALLOCATE(rhowfr) 861 ABI_DEALLOCATE(rhowfg) 862 rhor(:,:)=rhor(:,:)+nhat(:,:) 863 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 864 865 end if 866 867 MSG_COMMENT('New density rho(r) made from input wfs') 868 869 call timab(59,2,tsec) 870 871 ABI_DEALLOCATE(vlocal) 872 if (with_vxctau) then 873 ABI_DEALLOCATE(vxctaulocal) 874 end if 875 876 if (psps%usepaw==1) then 877 call pawcprj_free(cwaveprj) 878 ABI_DATATYPE_DEALLOCATE(cwaveprj) 879 if (mpi_enreg%paral_spinor==1) then 880 call pawcprj_free(cwaveprj_gat) 881 ABI_DATATYPE_DEALLOCATE(cwaveprj_gat) 882 else 883 nullify(cwaveprj_gat) 884 end if 885 end if 886 887 DBG_EXIT("COLL") 888 889 end subroutine energy