TABLE OF CONTENTS


ABINIT/ionion_realspace [ Functions ]

[ Top ] [ Functions ]

NAME

 ionion_realspace

FUNCTION

 Compute the ion/ion interaction energies and forces in real space
 case. Use ewald() instead if computations are done in reciprocal
 space since it also includes the correction for the shift done in
 potentials calculations and includes replica interactions.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  rmet(3,3)=metric tensor in real space (bohr^2)
  xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
  zion(ntypat)=charge on each type of atom (real number)

OUTPUT

  eew=final ion/ion energy in hartrees
  grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.

PARENTS

      setvtr

CHILDREN

      xred2xcart

SOURCE

 957 subroutine ionion_realSpace(dtset, eew, grewtn, rprimd, xred, zion)
 958 
 959 
 960 !This section has been created automatically by the script Abilint (TD).
 961 !Do not modify the following lines by hand.
 962 #undef ABI_FUNC
 963 #define ABI_FUNC 'ionion_realSpace'
 964 !End of the abilint section
 965 
 966  implicit none
 967 
 968 !Arguments ------------------------------------
 969 !scalars
 970  real(dp),intent(out) :: eew
 971  type(dataset_type),intent(in) :: dtset
 972 !arrays
 973  real(dp),intent(in) :: rprimd(3,3),zion(dtset%ntypat)
 974  real(dp),intent(in) :: xred(3,dtset%natom)
 975  real(dp),intent(out) :: grewtn(3,dtset%natom)
 976 
 977 !Local variables-------------------------------
 978 !scalars
 979  integer :: ia1,ia2,iatom,igeo
 980  real(dp) :: r
 981 !arrays
 982  real(dp),allocatable :: grew_cart(:,:),xcart(:,:)
 983 
 984 ! *************************************************************************
 985 
 986 !Store xcart for each atom
 987  ABI_ALLOCATE(xcart,(3, dtset%natom))
 988  call xred2xcart(dtset%natom, rprimd, xcart, xred)
 989 
 990 !Summing the interaction between ions.
 991  eew = 0._dp
 992  do ia1 = 1, dtset%natom, 1
 993    do ia2 = ia1 + 1, dtset%natom, 1
 994      r = sqrt((xcart(1, ia1) - xcart(1, ia2)) ** 2 + &
 995 &     (xcart(2, ia1) - xcart(2, ia2)) ** 2 + &
 996 &     (xcart(3, ia1) - xcart(3, ia2)) ** 2)
 997      eew = eew + zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / r
 998    end do
 999  end do
1000 
1001 !Allocate temporary array to store cartesian gradients.
1002  ABI_ALLOCATE(grew_cart,(3, dtset%natom))
1003 
1004 !Summing the forces for each atom
1005  do ia1 = 1, dtset%natom, 1
1006    grew_cart(:, ia1) = 0._dp
1007    do ia2 = 1, dtset%natom, 1
1008      if (ia1 /= ia2) then
1009        r = (xcart(1, ia1) - xcart(1, ia2)) ** 2 + &
1010 &       (xcart(2, ia1) - xcart(2, ia2)) ** 2 + &
1011 &       (xcart(3, ia1) - xcart(3, ia2)) ** 2
1012        do igeo = 1, 3, 1
1013          grew_cart(igeo, ia1) = grew_cart(igeo, ia1) - (xcart(igeo, ia1) - xcart(igeo, ia2)) * &
1014 &         zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / (r ** 1.5_dp)
1015        end do
1016      end if
1017    end do
1018  end do
1019 
1020  ABI_DEALLOCATE(xcart)
1021 
1022 !Transform cartesian gradients to reduced gradients.
1023  do iatom = 1, dtset%natom, 1
1024    do igeo = 1, 3, 1
1025      grewtn(igeo, iatom) = rprimd(1, igeo) * grew_cart(1, iatom) + &
1026 &     rprimd(2, igeo) * grew_cart(2, iatom) + &
1027 &     rprimd(3, igeo) * grew_cart(3, iatom)
1028    end do
1029  end do
1030  ABI_DEALLOCATE(grew_cart)
1031 
1032 end subroutine ionion_realSpace

ABINIT/ionion_surface [ Functions ]

[ Top ] [ Functions ]

NAME

 ionion_surface

FUNCTION

 Compute the ion/ion interaction energies and forces in real space
 case. Use ewald() instead if computations are done in reciprocal
 space since it also includes the correction for the shift done in
 potentials calculations and includes replica interactions.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  rmet(3,3)=metric tensor in real space (bohr^2)
  xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
  zion(ntypat)=charge on each type of atom (real number)

OUTPUT

  eew=final ion/ion energy in hartrees
  grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.

PARENTS

      setvtr

CHILDREN

      ionicenergyandforces,xred2xcart

SOURCE

