TABLE OF CONTENTS


ABINIT/m_odamix [ Modules ]

[ Top ] [ Modules ]

NAME

  m_odamix

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_odamix
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_xmpi
35  use m_xcdata
36 
37  use m_time,       only : timab
38  use m_geometry,   only : metric
39  use m_cgtools,    only : dotprod_vn
40  use m_pawang,     only : pawang_type
41  use m_pawrad,     only : pawrad_type
42  use m_pawtab,     only : pawtab_type
43  use m_paw_an,     only : paw_an_type
44  use m_paw_ij,     only : paw_ij_type
45  use m_pawfgrtab,  only : pawfgrtab_type
46  use m_pawrhoij,   only : pawrhoij_type
47  use m_paw_nhat,   only : pawmknhat
48  use m_paw_denpot, only : pawdenpot
49  use m_energies,   only : energies_type
50  use m_spacepar,   only : hartre
51  use m_rhotoxc,    only : rhotoxc
52  use m_fft,        only : fourdp
53 
54  implicit none
55 
56  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)
  [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

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