TABLE OF CONTENTS


ABINIT/energy [ Functions ]

[ Top ] [ 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.

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

221 subroutine energy(cg,compch_fft,dtset,electronpositron,&
222 & energies,eigen,etotal,gsqcut,indsym,irrzon,kg,mcg,mpi_enreg,my_natom,nfftf,ngfftf,nhat,&
223 & nhatgr,nhatgrdim,npwarr,n3xccc,occ,optene,paw_dmft,paw_ij,pawang,pawfgr,&
224 & pawfgrtab,pawrhoij,pawtab,phnons,ph1d,psps,resid,rhog,rhor,rprimd,strsxc,symrec,&
225 & taug,taur,usexcnhat,vhartr,vtrial,vpsp,vxc,vxctau,wfs,wvl,wvl_den,wvl_e,xccc3d,xred,ylm,&
226 & add_tfw) ! optional argument
227 
228 
229 !This section has been created automatically by the script Abilint (TD).
230 !Do not modify the following lines by hand.
231 #undef ABI_FUNC
232 #define ABI_FUNC 'energy'
233 !End of the abilint section
234 
235  implicit none
236 
237 !Arguments ------------------------------------
238 !scalars
239  integer,intent(in) :: mcg,my_natom,n3xccc,nfftf,nhatgrdim,optene,usexcnhat
240  logical,intent(in),optional :: add_tfw
241  real(dp),intent(in) :: gsqcut
242  real(dp),intent(out) :: compch_fft,etotal
243  type(MPI_type),intent(inout) :: mpi_enreg
244  type(dataset_type),intent(in) :: dtset
245  type(electronpositron_type),pointer :: electronpositron
246  type(energies_type),intent(inout) :: energies
247  type(paw_dmft_type), intent(inout) :: paw_dmft
248  type(pawang_type),intent(in) :: pawang
249  type(pawfgr_type),intent(in) :: pawfgr
250  type(pseudopotential_type),intent(in) :: psps
251  type(wvl_internal_type), intent(in) :: wvl
252  type(wvl_wf_type),intent(inout) :: wfs
253  type(wvl_denspot_type), intent(inout) :: wvl_den
254  type(wvl_energy_terms),intent(inout) ::wvl_e
255 !arrays
256 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
257  integer, intent(in) :: indsym(4,dtset%nsym,dtset%natom)
258  integer :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)),kg(3,dtset%mpw*dtset%mkmem)
259  integer, intent(in) :: ngfftf(18),npwarr(dtset%nkpt),symrec(3,3,dtset%nsym)
260  real(dp), intent(in) :: cg(2,mcg),eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
261  real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
262  real(dp), intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw)
263  real(dp),intent(in) :: nhatgr(nfftf,dtset%nspden,3*nhatgrdim)
264 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
265  real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
266  real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
267  real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden)
268  real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden)
269  real(dp), intent(out) :: strsxc(6)
270  real(dp), intent(in) :: rprimd(3,3),vpsp(nfftf),xccc3d(n3xccc),xred(3,dtset%natom)
271  real(dp), intent(out) :: vhartr(nfftf),vtrial(nfftf,dtset%nspden),vxc(nfftf,dtset%nspden)
272  real(dp),intent(out) :: vxctau(nfftf,dtset%nspden,4*dtset%usekden)
273  real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
274  type(paw_ij_type), intent(in) :: paw_ij(my_natom*psps%usepaw)
275  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
276  type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
277  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
278 
279 !Local variables-------------------------------
280 !scalars
281  integer :: bdtot_index,blocksize,choice,cplex,cplex_rf,cpopt,dimffnl
282  integer :: iband,iband_last,iblock,iblocksize,icg,ider,idir,ierr,ifft,ikg,ikpt,ilm
283  integer :: ipert,ipositron,iresid,ispden,isppol,istwf_k,izero
284  integer :: me_distrb,mpi_comm_sphgrid,my_ikpt,my_nspinor,n1,n2,n3,n4,n5,n6
285  integer :: nband_k,nblockbd,nfftotf,nkpg,nkxc,nk3xc,nnlout,npw_k,nspden_rhoij,option
286  integer :: option_rhoij,paw_opt,signs,spaceComm,tim_mkrho,tim_nonlop
287  logical :: add_tfw_,paral_atom,usetimerev,with_vxctau
288  logical :: wvlbigdft=.false.
289  logical :: non_magnetic_xc
290  real(dp) :: dotr,doti,eeigk,ekk,enlk,evxc,e_xcdc_vxctau,ucvol,ucvol_local,vxcavg
291  !character(len=500) :: message
292  type(gs_hamiltonian_type) :: gs_hamk
293  type(xcdata_type) :: xcdata
294 !arrays
295  integer,allocatable :: kg_k(:,:)
296  real(dp) :: gmet(3,3),gprimd(3,3),kpg_dum(0,0),kpoint(3),nonlop_out(1,1)
297  real(dp) :: qpt(3),rhodum(1),rmet(3,3),tsec(2),ylmgr_dum(1,1,1)
298  real(dp) :: vzeeman(4)
299  real(dp),allocatable :: buffer(:),cgrvtrial(:,:)
300  real(dp),allocatable :: cwavef(:,:),eig_k(:),enlout(:),ffnl(:,:,:,:),ffnl_sav(:,:,:,:)
301  real(dp),allocatable :: kinpw(:),kinpw_sav(:),kxc(:,:),occ_k(:),occblock(:)
302  real(dp),allocatable :: ph3d(:,:,:),ph3d_sav(:,:,:)
303  real(dp),allocatable :: resid_k(:),rhowfg(:,:),rhowfr(:,:),vlocal(:,:,:,:)
304  real(dp),allocatable :: vlocal_tmp(:,:,:),vxctaulocal(:,:,:,:,:),ylm_k(:,:),Vmagconstr(:,:)
305  type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
306  type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
307  type(pawcprj_type),pointer :: cwaveprj_gat(:,:)
308  type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
309 
310 ! *************************************************************************
311 
312  DBG_ENTER("COLL")
313 
314 ! Create variable for non_magnetic_xc
315  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
316 
317 !Check that usekden is not 0 if want to use vxctau
318  with_vxctau = (dtset%usekden/=0)
319 
320 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
321  nfftotf=PRODUCT(ngfftf(1:3))
322  if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
323    MSG_BUG('wrong values for nfft, nfftf!')
324  end if
325 
326 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
327  wvlbigdft=(dtset%usewvl==1 .and. dtset%wvl_bigdft_comp==1)
328 
329  call timab(59,1,tsec)
330 
331 !Data for parallelism
332  spaceComm=mpi_enreg%comm_cell
333  if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
334  if(mpi_enreg%paral_hf==1) spaceComm=mpi_enreg%comm_kpt
335  mpi_comm_sphgrid=mpi_enreg%comm_fft
336  if(dtset%usewvl==1) then
337    spaceComm=mpi_enreg%comm_wvl
338    mpi_comm_sphgrid=mpi_enreg%comm_wvl
339  end if
340  me_distrb=mpi_enreg%me_kpt
341  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
342  paral_atom=(my_natom/=dtset%natom)
343 
344 !Compute gmet, gprimd and ucvol from rprimd
345  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
346  if (dtset%usewvl == 0) then
347    ucvol_local = ucvol
348 #if defined HAVE_BIGDFT
349  else
350 !  We need to tune the volume when wavelets are used because, not
351 !  all FFT points are used.
352 !  ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
353    ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(product(wvl_den%denspot%dpbox%ndims), dp)
354 #endif
355  end if
356 
357 !Compute Hxc potential from density
358  option=1;nkxc=0
359  ipositron=electronpositron_calctype(electronpositron)
360  add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
361 
362  if (ipositron/=1) then
363 
364    if (dtset%icoulomb == 0) then
365 !    Use the periodic solver to compute Hxc.
366      call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr)
367      call xcdata_init(xcdata,dtset=dtset)
368      ABI_ALLOCATE(kxc,(1,nkxc))
369 !    to be adjusted for the call to rhotoxc
370      nk3xc=1
371      if (ipositron==0) then
372        call rhotoxc(energies%e_xc,kxc, &
373 &       mpi_enreg,nfftf,ngfftf,nhat,psps%usepaw,nhatgr,nhatgrdim, &
374 &       nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,rprimd,strsxc, &
375 &       usexcnhat,vxc,vxcavg,xccc3d,xcdata,taug=taug,taur=taur,vhartr=vhartr, &
376 &       vxctau=vxctau,exc_vdw_out=energies%e_xc_vdw,add_tfw=add_tfw_)
377      else
378        call rhotoxc(energies%e_xc,kxc, &
379 &       mpi_enreg,nfftf,ngfftf,nhat,psps%usepaw,nhatgr,nhatgrdim, &
380 &       nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,rprimd,strsxc, &
381 &       usexcnhat,vxc,vxcavg,xccc3d,xcdata, &
382 &       electronpositron=electronpositron,taug=taug,taur=taur,vhartr=vhartr, &
383 &       vxctau=vxctau,exc_vdw_out=energies%e_xc_vdw,add_tfw=add_tfw_)
384      end if
385      ABI_DEALLOCATE(kxc)
386    else if (dtset%usewvl == 0) then
387 !    Use the free boundary solver.
388      call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
389 &     dtset%icoulomb, dtset%ixc, mpi_enreg, nfftf, &
390 &     ngfftf,nhat,psps%usepaw,&
391 &     dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
392 &     usexcnhat,psps%usepaw,dtset%usewvl,&
393 &     vhartr, vxc, vxcavg, wvl,wvl_den,wvl_e,&
394 &     xccc3d,dtset%xclevel,dtset%xc_denpos)
395    end if
396  else
397    energies%e_xc=zero
398    call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfftf,ngfftf,nhat,nkxc,dtset%nspden,n3xccc,&
399 &   dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
400  end if
401  if (ipositron/=0) then
402    call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,nfftf,nfftotf,1,1,electronpositron%vha_ep,&
403 &   ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
404    vhartr=vhartr+electronpositron%vha_ep
405  end if
406 
407 !Total local potential (for either spin channel) is
408 !Hartree + local psp + Vxc(spin), minus its mean
409 !(Note : this potential should agree with the input vtrial)
410  do ispden=1,min(dtset%nspden,2)
411    do ifft=1,nfftf
412      vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
413    end do
414  end do
415  if (dtset%nspden==4) then
416    do ifft=1,nfftf
417      vtrial(ifft,3:4)=vxc(ifft,3:4)
418    end do
419  end if
420 
421 !Add the vzeeman pot in the trial pot
422 !Vzeeman might have to be allocated correctly --> to be checked
423  if (any(abs(dtset%zeemanfield(:))>tol8)) then
424    vzeeman(:) = zero
425    if(dtset%nspden==2)then
426      vzeeman(2) = -half*dtset%zeemanfield(3) ! For collinear ispden=2 is rho_up only
427    end if
428    if(dtset%nspden==4)then
429      vzeeman(1)=-half*dtset%zeemanfield(3)
430      vzeeman(2)= half*dtset%zeemanfield(3)
431      vzeeman(3)=-half*dtset%zeemanfield(1)
432      vzeeman(4)= half*dtset%zeemanfield(2)
433    end if
434    do ispden=1,dtset%nspden
435      do ifft=1,nfftf
436        vtrial(ifft,ispden)=vtrial(ifft,ispden)+vzeeman(ispden)
437      end do
438    end do
439  end if
440 
441 !Compute the constrained potential for the magnetic moments
442 !NOTE: here in energy.F90 rhor and vtrial are given on nfftf grid
443 !the values coming from mag_constr may be different from those calculated
444 !calling mag_constr with nfft in setvtr and rhotov
445  if (dtset%magconon==1.or.dtset%magconon==2) then
446    ABI_ALLOCATE(Vmagconstr, (nfftf,dtset%nspden))
447    Vmagconstr = zero
448    call mag_constr(dtset%natom, dtset%spinat, dtset%nspden, dtset%magconon, dtset%magcon_lambda, rprimd, &
449 &   mpi_enreg, nfftf, dtset%ngfft, dtset%ntypat, dtset%ratsph, rhor, &
450 &   dtset%typat, Vmagconstr, xred)
451    do ispden=1,dtset%nspden
452      do ifft=1,nfftf
453        vtrial(ifft,ispden)=vtrial(ifft,ispden)+Vmagconstr(ifft,ispden)
454      end do
455    end do
456    ABI_DEALLOCATE(Vmagconstr)
457  end if
458 
459 !Compute Hartree energy - use up+down rhor
460  if (ipositron/=1) then
461    call dotprod_vn(1,rhor,energies%e_hartree ,doti,nfftf,nfftotf,1,1,vhartr,&
462 &   ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
463    if (ipositron==0) energies%e_hartree=half*energies%e_hartree
464    if (ipositron==2) energies%e_hartree = half *(energies%e_hartree-electronpositron%e_hartree)
465  else
466    energies%e_hartree=zero
467  end if
468 
469 !Compute local psp energy - use up+down rhor
470  if (optene/=1) then
471    call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfftf,nfftotf,1,1,vpsp,&
472 &   ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
473  end if
474 
475 !Compute DC-xc energy - use up+down rhor
476  if (optene>0) then
477    if (ipositron/=1) then
478      if (psps%usepaw==0.or.usexcnhat/=0) then
479        call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfftf,nfftotf,dtset%nspden,1,vxc,&
480 &       ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
481        if (with_vxctau)then
482          call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfftf,nfftotf,dtset%nspden,1,vxctau(:,:,1),&
483 &         ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
484          energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
485        end if
486      else
487        ABI_ALLOCATE(rhowfr,(nfftf,dtset%nspden))
488        rhowfr=rhor-nhat
489        call dotprod_vn(1,rhowfr,energies%e_xcdc,doti,nfftf,nfftotf,dtset%nspden,1,vxc,&
490 &       ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
491        ABI_DEALLOCATE(rhowfr)
492      end if
493      if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
494    else
495      energies%e_xcdc=zero
496    end if
497  end if
498 
499  energies%e_eigenvalues=zero
500  energies%e_kinetic=zero
501  energies%e_nonlocalpsp=zero
502  bdtot_index=0
503  icg=0
504 
505  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
506  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
507 
508 !============================================
509 !==== Initialize most of the Hamiltonian ====
510 !============================================
511 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
512 !2) Perform the setup needed for the non-local factors:
513 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
514 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
515 
516  call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,&
517 & dtset%natom,dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
518 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
519 & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,&
520 & nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
521 
522  ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamk%nvloc))
523  if (with_vxctau) then
524    ABI_ALLOCATE(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4))
525  end if
526 
527 !PAW: additional initializations
528  if (psps%usepaw==1) then
529    ABI_DATATYPE_ALLOCATE(cwaveprj,(dtset%natom,my_nspinor))
530    call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
531    if (mpi_enreg%paral_spinor==1) then
532      ABI_DATATYPE_ALLOCATE(cwaveprj_gat,(dtset%natom,dtset%nspinor))
533      call pawcprj_alloc(cwaveprj_gat,0,gs_hamk%dimcprj)
534    else
535      cwaveprj_gat => cwaveprj
536    end if
537    if (paral_atom) then
538      ABI_DATATYPE_ALLOCATE(pawrhoij_unsym,(dtset%natom))
539      nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
540      call pawrhoij_alloc(pawrhoij_unsym,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,&
541 &     dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
542    else
543      pawrhoij_unsym => pawrhoij
544      call pawrhoij_init_unpacked(pawrhoij_unsym)
545    end if
546    option_rhoij=1
547    usetimerev=(dtset%kptopt>0.and.dtset%kptopt<3)
548  else
549    ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0))
550  end if
551 
552 !LOOP OVER SPINS
553  do isppol=1,dtset%nsppol
554    ikg=0
555 
556 !  Set up local potential vlocal with proper dimensioning, from vtrial
557 !  Also take into account the spin.
558    if(dtset%nspden/=4)then
559      if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then
560        call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal,2)
561        if(with_vxctau) then
562          do ispden=1,4
563            call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,&
564 &           vxctau(:,:,ispden),vxctaulocal(:,:,:,:,ispden),2)
565          end do
566        end if
567      else
568        ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
569        call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
570        call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal,2)
571        ABI_DEALLOCATE(cgrvtrial)
572      end if
573    else
574      ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6))
575      if (psps%usepaw==0) then
576        do ispden=1,dtset%nspden
577          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal_tmp,2)
578          vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
579        end do
580      else
581        ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
582        call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
583        do ispden=1,dtset%nspden
584          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal_tmp,2)
585          vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
586        end do
587        ABI_DEALLOCATE(cgrvtrial)
588      end if
589      ABI_DEALLOCATE(vlocal_tmp)
590    end if
591 
592 !  Continue Hamlitonian initializaton
593    call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.)
594    if (with_vxctau) then
595      call load_spin_hamiltonian(gs_hamk,isppol,vxctaulocal=vxctaulocal)
596    end if
597 
598 !  Loop over k points
599    do ikpt=1,dtset%nkpt
600      nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
601      istwf_k=dtset%istwfk(ikpt)
602      npw_k=npwarr(ikpt)
603 
604 !    Skip this k-point if not the proper processor
605      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
606        resid(1+bdtot_index : nband_k+bdtot_index) = zero
607        bdtot_index=bdtot_index+nband_k
608        cycle
609      end if
610 
611 !    Parallelism over FFT and/or bands: define sizes and tabs
612      if (mpi_enreg%paral_kgb==1) then
613        my_ikpt=mpi_enreg%my_kpttab(ikpt)
614        nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
615        my_bandfft_kpt => bandfft_kpt(my_ikpt)
616      else
617        my_ikpt=ikpt
618        nblockbd=nband_k/mpi_enreg%bandpp
619        !if (nband_k/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1
620      end if
621      blocksize=nband_k/nblockbd
622 
623      ABI_ALLOCATE(eig_k,(nband_k))
624      ABI_ALLOCATE(occ_k,(nband_k))
625      ABI_ALLOCATE(resid_k,(nband_k))
626      ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor*blocksize))
627      resid_k(:)=zero
628      kpoint(:)=dtset%kptns(:,ikpt)
629      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
630      eig_k(:)=eigen(1+bdtot_index:nband_k+bdtot_index)
631      if (minval(eig_k)>1.d100) eig_k=zero
632      cplex=2 ; if (istwf_k>1) cplex=1
633      eeigk=zero ; ekk=zero ; enlk=zero
634 
635      ABI_ALLOCATE(kg_k,(3,npw_k))
636      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
637 
638      ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
639      if (psps%useylm==1) then
640        do ilm=1,psps%mpsang*psps%mpsang
641          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
642        end do
643      end if
644 
645 !    Compute kinetic energy
646      ABI_ALLOCATE(kinpw,(npw_k))
647      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
648 
649 !    Compute kinetic energy of each band
650      do iblock=1,nblockbd
651        do iblocksize=1,blocksize
652          iband=(iblock-1)*blocksize+iblocksize
653          if (abs(occ_k(iband))>tol8) then
654            cwavef(1:2,1:npw_k*my_nspinor)= &
655 &           cg(:,1+(iband-1)*npw_k*my_nspinor+icg:iband*npw_k*my_nspinor+icg)
656            call meanvalue_g(dotr,kinpw,0,istwf_k,mpi_enreg,npw_k,my_nspinor,cwavef,cwavef,0)
657            energies%e_kinetic=energies%e_kinetic+dtset%wtk(ikpt)*occ_k(iband)*dotr
658          end if
659        end do
660      end do
661 
662 !    Compute nonlocal form factors ffnl at all (k+G):
663      ider=0;dimffnl=1;nkpg=0
664      ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
665      call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
666 &     gmet,gprimd,ider,ider,psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,&
667 &     psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
668 &     npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
669 &     psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
670 
671 !    Load k-dependent part in the Hamiltonian datastructure
672 !     - Compute 3D phase factors
673 !     - Prepare various tabs in case of band-FFT parallelism
674 !     - Load k-dependent quantities in the Hamiltonian
675      ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
676      call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,&
677 &     kinpw_k=kinpw,kg_k=kg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
678 &     compute_ph3d=.true.,compute_gbound=(mpi_enreg%paral_kgb/=1))
679 
680 !    Load band-FFT tabs (transposed k-dependent arrays)
681      if (mpi_enreg%paral_kgb==1) then
682        call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kinpw=kinpw_sav)
683        call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg)
684        call load_k_hamiltonian(gs_hamk,npw_fft_k=my_bandfft_kpt%ndatarecv, &
685 &       gbound_k =my_bandfft_kpt%gbound, &
686 &       kinpw_k  =my_bandfft_kpt%kinpw_gather, &
687 &       kg_k     =my_bandfft_kpt%kg_k_gather, &
688 &       ffnl_k   =my_bandfft_kpt%ffnl_gather, &
689 &       ph3d_k   =my_bandfft_kpt%ph3d_gather)
690      end if
691 
692 !    Setup gemm_nonlop
693      if (gemm_nonlop_use_gemm) then
694        gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt
695        call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax, &
696 &       gs_hamk%ntypat, gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, gs_hamk%ffnl_k,&
697 &       gs_hamk%ph3d_k)
698      end if
699 
700 #if defined HAVE_GPU_CUDA
701      if (dtset%use_gpu_cuda==1) then
702        call gpu_update_ffnl_ph3d(ph3d,size(ph3d),ffnl,size(ffnl))
703      end if
704 #endif
705 
706 !    Compute nonlocal psp energy (NCPP) or Rhoij (PAW)
707      ABI_ALLOCATE(enlout,(blocksize))
708      ABI_ALLOCATE(occblock,(blocksize))
709      do iblock=1,nblockbd
710        iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
711        if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
712 
713 !      Select occupied bands
714        occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
715        if(abs(maxval(occblock))>=tol8 ) then
716          cwavef(:,1:npw_k*my_nspinor*blocksize)=&
717 &         cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
718 
719          choice=1-gs_hamk%usepaw ; signs=1 ; idir=0 ; nnlout=blocksize
720          paw_opt=gs_hamk%usepaw;cpopt=gs_hamk%usepaw-1
721 
722          if (mpi_enreg%paral_kgb/=1) then
723            tim_nonlop=3
724            call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,(/zero/),mpi_enreg,blocksize,nnlout,&
725 &           paw_opt,signs,nonlop_out,tim_nonlop,cwavef,cwavef)
726          else
727            tim_nonlop=14
728            call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,(/zero/),blocksize,&
729 &           mpi_enreg,nnlout,paw_opt,signs,nonlop_out,tim_nonlop,cwavef,cwavef)
730          end if
731 
732          do iblocksize=1,blocksize
733            iband=(iblock-1)*blocksize+iblocksize
734            energies%e_eigenvalues=energies%e_eigenvalues+dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband)
735            energies%e_nonlocalpsp=energies%e_nonlocalpsp+dtset%wtk(ikpt)*occ_k(iband)*enlout(iblocksize)
736          end do
737 
738 !        PAW: accumulate rhoij
739          if (psps%usepaw==1) then
740            if (mpi_enreg%paral_spinor==1) then
741              call pawcprj_gather_spin(cwaveprj,cwaveprj_gat,dtset%natom,1,my_nspinor,dtset%nspinor,&
742 &             mpi_enreg%comm_spinor,ierr)
743              call pawaccrhoij(gs_hamk%atindx,cplex,cwaveprj_gat,cwaveprj_gat,0,isppol,dtset%natom,dtset%natom,&
744 &             dtset%nspinor,occ_k(iband),option_rhoij,pawrhoij_unsym,usetimerev,dtset%wtk(ikpt))
745            else
746              call pawaccrhoij(gs_hamk%atindx,cplex,cwaveprj,cwaveprj,0,isppol,dtset%natom,dtset%natom,&
747 &             dtset%nspinor,occ_k(iband),option_rhoij,pawrhoij_unsym,usetimerev,dtset%wtk(ikpt))
748            end if
749          end if
750 
751 !        End loop on bands
752        end if
753      end do
754 
755 !    Compute residual of each band (for informative purposes)
756      call mkresi(cg,eig_k,gs_hamk,icg,ikpt,isppol,mcg,mpi_enreg,nband_k,dtset%prtvol,resid_k)
757      resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
758 
759 !    Restore the bandfft tabs
760      if (mpi_enreg%paral_kgb==1) then
761        call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kinpw=kinpw_sav)
762      end if
763 
764 !    Incremente indexes
765      bdtot_index=bdtot_index+nband_k
766      if (dtset%mkmem/=0) then
767        icg=icg+npw_k*my_nspinor*nband_k
768        ikg=ikg+npw_k
769      end if
770 
771 #if defined HAVE_GPU_CUDA
772      if(dtset%use_gpu_cuda==1) then
773        call gpu_finalize_ffnl_ph3d()
774      end if
775 #endif
776 
777      ABI_DEALLOCATE(eig_k)
778      ABI_DEALLOCATE(occ_k)
779      ABI_DEALLOCATE(resid_k)
780      ABI_DEALLOCATE(enlout)
781      ABI_DEALLOCATE(occblock)
782      ABI_DEALLOCATE(ffnl)
783      ABI_DEALLOCATE(kinpw)
784      ABI_DEALLOCATE(ph3d)
785      ABI_DEALLOCATE(cwavef)
786      ABI_DEALLOCATE(kg_k)
787      ABI_DEALLOCATE(ylm_k)
788 
789 !    End loops on isppol and ikpt
790    end do
791  end do
792 
793  call destroy_hamiltonian(gs_hamk)
794 
795  if(xmpi_paral==1)then
796 !  Accumulate enl eeig and ek on all proc.
797    ABI_ALLOCATE(buffer,(3+dtset%mband*dtset%nkpt*dtset%nsppol))
798    buffer(1)=energies%e_nonlocalpsp ; buffer(2)=energies%e_kinetic ; buffer(3)=energies%e_eigenvalues
799    do iresid=1,dtset%mband*dtset%nkpt*dtset%nsppol
800      buffer(iresid+3)=resid(iresid)
801    end do
802    call timab(48,1,tsec)
803    call xmpi_sum(buffer,spaceComm,ierr)
804    call timab(48,2,tsec)
805    energies%e_nonlocalpsp=buffer(1) ; energies%e_kinetic=buffer(2) ; energies%e_eigenvalues=buffer(3)
806    do iresid=1,dtset%mband*dtset%nkpt*dtset%nsppol
807      resid(iresid)=buffer(iresid+3)
808    end do
809    ABI_DEALLOCATE(buffer)
810 !  Accumulate rhoij_
811    if (psps%usepaw==1) then
812      call pawrhoij_mpisum_unpacked(pawrhoij_unsym,spaceComm,comm2=mpi_enreg%comm_band)
813    end if
814  end if
815 
816 !Compute total (free) energy
817  if (optene==0.or.optene==2) then
818    etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc + &
819 &   energies%e_localpsp + energies%e_corepsp
820    if (psps%usepaw==0) etotal=etotal + energies%e_nonlocalpsp
821    if (psps%usepaw==1) etotal=etotal + energies%e_paw
822  else if (optene==1.or.optene==3) then
823    etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - &
824 &   energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc
825    if (psps%usepaw==1) etotal=etotal + energies%e_pawdc
826  end if
827  etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
828  if(dtset%occopt>=3 .and. dtset%occopt<=8) etotal=etotal-dtset%tsmear*energies%entropy
829 
830 !Additional stuff for electron-positron
831  if (dtset%positron/=0) then
832    if (ipositron==0) then
833      energies%e_electronpositron  =zero
834      energies%edc_electronpositron=zero
835    else
836      energies%e_electronpositron  =electronpositron%e_hartree+electronpositron%e_xc
837      energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc
838      if (psps%usepaw==1) then
839        energies%e_electronpositron  =energies%e_electronpositron  +electronpositron%e_paw
840        energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc
841      end if
842    end if
843    if (optene==0.or.optene==2) electronpositron%e0=etotal
844    if (optene==1.or.optene==3) electronpositron%e0=etotal-energies%edc_electronpositron
845    etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron
846  end if
847 
848 !Compute new charge density based on incoming wf
849 !Keep rhor and rhog intact for later use e.g. in stress. (? MT 08-12-2008: is that true now ?)
850 !=== Norm-conserving psps: simply compute rho from WFs
851  !paw_dmft%use_dmft=0 ! dmft not used here
852  !paw_dmft%use_sc_dmft=0 ! dmft not used here
853  if (psps%usepaw==0) then
854    tim_mkrho=3
855    call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
856 &   npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wfs)
857    if(dtset%usekden==1)then
858      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
859 &     npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl_den,wfs,option=1)
860    end if
861  else
862 !  === PAW case: symmetrize rhoij and add compensation charge density
863    tim_mkrho=3;option=1;choice=1
864    call symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,0,dtset%natom,dtset%nsym,&
865 &   dtset%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,&
866 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
867    call pawrhoij_free_unpacked(pawrhoij_unsym)
868    if (paral_atom) then
869      call pawrhoij_free(pawrhoij_unsym)
870      ABI_DATATYPE_DEALLOCATE(pawrhoij_unsym)
871    end if
872    ider=0;izero=0;cplex_rf=1;ipert=0;idir=0;qpt(:)=zero
873 
874    call pawmknhat(compch_fft,cplex_rf,ider,idir,ipert,izero,gprimd,&
875 &   my_natom,dtset%natom,nfftf,ngfftf,&
876 &   0,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,rhodum,nhat,pawrhoij,pawrhoij,&
877 &   pawtab,qpt,rprimd,ucvol_local,dtset%usewvl,xred,&
878 &   comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
879 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
880 &   distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
881 
882    ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
883    ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
884    rhowfr(:,:)=zero
885 
886    call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
887 &   npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol_local,wvl_den,wfs)
888 
889    call transgrid(1,mpi_enreg,dtset%nspden,+1,1,0,dtset%paral_kgb,pawfgr,rhowfg,rhodum,rhowfr,rhor)
890    ABI_DEALLOCATE(rhowfr)
891    ABI_DEALLOCATE(rhowfg)
892    rhor(:,:)=rhor(:,:)+nhat(:,:)
893    call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
894 
895  end if
896 
897  MSG_COMMENT('New density rho(r) made from input wfs')
898 
899  call timab(59,2,tsec)
900 
901  ABI_DEALLOCATE(vlocal)
902  if (with_vxctau) then
903    ABI_DEALLOCATE(vxctaulocal)
904  end if
905 
906  if (psps%usepaw==1) then
907    call pawcprj_free(cwaveprj)
908    ABI_DATATYPE_DEALLOCATE(cwaveprj)
909    if (mpi_enreg%paral_spinor==1) then
910      call pawcprj_free(cwaveprj_gat)
911      ABI_DATATYPE_DEALLOCATE(cwaveprj_gat)
912    else
913      nullify(cwaveprj_gat)
914    end if
915  end if
916 
917  DBG_EXIT("COLL")
918 
919 end subroutine energy