1064 subroutine ionion_surface(dtset, eew, grewtn, me, nproc, rprimd, wvl, wvl_den, xred)
1065 
1066 #if defined HAVE_BIGDFT
1067  use BigDFT_API, only: IonicEnergyandForces
1068 #endif
1069 
1070 !This section has been created automatically by the script Abilint (TD).
1071 !Do not modify the following lines by hand.
1072 #undef ABI_FUNC
1073 #define ABI_FUNC 'ionion_surface'
1074 !End of the abilint section
1075 
1076  implicit none
1077 
1078 !Arguments ------------------------------------
1079 !scalars
1080  integer, intent(in) :: me, nproc
1081  real(dp),intent(out) :: eew
1082  type(dataset_type),intent(in) :: dtset
1083  type(wvl_internal_type), intent(in) :: wvl
1084  type(wvl_denspot_type), intent(inout) :: wvl_den
1085 !arrays
1086  real(dp),intent(in) :: rprimd(3,3)
1087  real(dp),intent(in) :: xred(3,dtset%natom)
1088  real(dp),intent(out) :: grewtn(3,dtset%natom)
1089 
1090 !Local variables-------------------------------
1091 !scalars
1092  integer :: dispersion, iatom, igeo
1093  real(dp) :: psoffset
1094 !arrays
1095  real(dp),allocatable :: xcart(:,:)
1096  real(dp),pointer :: grew_cart(:,:),fdisp(:,:)
1097 #if defined HAVE_BIGDFT
1098  real(dp) :: edisp
1099  real(dp) :: ewaldstr(6)
1100 #endif
1101 
1102 ! *************************************************************************
1103 
1104 !Store xcart for each atom
1105  ABI_ALLOCATE(xcart,(3, dtset%natom))
1106  call xred2xcart(dtset%natom, rprimd, xcart, xred)
1107 
1108  nullify(fdisp)
1109  nullify(grew_cart)
1110  dispersion = 0
1111  psoffset = 0._dp
1112 #if defined HAVE_BIGDFT
1113  call IonicEnergyandForces(me, nproc, wvl_den%denspot%dpbox,&
1114 & wvl%atoms, dtset%efield, xcart, &
1115 & eew, grew_cart, dispersion, edisp, fdisp,&
1116 & ewaldstr,wvl%Glr%d%n1,wvl%Glr%d%n2,wvl%Glr%d%n3,&
1117 & wvl_den%denspot%V_ext, wvl_den%denspot%pkernel,psoffset)
1118 
1119  if (associated(fdisp)) then
1120    ABI_DEALLOCATE(fdisp)
1121  end if
1122 #endif
1123 
1124  ABI_DEALLOCATE(xcart)
1125 
1126 !Transform cartesian gradients to reduced gradients.
1127  do iatom = 1, dtset%natom, 1
1128    do igeo = 1, 3, 1
1129      grewtn(igeo, iatom) = -rprimd(1, igeo) * grew_cart(1, iatom) - &
1130 &     rprimd(2, igeo) * grew_cart(2, iatom) - &
1131 &     rprimd(3, igeo) * grew_cart(3, iatom)
1132    end do
1133  end do
1134  if (associated(grew_cart)) then
1135    ABI_DEALLOCATE(grew_cart)
1136  end if
1137 
1138 #if !defined HAVE_BIGDFT
1139  if (.false.) write(std_out,*) me,nproc,wvl%h(1),wvl_den%symObj
1140 #endif
1141 
1142 end subroutine ionion_surface

ABINIT/m_setvtr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_setvtr

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (XG, GMR, FJ, MT, EB, SPr)
  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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_setvtr
27 
28  use defs_basis
29  use defs_datatypes
30  use defs_abitypes
31  use defs_wvltypes
32  use m_abicore
33  use m_errors
34  use m_abi2big
35  use m_xmpi
36  use m_xcdata
37 
38  use m_time,              only : timab
39  use m_geometry,          only : xred2xcart
40  use m_cgtools,           only : dotprod_vn
41  use m_ewald,             only : ewald
42  use m_energies,          only : energies_type
43  use m_electronpositron,  only : electronpositron_type, electronpositron_calctype, rhohxcpositron
44  use libxc_functionals,   only : libxc_functionals_is_hybrid
45  use m_pawrad,            only : pawrad_type
46  use m_pawtab,            only : pawtab_type
47  use m_jellium,           only : jellium
48  use m_spacepar,          only : hartre
49  use m_dens,              only : mag_constr
50  use m_vdw_dftd2,         only : vdw_dftd2
51  use m_vdw_dftd3,         only : vdw_dftd3
52  use m_atm2fft,           only : atm2fft
53  use m_rhotoxc,           only : rhotoxc
54  use m_mklocl,            only : mklocl
55  use m_xchybrid,          only : xchybrid_ncpp_cc
56  use m_mkcore,            only : mkcore, mkcore_alt
57  use m_psolver,           only : psolver_rhohxc
58  use m_wvl_psi,          only : wvl_psitohpsi
59  use m_mkcore_wvl,       only : mkcore_wvl
60 
61 #if defined HAVE_BIGDFT
62  use BigDFT_API, only: denspot_set_history
63 #endif
64 
65  implicit none
66 
67  private

m_setvtr/setvtr [ Functions ]

[ Top ] [ m_setvtr ] [ Functions ]

NAME

 setvtr

FUNCTION

 Set up the trial potential and some energy terms

