TABLE OF CONTENTS


ABINIT/odamix [ Functions ]

[ Top ] [ Functions ]

NAME

 odamix

FUNCTION

 This routine is called to compute the total energy and various parts of it.
 The routine computes -if requested- the forces.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (FJ, MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.

INPUTS

  [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
  dtset <type(dataset_type)>=all input variables in this dataset
 berryopt  =  4/14: electric field is on -> add the contribution of the
                      -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j  terms to the total energy
   = 6/16, or 7/17: electric displacement field is on  -> add the contribution of the
                      Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j  terms to the total energy
   | berryopt  = 5: magnetic field is on -> add the contribution of the
   |                - \Omega B.M term to the total energy
   |          /= 5: magnetic field is off
   | bfield     = cartesian coordinates of the magnetic field in atomic units
   | dfield     = cartesian coordinates of the electric displacement field in atomic units (berryopt==6/7)
   | efield     = cartesian coordinates of the electric field in atomic units  (berryopt==4)
   | red_dfield = reduced the electric displacement field  (berryopt==16/17)
   | red_efieldbar = reduced the electric field (ebar)  (berryopt==14)
   | iatfix(3,natom)=1 for frozen atom along some direction, 0 for unfrozen
   | ionmov=governs the movement of atoms (see help file)
   | natom=number of atoms in cell.
   | nconeq=number of atomic constraint equations
   | nspden=number of spin-density components
   | nsym=number of symmetry elements in space group
   | occopt=option for occupancies
   | prtvol=integer controlling volume of printed output
   | tsmear=smearing energy or temperature (if metal)
   | wtatcon(3,natom,nconeq)=weights for atomic constraints
   | xclevel= XC functional level
  gprimd(3,3)=dimensional reciprocal space primitive translations
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  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
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  ntypat=number of types of atoms in unit cell.
  nvresid(nfft,nspden)=potential or density residual
  n3xccc=dimension of the xccc3d array (0 or nfft).
  optres=0 if residual array (nvresid) contains the potential residual
        =1 if residual array (nvresid) contains the density residual
  paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrad
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhog(2,nfft)=array for Fourier transform of electron density
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3
  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)
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vhartr(nfft)=array for holding Hartree potential
  vpsp(nfft)=array for holding local psp
  vxc(nfft,nspden)=array for holding XC potential
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  deltae=change in total energy
         between the previous and present SCF cycle
  etotal=total energy (hartree)

SIDE EFFECTS

 Input/Output:
  elast=previous value of the energy,
        needed to compute deltae, then updated.
  energies <type(energies_type)>=all part of total energy.
   | entropy(IN)=entropy due to the occupation number smearing (if metal)
   | e_localpsp(IN)=local psp energy (hartree)
   | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree)
   | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
   | e_ewald(IN)=Ewald energy (hartree)
   | e_vdw_dftd(IN)=VdW DFT-D energy
   | e_hartree(IN)=Hartree part of total energy (hartree units)
   | e_corepsp(IN)=psp core-core energy
   | e_kinetic(IN)=kinetic energy part of total energy.
   | e_nonlocalpsp(IN)=nonlocal pseudopotential part of total energy.
   | e_xc(IN)=exchange-correlation energy (hartree)
   | e_xcdc(IN)=exchange-correlation double-counting energy (hartree)
   | e_paw(IN)=PAW spherical part energy
   | e_pawdc(IN)=PAW spherical part double-counting energy
   | e_elecfield(OUT)=the term of the energy functional that depends explicitely
   |                  on the electric field:  enefield = -ucvol*E*P
   | e_magfield(OUT)=the term of the energy functional that depends explicitely
   |                  on the magnetic field:  enmagfield = -ucvol*B*M
   | e_entropy(OUT)=entropy energy due to the occupation number smearing (if metal)
   |                this value is %entropy * dtset%tsmear (hartree).
  kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
  [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to
      kinetic energy density (metaGGA cases) (optional output)
  ===== if psps%usepaw==1
   pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
    (gradients of rhoij for each atom with respect to atomic positions are computed here)

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

      scfcv

CHILDREN

      dotprod_vn,fourdp,hartre,metric,pawdenpot,pawmknhat,rhotoxc,timab
      xcdata_init,xmpi_sum

SOURCE

129 #if defined HAVE_CONFIG_H
130 #include "config.h"
131 #endif
132 
133 #include "abi_common.h"
134 
135 subroutine odamix(deltae,dtset,elast,energies,etotal,&
136 &          gprimd,gsqcut,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
137 &          nkxc,ntypat,nvresid,n3xccc,optres,paw_ij,&
138 &          paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,&
139 &          red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,usepaw,&
140 &          usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,&
141 &          taug,taur,vxctau,add_tfw) ! optional arguments
142 
143  use defs_basis
144  use defs_datatypes
145  use defs_abitypes
146  use m_profiling_abi
147  use m_errors
148  use m_xmpi
149  use m_xcdata
150 
151  use m_geometry,     only : metric
152  use m_cgtools,    only : dotprod_vn
153  use m_pawang, only : pawang_type
154  use m_pawrad, only : pawrad_type
155  use m_pawtab, only : pawtab_type
156  use m_paw_an, only : paw_an_type
157  use m_paw_ij, only : paw_ij_type
158  use m_pawfgrtab, only : pawfgrtab_type
159  use m_pawrhoij, only : pawrhoij_type
160  use m_energies, only : energies_type
161 
162 !This section has been created automatically by the script Abilint (TD).
163 !Do not modify the following lines by hand.
164 #undef ABI_FUNC
165 #define ABI_FUNC 'odamix'
166  use interfaces_18_timing
167  use interfaces_53_ffts
168  use interfaces_56_xc
169  use interfaces_65_paw
170 !End of the abilint section
171 
172  implicit none
173 
174 !Arguments ------------------------------------
175 !scalars
176  integer,intent(in) :: my_natom,n3xccc,nfft,nkxc,ntypat,optres
177  integer,intent(in) :: usepaw,usexcnhat
178  logical,intent(in),optional :: add_tfw
179  real(dp),intent(in) :: gsqcut,ucvol
180  real(dp),intent(inout) :: elast
181  real(dp),intent(out) :: deltae,etotal,vxcavg
182  type(MPI_type),intent(in) :: mpi_enreg
183  type(dataset_type),intent(in) :: dtset
184  type(energies_type),intent(inout) :: energies
185  type(pawang_type),intent(in) :: pawang
186  type(pseudopotential_type),intent(in) :: psps
187 !arrays
188  integer,intent(in) :: ngfft(18)
189  logical :: add_tfw_
190  real(dp),intent(in) :: gprimd(3,3)
191  real(dp),intent(in) :: red_ptot(3),rprimd(3,3),vpsp(nfft),xred(3,dtset%natom)
192  real(dp),intent(in),optional :: taug(2,nfft*dtset%usekden),taur(nfft,dtset%nspden*dtset%usekden)
193  real(dp),intent(inout) :: kxc(nfft,nkxc),nhat(nfft,dtset%nspden*usepaw)
194  real(dp),intent(inout) :: nvresid(nfft,dtset%nspden),rhog(2,nfft)
195  real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft)
196  real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden)
197  real(dp),intent(inout) :: xccc3d(n3xccc)
198  real(dp),intent(out) :: strsxc(6)
199  real(dp),intent(inout),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden)
200  type(paw_an_type),intent(inout) :: paw_an(my_natom)
201  type(paw_ij_type),intent(inout) :: paw_ij(my_natom)
202  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
203  type(pawrad_type),intent(in) :: pawrad(ntypat)
204  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
205  type(pawtab_type),intent(in) :: pawtab(ntypat)
206 
207 !Local variables-------------------------------
208 !scalars
209  integer :: cplex,iatom,ider,idir,ierr,ifft,ipert,irhoij,ispden,itypat,izero,iir,jjr,kkr
210  integer :: jrhoij,klmn,klmn1,kmix,nfftot,nhatgrdim,nselect,nzlmopt,nk3xc,option,optxc
211  logical :: with_vxctau
212  logical :: non_magnetic_xc
213  real(dp) :: alphaopt,compch_fft,compch_sph,doti,e1t10,e_ksnm1,e_xcdc_vxctau
214  real(dp) :: eenth,fp0,gammp1,ro_dlt,ucvol_local
215  character(len=500) :: message
216  type(xcdata_type) :: xcdata
217 !arrays
218  real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3)
219  real(dp) :: gmet(3,3),gprimdlc(3,3),qpt(3),rmet(3,3),tsec(2)
220  real(dp),allocatable :: nhatgr(:,:,:),rhoijtmp(:,:)
221 
222 ! *********************************************************************
223 
224 !DEBUG
225 !write(std_out,*)' odamix : enter'
226 !ENDDEBUG
227 
228  call timab(80,1,tsec)
229 
230 ! Initialise non_magnetic_xc for rhohxc
231  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
232 
233 !Check that usekden is not 0 if want to use vxctau
234  with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0))
235 
236 !To be adjusted for the call to rhotoxc
237  add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
238  nk3xc=1
239 
240 !faire un test sur optres=1, usewvl=0, nspden=1,nhatgrdim
241  if(optres/=1)then
242    write(message,'(a,i0,a)')' optres=',optres,', not allowed in oda => stop '
243    MSG_ERROR(message)
244  end if
245 
246  if(dtset%usewvl/=0)then
247    write(message,'(a,i0,a)')' usewvl=',dtset%usewvl,', not allowed in oda => stop '
248    MSG_ERROR(message)
249  end if
250 
251  if(dtset%nspden/=1)then
252    write(message,'(a,i0,a)')'  nspden=',dtset%nspden,', not allowed in oda => stop '
253    MSG_ERROR(message)
254  end if
255 
256  if (my_natom>0) then
257    if(paw_ij(1)%has_dijhat==0)then
258      message = ' dijhat variable must be allocated in odamix ! '
259      MSG_ERROR(message)
260    end if
261    if(paw_ij(1)%cplex==2)then
262      message = ' complex dij not allowed in odamix! '
263      MSG_ERROR(message)
264    end if
265  end if
266 
267 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
268 !!!!!!!!!! calculation of f'(0)= Eband_new-EH_old-E_xcdc_old-Ek_old-E_loc_old-E_nonloc_old
269 !!!!!!!!!! save previous energy E(rho_tild_n)
270 
271  fp0=energies%e_eigenvalues-energies%h0-two*energies%e_hartree-energies%e_xcdc
272  if (usepaw==1) then
273    do iatom=1,my_natom
274      itypat=pawrhoij(iatom)%itypat
275      do ispden=1,pawrhoij(iatom)%nspden
276        jrhoij=1
277        do irhoij=1,pawrhoij(iatom)%nrhoijsel
278          klmn=pawrhoij(iatom)%rhoijselect(irhoij)
279          ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn)
280          e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden))
281          jrhoij=jrhoij+pawrhoij(iatom)%cplex
282        end do
283        klmn1=1
284        do klmn=1,pawrhoij(iatom)%lmn2_size
285          ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,ispden)*pawtab(itypat)%dltij(klmn)
286          e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden))
287          klmn1=klmn1+pawrhoij(iatom)%cplex
288        end do
289      end do
290      if (paw_ij(iatom)%ndij>=2.and.pawrhoij(iatom)%nspden==1) then
291        jrhoij=1
292        do irhoij=1,pawrhoij(iatom)%nrhoijsel
293          klmn=pawrhoij(iatom)%rhoijselect(irhoij)
294          ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,1)*pawtab(itypat)%dltij(klmn)
295          e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2))
296          jrhoij=jrhoij+pawrhoij(iatom)%cplex
297        end do
298        klmn1=1
299        do klmn=1,pawrhoij(iatom)%lmn2_size
300          ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,1)*pawtab(itypat)%dltij(klmn)
301          e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2))
302          klmn1=klmn1+pawrhoij(iatom)%cplex
303        end do
304        e1t10=half*e1t10
305      end if
306    end do
307    if (mpi_enreg%nproc_atom>1) then
308      call xmpi_sum(e1t10,mpi_enreg%comm_atom,ierr)
309    end if
310    fp0=fp0-e1t10
311  end if
312  e_ksnm1=etotal
313 
314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
315 !!!!! Calculation of quantities that do not depend on rho_n+1
316 
317 !PAW: eventually recompute compensation density (and gradients)
318  nhatgrdim=0
319  if (usepaw==1) then
320    ider=-1;if (dtset%xclevel==2.or.usexcnhat==0) ider=0
321    if (dtset%xclevel==2.and.usexcnhat==1) ider=ider+2
322    if (ider>0) then
323      nhatgrdim=1
324      ABI_ALLOCATE(nhatgr,(nfft,dtset%nspden,3))
325    end if
326    if (ider>=0) then
327      ider=0;izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero
328      call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
329      nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,&
330 &     nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,&
331 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
332 &     comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
333 &     distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
334    end if
335  end if
336 
337 !------Compute Hartree and xc potentials----------------------------------
338  nfftot=PRODUCT(ngfft(1:3))
339 
340  call hartre(1,gsqcut,usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartr)
341 
342  call xcdata_init(xcdata,dtset=dtset) 
343 
344 !Compute xc potential (separate up and down if spin-polarized)
345  optxc=1
346  call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,&
347 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,dtset%paral_kgb,rhor,rprimd,strsxc,&
348 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
349 
350 !------Compute parts of total energy depending on potentials--------
351 
352  ucvol_local=ucvol
353 
354 !Compute Hartree energy energies%e_hartree
355  call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
356 & mpi_comm_sphgrid=mpi_enreg%comm_fft)
357  energies%e_hartree=half*energies%e_hartree
358 
359 
360 !Compute local psp energy energies%e_localpsp
361  call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,&
362 & mpi_comm_sphgrid=mpi_enreg%comm_fft)
363 
364 !Compute double-counting XC energy energies%e_xcdc
365  call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
366 & mpi_comm_sphgrid=mpi_enreg%comm_fft)
367  if (with_vxctau) then
368    call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,&
369 &   vxctau(:,:,1),ucvol_local,mpi_comm_sphgrid=mpi_enreg%comm_fft)
370    energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
371  end if
372 
373  if (usepaw/=0) then
374    nzlmopt=dtset%pawnzlm; option=2
375    do iatom=1,my_natom
376      itypat=paw_ij(iatom)%itypat
377      ABI_ALLOCATE(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size))
378      paw_ij(iatom)%has_dijhartree=1
379    end do
380    call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,0,dtset%ixc,my_natom,dtset%natom,dtset%nspden,ntypat,&
381 &   dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,&
382 &   pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,&
383 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
384    do iatom=1,my_natom
385      ABI_DEALLOCATE(paw_ij(iatom)%dijhartree)
386      paw_ij(iatom)%has_dijhartree=0
387    end do
388  end if
389 
390 
391 !When the finite-temperature VG broadening scheme is used,
392 !the total entropy contribution "tsmear*entropy" has a meaning,
393 !and gather the two last terms of Eq.8 of the VG paper
394 !Warning : might have to be changed for fixed moment calculations
395  if(dtset%occopt>=3 .and. dtset%occopt<=8) then
396    energies%e_entropy = - dtset%tsmear * energies%entropy
397  else
398    energies%e_entropy = zero
399  end if
400 !Turn it into an electric enthalpy,refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009), the missing volumce is added here
401  energies%e_elecfield = zero
402  if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then             !!HONG
403 
404    energies%e_elecfield = -dot_product(dtset%red_efieldbar,red_ptot)
405 
406    call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
407    eenth = zero
408    do iir=1,3
409      do jjr=1,3
410        eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)         !! HONG g^{-1})_ij ebar_i ebar_j
411      end do
412    end do
413    eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth
414    energies%e_elecfield = energies%e_elecfield + eenth
415 
416  end if
417 
418  energies%e_magfield = zero
419 !if (dtset%berryopt == 5) then
420 !emag = dot_product(mag_cart,dtset%bfield)
421 !energies%e_magfield = emag
422 !end if
423 
424 !HONG  Turn it into an internal enthalpy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009), but a little different: U=E_ks + (vol/8*pi) *  g^{-1})_ij ebar_i ebar_j
425  if (dtset%berryopt == 6 .or. dtset%berryopt == 16 )  then
426    energies%e_elecfield=zero
427    call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
428    do iir=1,3
429      do jjr=1,3
430        energies%e_elecfield = energies%e_elecfield + gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)     !! HONG g^{-1})_ij ebar_i ebar_j
431      end do
432    end do
433    energies%e_elecfield = ucvol_local/(8.0d0*pi)*energies%e_elecfield
434  end if
435 
436 !HONG  calculate internal energy and electric enthalpy for mixed BC case.
437  if ( dtset%berryopt == 17 ) then
438    energies%e_elecfield = zero
439    call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
440    A(:,:)=(4*pi/ucvol_local)*rmet(:,:)
441    A1(:,:)=A(:,:)
442    A_new(:,:)=A(:,:)
443    efield_new(:)=dtset%red_efield(:)
444    eenth = zero
445 
446    do kkr=1,3
447      if (dtset%jfielddir(kkr)==1) then    ! fixed ebar direction
448 
449 !      step 1 add -ebar*p
450        eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr)
451 
452 !      step 2  chang to e_new (change e to ebar)
453        efield_new(kkr)=dtset%red_efieldbar(kkr)
454 
455 !      step 3  chang matrix A to  A1
456 
457        do iir=1,3
458          do jjr=1,3
459            if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
460            if ((iir==kkr .and. jjr/=kkr) .or.  (iir/=kkr .and.  jjr==kkr)) &
461 &           A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
462            if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
463          end do
464        end do
465 
466        A(:,:)=A1(:,:)
467        A_new(:,:)=A1(:,:)
468      end if
469 
470    end do  ! end fo kkr
471 
472 
473    do iir=1,3
474      do jjr=1,3
475        eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
476      end do
477    end do
478 
479    energies%e_elecfield=energies%e_elecfield+eenth
480 
481  end if   ! berryopt==17
482 
483 
484 
485 
486 
487  etotal = energies%e_kinetic+ energies%e_hartree + energies%e_xc + &
488 & energies%e_localpsp + energies%e_corepsp + &
489 & energies%e_entropy + energies%e_elecfield + energies%e_magfield
490 !etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - &
491 !& energies%e_xcdc + energies%e_corepsp + &
492 !& energies%e_entropy + energies%e_elecfield
493  etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
494  if (usepaw==0) then
495    etotal = etotal + energies%e_nonlocalpsp
496  else
497    etotal = etotal + energies%e_paw
498  end if
499 
500 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
501 !!!!!!!!!!!!!! now, compute mixed densities
502 
503  gammp1=etotal-e_ksnm1-fp0
504  if (fp0>0.d0) then
505    write(std_out,*) "fp0 est positif"
506 !  stop
507  end if
508  write(std_out,*) "fp0 ",fp0
509  alphaopt=-fp0/two/gammp1
510 
511  if (alphaopt>one.or.alphaopt<0.d0) alphaopt=one
512  if (abs(energies%h0)<=tol10) alphaopt=one
513  write(std_out,*) " alphaopt",alphaopt
514 
515  energies%h0=(one-alphaopt)*energies%h0 + alphaopt*(energies%e_kinetic+energies%e_localpsp)
516  if (usepaw==0) energies%h0=energies%h0 + alphaopt*energies%e_nonlocalpsp
517 
518  rhor= rhor+(alphaopt-one)*nvresid
519  call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
520 
521  if (usepaw==1) then
522    do iatom=1,my_natom
523      ABI_ALLOCATE(rhoijtmp,(pawrhoij(iatom)%cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
524      rhoijtmp=zero
525      if (pawrhoij(iatom)%cplex==1) then
526        if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
527          do ispden=1,pawrhoij(iatom)%nspden
528            do irhoij=1,pawrhoij(iatom)%nrhoijsel
529              klmn=pawrhoij(iatom)%rhoijselect(irhoij)
530              rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden)
531            end do
532          end do
533        end if
534        do ispden=1,pawrhoij(iatom)%nspden
535          do kmix=1,pawrhoij(iatom)%lmnmix_sz
536            klmn=pawrhoij(iatom)%kpawmix(kmix)
537            rhoijtmp(klmn,ispden)=rhoijtmp(klmn,ispden)+(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn,ispden)
538          end do
539        end do
540        nselect=0
541        do klmn=1,pawrhoij(iatom)%lmn2_size
542          if (any(abs(rhoijtmp(klmn,:))>tol10)) then
543            nselect=nselect+1
544            pawrhoij(iatom)%rhoijselect(nselect)=klmn
545            do ispden=1,pawrhoij(iatom)%nspden
546              pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden)
547            end do
548          end if
549        end do
550        pawrhoij(iatom)%nrhoijsel=nselect
551        ABI_DEALLOCATE(rhoijtmp)
552      else
553        if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
554          jrhoij=1
555          do ispden=1,pawrhoij(iatom)%nspden
556            do irhoij=1,pawrhoij(iatom)%nrhoijsel
557              klmn=2*pawrhoij(iatom)%rhoijselect(irhoij)-1
558              rhoijtmp(klmn:klmn+1,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+1,ispden)
559              jrhoij=jrhoij+2
560            end do
561          end do
562        end if
563        do ispden=1,pawrhoij(iatom)%nspden
564          do kmix=1,pawrhoij(iatom)%lmnmix_sz
565            klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
566            rhoijtmp(klmn:klmn+1,ispden)=rhoijtmp(klmn:klmn+1,ispden) &
567 &           +(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
568          end do
569        end do
570        nselect=0
571        do klmn=1,pawrhoij(iatom)%lmn2_size
572          if (any(abs(rhoijtmp(2*klmn-1:2*klmn,:))>tol10)) then
573            nselect=nselect+1
574            pawrhoij(iatom)%rhoijselect(nselect)=klmn
575            do ispden=1,pawrhoij(iatom)%nspden
576              pawrhoij(iatom)%rhoijp(2*nselect-1:2*nselect,ispden)=rhoijtmp(2*klmn-1:2*klmn,ispden)
577            end do
578          end if
579        end do
580        pawrhoij(iatom)%nrhoijsel=nselect
581        ABI_DEALLOCATE(rhoijtmp)
582      end if
583 
584    end do
585  end if
586 
587 
588 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
589 !!!!! Calcul des quantites qui dependent de rho_tilde_n+1 (rho apres mixing)
590 
591  if (usepaw==1) then
592    if (ider>=0) then
593      izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero
594      call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
595 &     nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,nhatgr,&
596 &     nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,&
597 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
598 &     comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
599 &     distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
600    end if
601  end if
602 
603 !------Compute Hartree and xc potentials----------------------------------
604 
605  call hartre(1,gsqcut,usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartr)
606 
607 !Compute xc potential (separate up and down if spin-polarized)
608  optxc=1;if (nkxc>0) optxc=2
609  call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,&
610 & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,dtset%paral_kgb,rhor,rprimd,strsxc,&
611 & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
612 
613  if (nhatgrdim>0)  then
614    ABI_DEALLOCATE(nhatgr)
615  end if
616 
617 !------Compute parts of total energy depending on potentials--------
618 
619  ucvol_local = ucvol
620 
621 !Compute Hartree energy energies%e_hartree
622  call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
623 & mpi_comm_sphgrid=mpi_enreg%comm_fft)
624  energies%e_hartree=half*energies%e_hartree
625 
626 !Compute double-counting XC energy energies%e_xcdc
627  call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
628 & mpi_comm_sphgrid=mpi_enreg%comm_fft)
629 
630  etotal=energies%h0+energies%e_hartree+energies%e_xc+energies%e_corepsp + &
631 & energies%e_entropy + energies%e_elecfield + energies%e_magfield
632  etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
633  if (usepaw==1) then
634    do iatom=1,my_natom
635      itypat=paw_ij(iatom)%itypat
636      ABI_ALLOCATE(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size))
637      paw_ij(iatom)%has_dijhartree=1
638    end do
639    call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,0,dtset%ixc,my_natom,dtset%natom, &
640 &   dtset%nspden,ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang, &
641 &   dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,&
642 &   dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,&
643 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
644    do iatom=1,my_natom
645      ABI_DEALLOCATE(paw_ij(iatom)%dijhartree)
646      paw_ij(iatom)%has_dijhartree=0
647    end do
648    etotal=etotal+energies%e_paw
649  end if
650 !Compute energy residual
651  deltae=etotal-elast
652  elast=etotal
653 
654  do ispden=1,min(dtset%nspden,2)
655 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vpsp,vxc)
656    do ifft=1,nfft
657      vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
658    end do
659  end do
660  if(dtset%nspden==4) vtrial(:,3:4)=vxc(:,3:4)
661 
662  call timab(80,2,tsec)
663 
664 !DEBUG
665 !write(std_out,*) 'eeig-ehart+enxc-enxcdc+eew+eii+eent+enefield+epawdc',eeig,ehart,enxc,enxcdc,eew,eii,eent,enefield,epawdc
666 !ENDEBUG
667 
668 end subroutine odamix