ABINIT/m_dft_energy [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dft_energy

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_dft_energy
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_abicore
34  use m_hamiltonian
35  use m_errors
36  use m_xmpi
37  use m_gemm_nonlop
38  use m_xcdata
39  use m_cgtools
40 
41  use m_time,             only : timab
42  use m_geometry,         only : metric
43  use m_kg,               only : mkkin
44  use m_energies,         only : energies_type
45  use m_electronpositron, only : electronpositron_type, electronpositron_calctype, rhohxcpositron
46  use m_bandfft_kpt,      only : bandfft_kpt, bandfft_kpt_type, prep_bandfft_tabs, &
47                                 bandfft_kpt_savetabs, bandfft_kpt_restoretabs
48  use m_pawang,           only : pawang_type
49  use m_pawtab,           only : pawtab_type
50  use m_paw_ij,           only : paw_ij_type
51  use m_pawfgrtab,        only : pawfgrtab_type
52  use m_pawrhoij,         only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_init_unpacked, &
53                                 pawrhoij_mpisum_unpacked, pawrhoij_free_unpacked, pawrhoij_get_nspden, symrhoij
54  use m_pawcprj,          only : pawcprj_type,pawcprj_alloc,pawcprj_free,pawcprj_gather_spin
55  use m_pawfgr,           only : pawfgr_type
56  use m_paw_dmft,         only : paw_dmft_type
57  use m_paw_nhat,         only : pawmknhat
58  use m_paw_occupancies,  only : pawaccrhoij
59  use m_fft,              only : fftpac, fourdp
60  use m_spacepar,         only : meanvalue_g, hartre
61  use m_dens,             only : mag_constr
62  use m_mkrho,            only : mkrho
63  use m_mkffnl,           only : mkffnl
64  use m_getghc,           only : getghc
65  use m_rhotoxc,          only : rhotoxc
66  use m_mpinfo,           only : proc_distrb_cycle
67  use m_nonlop,           only : nonlop
68  use m_fourier_interpol, only : transgrid
69  use m_prep_kgb,         only : prep_getghc, prep_nonlop
70  use m_psolver,          only : psolver_rhohxc
71 
72  implicit none
73 
74  private

ABINIT/mkresi [ Functions ]

[ Top ] [ Functions ]

NAME

 mkresi

FUNCTION

 Make residuals from knowledge of wf in G space and application of Hamiltonian.

INPUTS

  cg(2,mcg)=<G|Cnk>=Fourier coefficients of wavefunction
  gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
  icg=shift to be applied on the location of data in the array cg
  ikpt=index of k-point
  isppol=index of spin
  mcg=second dimension of the cg array
  mpi_enreg=information about MPI parallelization
  nband=number of bands involved in subspace matrix.
  npw=number of planewaves in basis sphere at this k point.
  prtvol=control print volume and debugging output
  usepaw= 0 for non paw calculation; =1 for paw calculation

OUTPUT

  eig_k(nband)$= \langle C_n \mid H \mid C_n \rangle $ for each band.
  resid_k(nband)=residual for each band
   $= \langle C_n \mid H H \mid C_n \rangle- \langle C_n \mid H \mid C_n \rangle^2 $.

PARENTS

      energy

CHILDREN

      dotprod_g,getghc,prep_getghc,sqnorm_g,timab

SOURCE

 955 subroutine mkresi(cg,eig_k,gs_hamk,icg,ikpt,isppol,mcg,mpi_enreg,nband,prtvol,resid_k)
 956 
 957 
 958 !This section has been created automatically by the script Abilint (TD).
 959 !Do not modify the following lines by hand.
 960 #undef ABI_FUNC
 961 #define ABI_FUNC 'mkresi'
 962 !End of the abilint section
 963 
 964  implicit none
 965 
 966 !Arguments ------------------------------------
 967 !scalars
 968  integer,intent(in) :: icg,ikpt,isppol,mcg,nband,prtvol
 969  type(MPI_type),intent(inout) :: mpi_enreg
 970  type(gs_hamiltonian_type),intent(inout) :: gs_hamk
 971 !arrays
 972  real(dp),intent(in) :: cg(2,mcg)
 973  real(dp),intent(out) :: eig_k(nband),resid_k(nband)
 974 
 975 !Local variables-------------------------------
 976 !scalars
 977  integer,parameter :: tim_getghc=3
 978  integer :: blocksize,cpopt,iband,iband_last,iblock,iblocksize,ipw,ipw_shift
 979  integer :: my_nspinor,nblockbd,npw_k
 980  real(dp) :: doti,dotr
 981 !arrays
 982  real(dp) :: tsec(2)
 983  real(dp),allocatable,target :: cwavef(:,:),ghc(:,:),gsc(:,:),gvnlc(:,:)
 984  real(dp), ABI_CONTIGUOUS pointer :: cwavef_ptr(:,:),ghc_ptr(:,:),gsc_ptr(:,:)
 985  type(pawcprj_type) :: cwaveprj(0,0)
 986 
 987 ! *************************************************************************
 988 
 989 !Keep track of total time spent in mkresi
 990  call timab(13,1,tsec)
 991 
 992 !Parallelism over FFT and/or bands: define sizes and tabs
 993  my_nspinor=max(1,gs_hamk%nspinor/mpi_enreg%nproc_spinor)
 994  if (mpi_enreg%paral_kgb==1) then
 995    nblockbd=nband/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
 996  else
 997    nblockbd=nband/mpi_enreg%nproc_fft
 998    if (nband/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1
 999  end if
1000  blocksize=nband/nblockbd
1001 
1002  npw_k=gs_hamk%npw_k
1003  ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor))
1004  ABI_ALLOCATE(ghc,(2,npw_k*my_nspinor))
1005  ABI_ALLOCATE(gvnlc,(2,npw_k*my_nspinor))
1006  if (gs_hamk%usepaw==1)  then
1007    ABI_ALLOCATE(gsc,(2,npw_k*my_nspinor))
1008  else
1009    ABI_ALLOCATE(gsc,(0,0))
1010  end if
1011 
1012 !Loop over (blocks of) bands
1013  do iblock=1,nblockbd
1014    iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband)
1015    if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,mpi_enreg%me_kpt)) cycle
1016 
1017 !  Load |Cn>
1018    ipw_shift=(iblock-1)*npw_k*my_nspinor*blocksize+icg
1019 !$OMP PARALLEL DO
1020    do ipw=1,npw_k*my_nspinor*blocksize
1021      cwavef(1,ipw)=cg(1,ipw+ipw_shift)
1022      cwavef(2,ipw)=cg(2,ipw+ipw_shift)
1023    end do
1024 
1025 !  Compute H|Cn>
1026    cpopt=-1
1027    if (mpi_enreg%paral_kgb==0) then
1028      call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamk,gvnlc,zero,mpi_enreg,1,&
1029 &     prtvol,gs_hamk%usepaw,tim_getghc,0)
1030    else
1031      call prep_getghc(cwavef,gs_hamk,gvnlc,ghc,gsc,zero,nband,mpi_enreg,&
1032 &     prtvol,gs_hamk%usepaw,cpopt,cwaveprj,&
1033 &     already_transposed=.false.)
1034    end if
1035 
1036 !  Compute the residual, <Cn|(H-<Cn|H|Cn>)**2|Cn>:
1037    do iblocksize=1,blocksize
1038      iband=(iblock-1)*blocksize+iblocksize
1039      ipw_shift=(iblocksize-1)*npw_k*my_nspinor
1040      cwavef_ptr => cwavef(:,1+ipw_shift:npw_k*my_nspinor+ipw_shift)
1041      ghc_ptr    => ghc   (:,1+ipw_shift:npw_k*my_nspinor+ipw_shift)
1042 
1043 !    First get eigenvalue <Cn|H|Cn>:
1044      call dotprod_g(dotr,doti,gs_hamk%istwf_k,npw_k*my_nspinor,1,cwavef_ptr,ghc_ptr,&
1045 &     mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1046      eig_k(iband)=dotr
1047 
1048 !    Next need <G|(H-S<Cn|H|Cn>)|Cn> (in ghc):
1049      if (gs_hamk%usepaw==0) then
1050 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef_ptr,ghc_ptr,eig_k,iband,npw_k,my_nspinor)
1051        do ipw=1,npw_k*my_nspinor
1052          ghc_ptr(1,ipw)=ghc_ptr(1,ipw)-eig_k(iband)*cwavef_ptr(1,ipw)
1053          ghc_ptr(2,ipw)=ghc_ptr(2,ipw)-eig_k(iband)*cwavef_ptr(2,ipw)
1054        end do
1055      else
1056        gsc_ptr => gsc(:,1+ipw_shift:npw_k*my_nspinor+ipw_shift)
1057 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gsc_ptr,ghc_ptr,eig_k,iband,npw_k,my_nspinor)
1058        do ipw=1,npw_k*my_nspinor
1059          ghc_ptr(1,ipw)=ghc_ptr(1,ipw)-eig_k(iband)*gsc_ptr(1,ipw)
1060          ghc_ptr(2,ipw)=ghc_ptr(2,ipw)-eig_k(iband)*gsc_ptr(2,ipw)
1061        end do
1062      end if
1063 
1064 !    Then simply square the result:
1065      call sqnorm_g(dotr,gs_hamk%istwf_k,npw_k*my_nspinor,ghc_ptr,&
1066 &     mpi_enreg%me_g0,mpi_enreg%comm_fft)
1067      resid_k(iband)=dotr
1068 
1069    end do ! iblocksize
1070 
1071  end do ! iblock
1072 
1073  ABI_DEALLOCATE(cwavef)
1074  ABI_DEALLOCATE(ghc)
1075  ABI_DEALLOCATE(gvnlc)
1076  ABI_DEALLOCATE(gsc)
1077 
1078  call timab(13,2,tsec)
1079 
1080 end subroutine mkresi