INPUTS

  [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
  atindx1(dtset%natom)=index table for atoms, inverse of atindx
  dtset <type(dataset_type)>=all input variables in this dataset
   |       if =0,1 no xc kernel, =2 spin-averaged (LDA) kernel
   | densfor_pred=govern the choice of preconditioner for the SCF cycle
   | iscf=determines the way the SCF cycle is handled
   | natom=number of atoms in cell.
   | nspden=number of spin-density components
   | qprtrb(3)= integer wavevector of possible perturbing potential
   |            in basis of reciprocal lattice translations
   | typat(natom)=type integer for each atom in cell
   | vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
   |  perturbing potential is added of the form
   |  V(G)=(vprtrb(1)+I*vprtrb(2))/2 at the values G=qprtrb and
   |  (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb (integers)
   |  for each type of atom, from psp (used in norm-conserving only)
  gmet(3,3)=metric tensor for G vecs (in bohr**-2)
  gprimd(3,3)=dimensional reciprocal space primitive translations
  gsqcut=cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential)
  istep=step number in the main loop of scfcv
  mgfft=maximum size of 1D FFTs
  moved_rhor=1 if the density was moved just before
  mpi_enreg=information about MPI parallelization
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat(nfft,nspden*usepaw)= -PAW only- compensation density
  nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
  nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
  nkxc=second dimension of the array kxc
  ntypat=number of types of atoms in unit cell.
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  n3xccc=dimension of the xccc3d array (0 or nfft).
  optene=>0 if some additional energies have to be computed
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=phase (structure factor) information.
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhog(2,nfft)=Fourier transform of electron density
  rhor(nfft,nspden)=electron density in electrons/bohr**3.
   | definition for spin components:
   | case of nspden = 2
   |      rhor(:,1) => rho_up + rho_dwn
   |      rhor(:,2) => rho_up
   | case of nspden = 4
   |      rhor(:,1)   => rho_upup + rho_dwndwn
   |      rhor(:,2:4) => {m_x,m_y,m_z}
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  [taug(2,nfftf*dtset%usekden)]=array for Fourier transform of kinetic energy density
  [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density
  ucvol = unit cell volume (bohr^3)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  energies <type(energies_type)>=all part of total energy.
   | e_xc=exchange-correlation energy (hartree)
   | In case of hybrid compensation algorithm:
   | 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 density
  ==== if optene==2 or 4
   | e_localpsp=local psp energy (hartree)
  ==== if dtset%icoulomb == 0
   | e_ewald=Ewald energy (hartree)
  ==== if optene>=1
   | e_hartree=Hartree part of total energy (hartree)
  ==== if optene==3 or 4
   | e_xcdc=exchange-correlation double-counting energy (hartree)
  ==== if dtset%vdw_xc == 5 or 6 or 7
   | e_vdw_dftd=Dispersion energy from DFT-D Van der Waals correction (hartree)
  grchempottn(3,natom)=grads of spatially-varying chemical energy (hartree)
  grewtn(3,natom)=grads of Ewald energy (hartree)
  grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree)
  kxc(nfft,nkxc)=exchange-correlation kernel, will be computed if nkxc/=0 .
                 see routine rhotoxc for a more complete description
  strsxc(6)=xc contribution to stress tensor (hartree/bohr^3)
  vxcavg=mean of the vxc potential

SIDE EFFECTS

  [electronpositron <type(electronpositron_type)>]=quantities for the electron-positron annihilation (optional argument)
  moved_atm_inside=1 if the atomic positions were moved inside the SCF loop.
  vhartr(nfft)=Hartree potential (Hartree)
  vpsp(nfft)=local psp (Hartree)
  vtrial(nfft,nspden)= trial potential (Hartree)
  vxc(nfft,nspden)= xc potential (Hartree)
  [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output
       i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density
  [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to
    kinetic energy density (metaGGA cases) (optional output)
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3

NOTES

  In case of PAW calculations:
    All computations are done on the fine FFT grid.
    All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
    All arrays (densities/potentials...) are computed on this fine FFT grid.
  ! Developpers have to be careful when introducing others arrays:
      they have to be stored on the fine FFT grid.
  In case of norm-conserving calculations the FFT grid is the usual FFT grid.

PARENTS

      bethe_salpeter,scfcv,screening,sigma

CHILDREN

      atm2fft,denspot_set_history,dotprod_vn,ewald,hartre,ionion_realspace
      ionion_surface,jellium,mag_constr,mkcore,mkcore_alt,mkcore_wvl,mklocl
      psolver_rhohxc,rhohxcpositron,rhotoxc,spatialchempot,timab,vdw_dftd2
      vdw_dftd3,wvl_psitohpsi,wvl_vtrial_abi2big,xcdata_init,xchybrid_ncpp_cc
      xred2xcart

SOURCE

198 subroutine setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,&
199 &  istep,kxc,mgfft,moved_atm_inside,moved_rhor,mpi_enreg,&
200 &  nattyp,nfft,ngfft,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,ntypat,n1xccc,n3xccc,&
201 &  optene,pawrad,pawtab,ph1d,psps,rhog,rhor,rmet,rprimd,strsxc,&
202 &  ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,xccc3d,xred,&
203 &  electronpositron,taug,taur,vxc_hybcomp,vxctau,add_tfw) ! optionals arguments
204 
205 
206 !This section has been created automatically by the script Abilint (TD).
207 !Do not modify the following lines by hand.
208 #undef ABI_FUNC
209 #define ABI_FUNC 'setvtr'
210 !End of the abilint section
211 
212  implicit none
213 
214 !Arguments ------------------------------------
215 !scalars
216  integer,intent(in) :: istep,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nhatgrdim,nkxc,ntypat
217  integer,intent(in) :: optene,usexcnhat
218  integer,intent(inout) :: moved_atm_inside,moved_rhor
219  logical,intent(in),optional :: add_tfw
220  real(dp),intent(in) :: gsqcut,ucvol
221  real(dp),intent(out) :: vxcavg
222  type(MPI_type),intent(in) :: mpi_enreg
223  type(dataset_type),intent(inout) :: dtset
224  type(electronpositron_type),pointer,optional :: electronpositron
225  type(energies_type),intent(inout) :: energies
226  type(pseudopotential_type),intent(in) :: psps
227  type(wvl_data), intent(inout) :: wvl
228 !arrays
229  integer, intent(in) :: atindx1(dtset%natom),nattyp(ntypat),ngfft(18)
230  real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
231  real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw)
232  real(dp),intent(in) :: nhatgr(:,:,:) !(nfft,dtset%nspden,3*nhatgrdim)
233  real(dp),intent(in) :: rhog(2,nfft)
234  real(dp),intent(in),optional :: taug(2,nfft*dtset%usekden)
235  real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
236  real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
237  real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft)
238  real(dp),intent(inout),optional :: taur(nfft,dtset%nspden*dtset%usekden)
239  real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden)
240  real(dp),intent(out),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden)
241  real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden)
242  real(dp),intent(inout) :: xccc3d(n3xccc)
243  real(dp),intent(in) :: xred(3,dtset%natom)
244  real(dp),intent(out) :: grchempottn(3,dtset%natom)
245  real(dp),intent(out) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc),strsxc(6)
246  type(pawtab_type),intent(in) :: pawtab(ntypat*dtset%usepaw)
247  type(pawrad_type),intent(in) :: pawrad(ntypat*dtset%usepaw)
248 
249 !Local variables-------------------------------
250 !scalars
251  integer :: coredens_method,mpi_comm_sphgrid,nk3xc
252  integer :: iatom,ifft,ipositron,ispden,nfftot
253  integer :: optatm,optdyfr,opteltfr,optgr,option,option_eff,optn,optn2,optstr,optv,vloc_method
254  real(dp) :: doti,e_xcdc_vxctau,ebb,ebn,evxc,ucvol_local,rpnrm
255  logical :: add_tfw_,is_hybrid_ncpp,non_magnetic_xc,with_vxctau,wvlbigdft
256  real(dp), allocatable :: xcart(:,:)
257  character(len=500) :: message
258  type(xcdata_type) :: xcdata,xcdatahyb
259 !arrays
260  real(dp),parameter :: identity(1:4)=(/1._dp,1._dp,0._dp,0._dp/)
261  real(dp) :: dummy6(6),tsec(2)
262  real(dp) :: grewtn_fake(3,1)
263  real(dp) :: dummy_in(0)
264  real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0)
265  real(dp) :: strn_dummy6(6), strv_dummy6(6)
266  real(dp) :: vzeeman(4)
267  real(dp),allocatable :: grtn(:,:),dyfr_dum(:,:,:),gr_dum(:,:)
268  real(dp),allocatable :: rhojellg(:,:),rhojellr(:),rhowk(:,:),vjell(:)
269  real(dp),allocatable :: Vmagconstr(:,:),rhog_dum(:,:)
270 
271 ! *********************************************************************
272 
273  call timab(91,1,tsec)
274 
275 ! Initialise non_magnetic_xc for rhohxc
276  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
277 
278 !Check that usekden is not 0 if want to use vxctau
279  with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0))
280 
281 !Check if we're in hybrid norm conserving pseudopotential with a core correction
282  is_hybrid_ncpp=(dtset%usepaw==0 .and. n3xccc/=0 .and. &
283 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
284 
285 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
286  wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
287 
288 !Get size of FFT grid
289  nfftot=PRODUCT(ngfft(1:3))
290 
291 !mpi
292  mpi_comm_sphgrid=mpi_enreg%comm_fft
293  if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
294 
295 !Test electron-positron case
296  ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron)
297 
298 !Test addition of Weiszacker gradient correction to Thomas-Fermi kin energy
299  add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
300 
301 !Get Ewald energy and Ewald forces, as well as vdW-DFTD energy and forces, and chemical potential energy and forces.
302 !-------------------------------------------------------------------------------------------------------------------
303  call timab(5,1,tsec)
304  if (ipositron/=1) then
305    if (dtset%icoulomb == 0 .or. (dtset%usewvl == 0 .and. dtset%icoulomb == 2)) then
306 !    Periodic system, need to compute energy and forces due to replica and
307 !    to correct the shift in potential calculation.
308      call ewald(energies%e_ewald,gmet,grewtn,dtset%natom,ntypat,rmet,dtset%typat,ucvol,xred,psps%ziontypat)
309 !    For a periodic system bearing a finite charge, the monopole correction to the
310 !    energy is relevant.
311 !    See Leslie and Gillan, JOURNAL OF PHYSICS C-SOLID STATE PHYSICS 18, 973 (1985)
312      if(abs(dtset%charge)>tol8) then
313        call ewald(energies%e_monopole,gmet,grewtn_fake,1,1,rmet,(/1/),ucvol,(/0.0_dp,0.0_dp,0.0_dp/),(/dtset%charge/))
314        energies%e_monopole=-energies%e_monopole
315      end if
316    else if (dtset%icoulomb == 1) then
317 !    In a non periodic system (real space computation), the G=0 divergence
318 !    doesn't occur and ewald is not needed. Only the ion/ion interaction
319 !    energy is relevant and used as Ewald energy and gradient.
320      call ionion_realSpace(dtset, energies%e_ewald, grewtn, rprimd, xred, psps%ziontypat)
321    else if (dtset%icoulomb == 2) then
322      call ionion_surface(dtset, energies%e_ewald, grewtn, mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, rprimd, &
323 &     wvl%descr, wvl%den, xred)
324    end if
325    if (dtset%nzchempot>0) then
326      call spatialchempot(energies%e_chempot,dtset%chempot,grchempottn,dtset%natom,ntypat,dtset%nzchempot,dtset%typat,xred)
327    end if
328    if (dtset%vdw_xc==5.and.ngrvdw==dtset%natom) then
329      call vdw_dftd2(energies%e_vdw_dftd,dtset%ixc,dtset%natom,ntypat,1,dtset%typat,rprimd,&
330 &     dtset%vdw_tol,xred,psps%znucltypat,fred_vdw_dftd2=grvdw)
331    end if
332    if ((dtset%vdw_xc==6.or.dtset%vdw_xc==7).and.ngrvdw==dtset%natom) then
333      call vdw_dftd3(energies%e_vdw_dftd,dtset%ixc,dtset%natom,&
334 &     ntypat,1,dtset%typat,rprimd,dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,&
335 &     xred,psps%znucltypat,fred_vdw_dftd3=grvdw)
336    end if
337  else
338    energies%e_ewald=zero
339    energies%e_chempot=zero
340    grchempottn=zero
341    grewtn=zero
342    energies%e_vdw_dftd=zero
343    if (ngrvdw>0) grvdw=zero
344  end if
345  call timab(5,2,tsec)
346 
347 !Compute parts of total energy depending on potentials
348 !--------------------------------------------------------------
349  if (dtset%usewvl == 0) then
350    ucvol_local = ucvol
351 #if defined HAVE_BIGDFT
352  else
353 !  We need to tune the volume when wavelets are used because, not all FFT points are used.
354 !  ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
355    ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp)
356 #endif
357  end if
358 
359 !Determine by which method the local ionic potential and/or the pseudo core charge density
360 ! have to be computed
361 !Local ionic potential:
362 ! Method 1: PAW
363 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
364  vloc_method=1;if (psps%usepaw==0) vloc_method=2
365  if (dtset%icoulomb>0) vloc_method=2
366  if (psps%usewvl==1) vloc_method=2
367 !Pseudo core charge density:
368 ! Method 1: PAW, nc_xccc_gspace
369 ! Method 2: Norm-conserving PP, wavelets
370  coredens_method=1;if (psps%usepaw==0) coredens_method=2
371  if (psps%nc_xccc_gspace==1) coredens_method=1
372  if (psps%nc_xccc_gspace==0) coredens_method=2
373  if (psps%usewvl==1) coredens_method=2
374 
375 !Local ionic potential and/or pseudo core charge by method 1
376  if (vloc_method==1.or.coredens_method==1) then
377    call timab(552,1,tsec)
378    optv=0;if (vloc_method==1) optv=1
379    optn=0;if (coredens_method==1) optn=n3xccc/nfft
380    optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
381    call atm2fft(atindx1,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,&
382 &   gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
383 &   optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
384 &   dummy_in,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
385 &   comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
386 &   paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
387    call timab(552,2,tsec)
388  end if
389 
390 !Local ionic potential by method 2
391  if (vloc_method==2) then
392    option=1
393    ABI_ALLOCATE(gr_dum,(3,dtset%natom))
394    ABI_ALLOCATE(dyfr_dum,(3,3,dtset%natom))
395    ABI_ALLOCATE(rhog_dum,(2,nfft))
396    call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,&
397 &   gr_dum,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,&
398 &   nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,&
399 &   dtset%qprtrb,rhog_dum,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
400    ABI_DEALLOCATE(gr_dum)
401    ABI_DEALLOCATE(dyfr_dum)
402    ABI_DEALLOCATE(rhog_dum)
403  end if
404 
405 !3D pseudo core electron density xccc3d by method 2
406  if (coredens_method==2.and.n1xccc/=0) then
407    call timab(91,2,tsec)
408    call timab(92,1,tsec)
409    option=1
410    ABI_ALLOCATE(gr_dum,(3,dtset%natom))
411    ABI_ALLOCATE(dyfr_dum,(3,3,dtset%natom))
412    if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then
413      call mkcore(dummy6,dyfr_dum,gr_dum,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
414 &     ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,&
415 &     vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred)
416    else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then
417      call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,&
418 &     nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
419 &     ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred,pawrad,pawtab,psps%usepaw)
420    else if (psps%usewvl==1.and.psps%usepaw==1) then
421 #if defined HAVE_BIGDFT
422 !      call mkcore_wvl_old(atindx1,dummy6,dyfr_dum,wvl%descr%atoms%astruct%geocode,gr_dum,wvl%descr%h,&
423 ! &         dtset%natom,nattyp,nfft,wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),&
424 ! &         dtset%nspden,ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,&
425 ! &         wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,wvl%den%denspot%dpbox%n3pi,n3xccc,option,&
426 ! &         pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d,xred,&
427 ! &         mpi_comm_wvl=mpi_enreg%comm_wvl)
428      call mkcore_wvl(atindx1,dummy6,gr_dum,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,&
429 &     n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d,&
430 &     psps%xcccrc,xred,wvl%den,wvl%descr,mpi_comm_wvl=mpi_enreg%comm_wvl)
431 #endif
432    end if
433    ABI_DEALLOCATE(gr_dum)
434    ABI_DEALLOCATE(dyfr_dum)
435    call timab(92,2,tsec)
436    call timab(91,1,tsec)
437  end if
438 
439 !Adds the jellium potential to the local part of ionic potential
440  if (dtset%jellslab/=0) then
441    ABI_ALLOCATE(vjell,(nfft))
442    ABI_ALLOCATE(rhojellg,(2,nfft))
443    ABI_ALLOCATE(rhojellr,(nfft))
444    option=1
445    call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,dtset%nspden,option,mpi_enreg%paral_kgb,&
446 &   dtset%slabwsrad,rhojellg,rhojellr,rprimd,vjell,dtset%slabzbeg,dtset%slabzend)
447 !  Compute background-background energy
448    call dotprod_vn(1,rhojellr,ebb,doti,nfft,nfftot,1,1,vjell,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
449    ebb=half*ebb
450 !  Compute electrostatic energy between background and nuclei before adding vjell to vpsp
451    call dotprod_vn(1,rhojellr,ebn,doti,nfft,nfftot,1,1,vpsp,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
452 !  Update e_ewald with ebb and ebn
453    energies%e_ewald=energies%e_ewald+ebb+ebn
454 !  Compute gradient of ebn wrt tn
455 !  This is not yet coded for usewvl or icoulomb=1
456    if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
457      write(message,'(3a)')&
458 &     'The computation of forces due to jellium background',ch10,&
459 &     'has to be verified in the PAW formalism.'
460      MSG_WARNING(message)
461 
462      ABI_ALLOCATE(grtn,(3,dtset%natom))
463      optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=1;optn=0;optn2=1
464      call atm2fft(atindx1,dummy_out1,vpsp,dummy_out2,dummy_out3,dummy_out4,dummy_in,&
465 &     gmet,gprimd,dummy_out5,grtn,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
466 &     optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
467 &     rhojellg,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
468 &     comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
469 &     paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
470 
471 !    Update grewtn with gradient of ebn wrt tn
472      do iatom=1,dtset%natom
473        grewtn(1:3,iatom)=grewtn(1:3,iatom)+grtn(1:3,iatom)
474      end do
475      ABI_DEALLOCATE(grtn)
476    else ! of usepaw==1
477      option=2
478      ABI_ALLOCATE(dyfr_dum,(3,3,dtset%natom))
479      ABI_ALLOCATE(grtn,(3,dtset%natom))
480      call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,&
481 &     grtn,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,&
482 &     nfft,ngfft,1,ntypat,option,pawtab,ph1d,psps,dtset%qprtrb,rhojellg,&
483 &     rhojellr,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
484 !    Update grewtn with gradient of ebn wrt tn (reestablish order of atoms)
485      do iatom=1,dtset%natom
486        grewtn(1:3,atindx1(iatom))=grewtn(1:3,atindx1(iatom))+grtn(1:3,iatom)
487      end do
488      ABI_DEALLOCATE(dyfr_dum)
489      ABI_DEALLOCATE(grtn)
490    end if ! of usepaw==1
491    vpsp(:)=vpsp(:)+vjell(:)
492    ABI_DEALLOCATE(vjell)
493    ABI_DEALLOCATE(rhojellg)
494    ABI_DEALLOCATE(rhojellr)
495  end if
496 
497 !Additional stuff for electron-positron calculation
498 !Compute the electronic/positronic local (Hartree) potential
499  if (ipositron==1) vpsp=-vpsp
500 
501 !If we are at the initialisation, or
502 !if the atom positions has changed and the non-linear core correction
503 !is included, or the rhor has changed, one needs to compute the xc stuff.
504 !One needs also to compute the Hartree stuff if the density changed,
505 !or at initialisation.
506 !--------------------------------------------------------------
507 
508  if(istep==1 .or. n1xccc/=0 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) then
509 
510    option=0
511    if(istep==1 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) option=1
512    if (nkxc>0) option=2
513    if (dtset%iscf==-1) option=-2
514    if (ipositron/=1) then
515      if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then
516        if(option/=0 .and. option/=10)then
517          call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartr)
518        end if
519        call xcdata_init(xcdata,dtset=dtset)
520        if(mod(dtset%fockoptmix,100)==11)then
521          xcdatahyb=xcdata
522 !        Setup the auxiliary xc functional information
523          call xcdata_init(xcdata,dtset=dtset,auxc_ixc=0,ixc=dtset%auxc_ixc)
524        end if
525 !      Use the periodic solver to compute Hxc
526        nk3xc=1
527        if (ipositron==0) then
528 
529 !        Compute energies%e_xc and associated quantities
530          if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then
531 !          Not yet able to deal fully with the full XC kernel in case of GGA + spin
532            option_eff=option
533            if(xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12
534            call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,&
535 &           nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
536 &           option_eff,dtset%paral_kgb,rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,&
537 &           taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
538          else
539 !          Only when is_hybrid_ncpp, and moreover, the xc functional is not the auxiliary xc functional, then call xchybrid_ncpp_cc
540            call xchybrid_ncpp_cc(dtset,energies%e_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
541 &           strsxc,vxcavg,xccc3d,vxc=vxc)
542          end if
543 
544 !        Possibly compute energies%e_hybcomp_E0
545          if(mod(dtset%fockoptmix,100)==11)then
546 !          This call to rhotoxc uses the hybrid xc functional
547            if(.not.is_hybrid_ncpp)then
548 !            Not yet able to deal fully with the full XC kernel in case of GGA + spin
549              option_eff=option
550              if(xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12
551              call rhotoxc(energies%e_hybcomp_E0,kxc,mpi_enreg,nfft,ngfft,&
552 &             nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
553 &             option,dtset%paral_kgb,rhor,rprimd,strsxc,usexcnhat,vxc_hybcomp,vxcavg,xccc3d,xcdatahyb,&
554 &             taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
555            else
556              call xchybrid_ncpp_cc(dtset,energies%e_hybcomp_E0,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
557 &             strsxc,vxcavg,xccc3d,vxc=vxc_hybcomp)
558            end if
559 
560 !          Combine hybrid and auxiliary quantities
561            energies%e_xc=energies%e_xc*dtset%auxc_scal
562            energies%e_hybcomp_E0=energies%e_hybcomp_E0-energies%e_xc
563            vxc(:,:)=vxc(:,:)*dtset%auxc_scal
564            vxc_hybcomp(:,:)=vxc_hybcomp(:,:)-vxc(:,:)
565 
566          end if
567 
568        else if (ipositron==2) then
569 !        Not yet able to deal fully with the full XC kernel in case of GGA + spin
570          option_eff=option
571          if(xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12
572          call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,&
573 &         nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
574 &         option,dtset%paral_kgb,rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,&
575 &         taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,&
576 &         electronpositron=electronpositron)
577        end if
578 
579      elseif(.not. wvlbigdft) then
580 !      Use the free boundary solver
581        call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
582 &       dtset%icoulomb, dtset%ixc, &
583 &       mpi_enreg, nfft, ngfft,&
584 &       nhat,psps%usepaw,&
585 &       dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
586 &       usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, &
587 &       vxcavg,wvl%descr,wvl%den,wvl%e,&
588 &       xccc3d,dtset%xclevel,dtset%xc_denpos)
589      end if
590    else
591      energies%e_xc=zero
592      call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,&
593 &     dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
594    end if
595    if (ipositron/=0) then
596      if (optene>=1) then
597        call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,&
598 &       nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
599      end if
600      vhartr=vhartr+electronpositron%vha_ep
601    end if
602  end if
603 
604 !Compute the trial potential
605 !-------------------------------------------------------------
606  if (.not. wvlbigdft) then
607 !  Now, compute trial Hxc potential. Local psp potential will be added later.
608    if(moved_atm_inside==0 .or.dtset%iscf>=10) then
609 
610 !    Compute starting Hxc potential.
611 !    Multiply by identity, should not change anything if nspden /= 4
612      do ispden=1,dtset%nspden
613        vtrial(:,ispden)=vhartr(:)*identity(ispden)+vxc(:,ispden)
614      end do
615 
616    else
617 
618 !    One should be here only when moved_atm_inside==1
619 !    The (H)xc now added corrects the previous one.
620      if(dtset%densfor_pred==1)then
621 !      xc was substracted off. This should be rationalized later
622        do ispden=1,dtset%nspden
623          vtrial(:,ispden)=vtrial(:,ispden)+vxc(:,ispden)
624        end do
625      else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then
626 !      Hxc was substracted off. This should be rationalized later
627        do ispden=1,dtset%nspden
628          vtrial(:,ispden)=vtrial(:,ispden)+vhartr(:)*identity(ispden)+vxc(:,ispden)
629        end do
630      end if
631    end if
632 
633 !  Adds the local part of the potential
634    if ((moved_atm_inside==0).or.(dtset%densfor_pred/=3)) then
635      do ispden=1,min(2,dtset%nspden)
636        do ifft=1,nfft
637          vtrial(ifft,ispden)=vtrial(ifft,ispden)+vpsp(ifft)
638        end do
639      end do
640    end if
641 
642 !  Adds the compensating vxc for hybrids
643    if(mod(dtset%fockoptmix,100)==11)then
644      vtrial(:,:)=vtrial(:,:)+vxc_hybcomp(:,:)
645    end if
646 
647    if(dtset%usewvl==1) then
648      call wvl_vtrial_abi2big(1,vtrial,wvl%den)
649    end if
650 
651  else
652 
653 !  Compute with covering comms the different part of the potential.
654 #if defined HAVE_BIGDFT
655 !  Copy e_ewald.
656    wvl%e%energs%eion = energies%e_ewald
657 !  Setup the mixing, if necessary
658    call denspot_set_history(wvl%den%denspot,dtset%iscf,dtset%nsppol, &
659 &   wvl%den%denspot%dpbox%ndims(1),wvl%den%denspot%dpbox%ndims(2))
660 #endif
661    ABI_ALLOCATE(xcart,(3, dtset%natom))
662    call xred2xcart(dtset%natom, rprimd, xcart, xred)
663    call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
664 &   energies%e_kinetic, energies%e_localpsp, energies%e_nonlocalpsp, energies%e_sicdc, &
665 &   istep, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl, dtset%nspden, &
666 &   rpnrm, .true.,evxc, wvl,.true., xcart, strsxc,&
667 &   vtrial, vxc)
668    if (optene==3.or.optene==4) energies%e_xcdc=evxc
669    ABI_DEALLOCATE(xcart)
670 
671  end if
672 
673 !Add the zeeman field to vtrial
674  if (any(abs(dtset%zeemanfield(:))>tol8)) then
675    vzeeman(:) = zero                            ! vzeeman_ij = -1/2*sigma_ij^alpha*B_alpha
676    if(dtset%nspden==2)then
677      vzeeman(1) = -half*dtset%zeemanfield(3)   ! v_dwndwn = -1/2*B_z
678      vzeeman(2) =  half*dtset%zeemanfield(3)   ! v_upup   =  1/2*B_z
679      do ifft=1,nfft
680        vtrial(ifft,1) = vtrial(ifft,1) + vzeeman(1) !SPr: added 1st component
681        vtrial(ifft,2) = vtrial(ifft,2) + vzeeman(2)
682      end do !ifft
683    end if
684    if(dtset%nspden==4)then
685      vzeeman(1)=-half*dtset%zeemanfield(3)     ! v_dwndwn                  => v_11
686      vzeeman(2)= half*dtset%zeemanfield(3)     ! v_upup                    => v_22
687      vzeeman(3)=-half*dtset%zeemanfield(1)     ! Re(v_dwnup) = Re(v_updwn) => Re(v_12)
688      vzeeman(4)= half*dtset%zeemanfield(2)     ! Im(v_dwnup) =-Im(v_dwnup) => Im(v_12)
689      do ispden=1,dtset%nspden
690        do ifft=1,nfft
691          vtrial(ifft,ispden) = vtrial(ifft,ispden) + vzeeman(ispden)
692        end do
693      end do
694    end if
695  end if
696 
697 !Compute the constrained potential for the magnetic moments
698  if (dtset%magconon==1.or.dtset%magconon==2) then
699    ABI_ALLOCATE(Vmagconstr, (nfft,dtset%nspden))
700    Vmagconstr = zero
701    call mag_constr(dtset%natom,dtset%spinat,dtset%nspden,dtset%magconon,dtset%magcon_lambda,rprimd, &
702 &   mpi_enreg,nfft,ngfft,dtset%ntypat,dtset%ratsph,rhor,dtset%typat,Vmagconstr,xred)
703    if(dtset%nspden==4)then
704      do ispden=1,dtset%nspden ! (SPr: both components should be used? EB: Yes it should be the case, corrected now)
705        do ifft=1,nfft
706          vtrial(ifft,ispden) = vtrial(ifft,ispden) + Vmagconstr(ifft,ispden)
707        end do !ifft
708      end do !ispden
709    else if(dtset%nspden==2)then
710      do ifft=1,nfft
711 !      TODO : MJV: check that magnetic constraint works also for nspden 2 or add input variable condition
712 !              EB: ispden=2 is rho_up only: to be tested
713 !             SPr: for ispden=2, both components should be used (e.g. see definition for vzeeman)?
714        vtrial(ifft,1) = vtrial(ifft,1) + Vmagconstr(ifft,1) !SPr: added the first component here
715        vtrial(ifft,2) = vtrial(ifft,2) + Vmagconstr(ifft,2)
716      end do !ifft
717    end if
718    ABI_DEALLOCATE(Vmagconstr)
719  end if
720 
721 !Compute parts of total energy depending on potentials
722 !--------------------------------------------------------------
723 
724 !For icoulomb==0 or usewvl Ehartree is calculated in psolver_rhohxc().
725 !For PAW we recalculate this since nhat was not taken into account
726 !in psolver_rhohxc: E_H= int v_H (n+nhat) dr
727 
728  if (optene>=1 .and. .not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then
729 !  Compute Hartree energy ehart
730 !  Already available in the Psolver case through psolver_rhohxc().
731    if (ipositron/=1) then
732      call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
733 &     mpi_comm_sphgrid=mpi_comm_sphgrid)
734      if (ipositron==0) energies%e_hartree = half * energies%e_hartree
735      if (ipositron==2) energies%e_hartree = half * (energies%e_hartree-electronpositron%e_hartree)
736    else
737      energies%e_hartree=zero
738    end if
739  end if
740 
741  if(mod(dtset%fockoptmix,100)==11)then
742    if (.not. wvlbigdft) then
743      call dotprod_vn(1,rhor,energies%e_hybcomp_v0,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol_local,&
744 &     mpi_comm_sphgrid=mpi_comm_sphgrid)
745      energies%e_hybcomp_v=energies%e_hybcomp_v0
746    end if
747  end if
748 
749  if (optene==2.or.optene==4 .and. .not. wvlbigdft) then
750 !  Compute local psp energy eei
751    call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,&
752 &   mpi_comm_sphgrid=mpi_comm_sphgrid)
753  end if
754 
755  if (optene==3.or.optene==4 .and. .not. wvlbigdft) then
756 !  Compute double-counting XC energy enxcdc
757    if (ipositron/=1) then
758      if (dtset%usepaw==0.or.usexcnhat/=0) then
759        call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
760 &       mpi_comm_sphgrid=mpi_comm_sphgrid)
761        if (with_vxctau) then
762          call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),&
763 &         ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
764          energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
765        end if
766      else
767        ABI_ALLOCATE(rhowk,(nfft,dtset%nspden))
768        rhowk=rhor-nhat
769        call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
770 &       mpi_comm_sphgrid=mpi_comm_sphgrid)
771        ABI_DEALLOCATE(rhowk)
772      end if
773      if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
774    else
775      energies%e_xcdc=zero
776    end if
777  end if
778 
779 !--------------------------------------------------------------
780 
781 !The initialisation for the new atomic positions has been done
782  moved_atm_inside=0
783 
784  call timab(91,2,tsec)
785 
786 end subroutine setvtr

