TABLE OF CONTENTS


ABINIT/m_odamix [ Modules ]

[ Top ] [ Modules ]

NAME

  m_odamix

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_odamix
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_xmpi
28  use m_xcdata
29  use m_dtset
30 
31 
32  use defs_datatypes, only : pseudopotential_type
33  use defs_abitypes,  only : MPI_type
34  use m_time,       only : timab
35  use m_geometry,   only : metric
36  use m_cgtools,    only : dotprod_vn
37  use m_pawang,     only : pawang_type
38  use m_pawrad,     only : pawrad_type
39  use m_pawtab,     only : pawtab_type
40  use m_paw_an,     only : paw_an_type
41  use m_paw_ij,     only : paw_ij_type
42  use m_pawfgrtab,  only : pawfgrtab_type
43  use m_pawrhoij,   only : pawrhoij_type,pawrhoij_filter
44  use m_paw_nhat,   only : pawmknhat
45  use m_paw_denpot, only : pawdenpot
46  use m_energies,   only : energies_type
47  use m_spacepar,   only : hartre
48  use m_rhotoxc,    only : rhotoxc
49  use m_fft,        only : fourdp
50 
51  implicit none
52 
53  private

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.

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)
  [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_nucdip(IN)=energy of nuclear dipole array
   | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE 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.

SOURCE

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