m_setvtr/spatialchempot [ Modules ]

[ Top ] [ m_setvtr ] [ Modules ]

NAME

  spatialchempot

FUNCTION

  Treat spatially varying chemical potential.
  Compute energy and derivative with respect to dimensionless reduced atom coordinates of the
  spatially varying chemical potential. No contribution to stresses.

INPUTS

 chempot(3,nzchempot,ntypat)=input array with information about the chemical potential (see input variable description)
 natom=number of atoms in unit cell
 ntypat=number of type of atoms
 nzchempot=number of limiting planes for chemical potential
 typat(natom)=integer label of each type of atom (1,2,...)
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)

OUTPUT

 e_chempot=chemical potential energy in hartrees
 grchempottn(3,natom)=grads of e_chempot wrt xred(3,natom), hartrees.

PARENTS

      setvtr

CHILDREN

SOURCE

817 subroutine spatialchempot(e_chempot,chempot,grchempottn,natom,ntypat,nzchempot,typat,xred)
818 
819 
820 !This section has been created automatically by the script Abilint (TD).
821 !Do not modify the following lines by hand.
822 #undef ABI_FUNC
823 #define ABI_FUNC 'spatialchempot'
824 !End of the abilint section
825 
826  implicit none
827 
828 !Arguments ------------------------------------
829 !scalars
830  integer,intent(in) :: natom,ntypat,nzchempot
831  real(dp),intent(out) :: e_chempot
832 !arrays
833  integer,intent(in) :: typat(natom)
834  real(dp),intent(in) :: chempot(3,nzchempot,ntypat),xred(3,natom)
835  real(dp),intent(out) :: grchempottn(3,natom)
836 
837 !Local variables-------------------------------
838 !scalars
839  integer :: iatom,itypat,iz
840  real(dp) :: a_2,a_3,cp0,cp1,dcp0,dcp1,ddz,deltaz,deltaziz
841  real(dp) :: dqz,dz1,qz,zred,z0
842 !character(len=500) :: message
843 
844 ! *************************************************************************
845 
846 !DEBUG
847 !write(std_out,'(a)')' spatialchempot : enter '
848 !write(std_out,'(a,2i6)')' nzchempot,ntypat',nzchempot,ntypat
849 !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,1)=',chempot(1:3,1:2,1)
850 !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,2)=',chempot(1:3,1:2,2)
851 !ENDDEBUG
852 
853  e_chempot=zero
854  grchempottn(:,:)=zero
855 
856 !Loop on the different atoms
857  do iatom=1,natom
858 
859    itypat=typat(iatom)
860    zred=xred(3,iatom)
861 
862 !  Determine the delimiting plane just lower to zred
863 !  First compute the relative zred with respect to the first delimiting plane
864 !  Take into account a tolerance :
865    deltaz=zred-chempot(1,1,itypat)
866    deltaz=modulo(deltaz+tol12,1.0d0)-tol12
867 !  deltaz is positive (or higher than -tol12), and lower than one-tol12.
868    do iz=2,nzchempot+1
869      if(iz/=nzchempot+1)then
870        deltaziz=chempot(1,iz,itypat)-chempot(1,1,itypat)
871      else
872        deltaziz=one
873      end if
874      if(deltaziz>deltaz)exit
875    end do
876 
877 !  Defines coordinates and values inside the delimiting interval,
878 !  with respect to the lower delimiting plane
879    z0=chempot(1,iz-1,itypat)-chempot(1,1,itypat) ; cp0=chempot(2,iz-1,itypat) ; dcp0=chempot(3,iz-1,itypat)
880    if(iz/=nzchempot+1)then
881      dz1=chempot(1,iz,itypat)-chempot(1,iz-1,itypat) ; cp1=chempot(2,iz,itypat) ; dcp1=chempot(3,iz,itypat)
882    else
883      dz1=(chempot(1,1,itypat)+one)-chempot(1,nzchempot,itypat) ; cp1=chempot(2,1,itypat) ; dcp1=chempot(3,1,itypat)
884    end if
885    ddz=deltaz-z0
886 
887 !DEBUG
888 !  write(std_out,'(a,2i5)')' Delimiting planes, iz-1 and iz=',iz-1,iz
889 !  write(std_out,'(a,2es13.3)')' z0,  dz1= :',z0,dz1
890 !  write(std_out,'(a,2es13.3)')' cp0, cp1= :',cp0,cp1
891 !  write(std_out,'(a,2es13.3)')' dcp0, dcp1= :',dcp0,dcp1
892 !  write(std_out,'(a,2es13.3)')' deltaz,ddz=',deltaz,ddz
893 !ENDDEBUG
894 
895 !  Determine the coefficient of the third-order polynomial taking z0 as origin
896 !  P(dz=z-z0)= a_3*dz**3 + a_2*dz**2 + a_1*dz + a_0 ; obviously a_0=cp0 and a_1=dcp0
897 !  Define qz=a_3*dz + a_2 and dqz=3*a_3*dz + 2*a_2
898    qz=((cp1-cp0)-dcp0*dz1)/dz1**2
899    dqz=(dcp1-dcp0)/dz1
900    a_3=(dqz-two*qz)/dz1
901    a_2=three*qz-dqz
902 
903 !  Compute value and gradient of the chemical potential, at ddz wrt to lower delimiting plane
904    e_chempot=e_chempot+(a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0)
905    grchempottn(3,iatom)=three*a_3*ddz**2 + two*a_2*ddz + dcp0
906 
907 !DEBUG
908 !  write(std_out,'(a,4es16.6)')' qz,dqz=',qz,dqz
909 !  write(std_out,'(a,4es16.6)')' cp0,dcp0,a_2,a_3=',cp0,dcp0,a_2,a_3
910 !  write(std_out,'(a,2es13.3)')' dcp0*ddz + cp0=',dcp0*ddz + cp0
911 !  write(std_out,'(a,2es13.3)')' a_2*ddz**2=',a_2*ddz**2
912 !  write(std_out,'(a,2es13.3)')' a_3*ddz**3=',a_3*ddz**3
913 !  write(std_out,'(a,2es13.3)')' contrib=',a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0
914 !  write(std_out,'(a,2es13.3)')' e_chempot=',e_chempot
915 !  write(std_out,'(a,3es20.10)')' grchempottn=',grchempottn(:,iatom)
916 !ENDDEBUG
917 
918  end do
919 
920 !DEBUG
921 !write(std_out,'(a)')' spatialchempot : exit '
922 !write(std_out,'(a,es16.6)') ' e_chempot=',e_chempot
923 !ENDDEBUG
924 
925 end subroutine spatialchempot