Line data Source code
1 : !!****m* ABINIT/m_odamix
2 : !! NAME
3 : !! m_odamix
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (FJ, MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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 : use m_xc_tb09, only : xc_tb09_update_c
51 : use m_dft_energy, only : entropy
52 :
53 : implicit none
54 :
55 : private
56 : !!***
57 :
58 : public :: odamix
59 : !!***
60 :
61 : contains
62 : !!***
63 :
64 : !!****f* ABINIT/odamix
65 : !! NAME
66 : !! odamix
67 : !!
68 : !! FUNCTION
69 : !! This routine is called to compute the total energy and various parts of it.
70 : !! The routine computes -if requested- the forces.
71 : !!
72 : !! INPUTS
73 : !! [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
74 : !! dtset <type(dataset_type)>=all input variables in this dataset
75 : !! berryopt = 4/14: electric field is on -> add the contribution of the
76 : !! -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy
77 : !! = 6/16, or 7/17: electric displacement field is on -> add the contribution of the
78 : !! Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy
79 : !! | berryopt = 5: magnetic field is on -> add the contribution of the
80 : !! | - \Omega B.M term to the total energy
81 : !! | /= 5: magnetic field is off
82 : !! | bfield = cartesian coordinates of the magnetic field in atomic units
83 : !! | dfield = cartesian coordinates of the electric displacement field in atomic units (berryopt==6/7)
84 : !! | efield = cartesian coordinates of the electric field in atomic units (berryopt==4)
85 : !! | red_dfield = reduced the electric displacement field (berryopt==16/17)
86 : !! | red_efieldbar = reduced the electric field (ebar) (berryopt==14)
87 : !! | iatfix(3,natom)=1 for frozen atom along some direction, 0 for unfrozen
88 : !! | ionmov=governs the movement of atoms (see help file)
89 : !! | natom=number of atoms in cell.
90 : !! | nconeq=number of atomic constraint equations
91 : !! | nspden=number of spin-density components
92 : !! | nsym=number of symmetry elements in space group
93 : !! | occopt=option for occupancies
94 : !! | prtvol=integer controlling volume of printed output
95 : !! | tsmear=smearing energy or temperature (if metal)
96 : !! | wtatcon(3,natom,nconeq)=weights for atomic constraints
97 : !! | xclevel= XC functional level
98 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
99 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
100 : !! mpi_enreg=information about MPI parallelization
101 : !! my_natom=number of atoms treated by current processor
102 : !! nfft=(effective) number of FFT grid points (for this processor)
103 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
104 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
105 : !! ntypat=number of types of atoms in unit cell.
106 : !! nvresid(nfft,nspden)=potential or density residual
107 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
108 : !! optres=0 if residual array (nvresid) contains the potential residual
109 : !! =1 if residual array (nvresid) contains the density residual
110 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
111 : !! paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
112 : !! pawang <type(pawang_type)>=paw angular mesh and related data
113 : !! pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
114 : !! pawrad
115 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
116 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
117 : !! rhog(2,nfft)=array for Fourier transform of electron density
118 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3
119 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
120 : !! [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density
121 : !! ucvol = unit cell volume (Bohr**3)
122 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
123 : !! vhartr(nfft)=array for holding Hartree potential
124 : !! vpsp(nfft)=array for holding local psp
125 : !! vxc(nfft,nspden)=array for holding XC potential
126 : !! xred(3,natom)=reduced dimensionless atomic coordinates
127 : !!
128 : !! OUTPUT
129 : !! deltae=change in total energy
130 : !! between the previous and present SCF cycle
131 : !! etotal=total energy (hartree)
132 : !!
133 : !! SIDE EFFECTS
134 : !! Input/Output:
135 : !! elast=previous value of the energy,
136 : !! needed to compute deltae, then updated.
137 : !! energies <type(energies_type)>=all part of total energy.
138 : !! | entropy(IN)=entropy due to the occupation number smearing (if metal)
139 : !! | e_localpsp(IN)=local psp energy (hartree)
140 : !! | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree)
141 : !! | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
142 : !! | e_ewald(IN)=Ewald energy (hartree)
143 : !! | e_vdw_dftd(IN)=VdW DFT-D energy
144 : !! | e_hartree(IN)=Hartree part of total energy (hartree units)
145 : !! | e_corepsp(IN)=psp core-core energy
146 : !! | e_kinetic(IN)=kinetic energy part of total energy.
147 : !! | e_nucdip(IN)=energy of nuclear dipole array
148 : !! | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy.
149 : !! | e_xc(IN)=exchange-correlation energy (hartree)
150 : !! | e_xcdc(IN)=exchange-correlation double-counting energy (hartree)
151 : !! | paw%epaw(IN)=PAW spherical part energy
152 : !! | paw%epaw_dc(IN)=PAW spherical part double-counting energy
153 : !! | e_elecfield(OUT)=the term of the energy functional that depends explicitely
154 : !! | on the electric field: enefield = -ucvol*E*P
155 : !! | e_magfield(OUT)=the term of the energy functional that depends explicitely
156 : !! | on the magnetic field: enmagfield = -ucvol*B*M
157 : !! entropy=entropy due to the occupation number smearing (if metal)
158 : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0
159 : !! [vxctau(nfftf,dtset%nspden,4*usevxctau)]=derivative of XC energy density with respect to
160 : !! kinetic energy density (metaGGA cases) (optional output)
161 : !! ===== if psps%usepaw==1
162 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
163 : !! (gradients of rhoij for each atom with respect to atomic positions are computed here)
164 : !!
165 : !! NOTES
166 : !! In case of PAW calculations:
167 : !! All computations are done on the fine FFT grid.
168 : !! All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
169 : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
170 : !! ! Developpers have to be careful when introducing others arrays:
171 : !! they have to be stored on the fine FFT grid.
172 : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
173 : !!
174 : !! SOURCE
175 :
176 0 : subroutine odamix(deltae,dtset,elast,energies,etotal,&
177 0 : & gprimd,gsqcut,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,&
178 0 : & nkxc,ntypat,nvresid,n3xccc,optres,paw_ij,&
179 0 : & paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,&
180 0 : & red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,usepaw,&
181 0 : & usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,&
182 0 : & taur,vxctau,add_tfw) ! optional arguments
183 :
184 : !Arguments ------------------------------------
185 : !scalars
186 : integer,intent(in) :: my_natom,n3xccc,nfft,nkxc,ntypat,optres
187 : integer,intent(in) :: usepaw,usexcnhat
188 : logical,intent(in),optional :: add_tfw
189 : real(dp),intent(in) :: gsqcut,ucvol
190 : real(dp),intent(inout) :: elast
191 : real(dp),intent(out) :: deltae,etotal,vxcavg
192 : type(MPI_type),intent(in) :: mpi_enreg
193 : type(dataset_type),intent(in) :: dtset
194 : type(energies_type),intent(inout) :: energies
195 : type(pawang_type),intent(in) :: pawang
196 : type(pseudopotential_type),intent(in) :: psps
197 : !arrays
198 : integer,intent(in) :: ngfft(18)
199 : logical :: add_tfw_
200 : real(dp),intent(in) :: gprimd(3,3)
201 : real(dp),intent(in) :: red_ptot(3),rprimd(3,3),vpsp(nfft),xred(3,dtset%natom)
202 : real(dp),intent(in),optional :: taur(nfft,dtset%nspden*dtset%usekden)
203 : real(dp),intent(inout) :: kxc(nfft,nkxc),nhat(nfft,dtset%nspden*usepaw)
204 : real(dp),intent(inout) :: nvresid(nfft,dtset%nspden),rhog(2,nfft)
205 : real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft)
206 : real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden)
207 : real(dp),intent(inout) :: xccc3d(n3xccc)
208 : real(dp),intent(out) :: strsxc(6)
209 : real(dp),intent(inout),optional :: vxctau(:,:,:) !vxctau(nfft,dtset%nspden,4*dtset%usekden)
210 : type(paw_an_type),intent(inout) :: paw_an(my_natom)
211 : type(paw_ij_type),intent(inout) :: paw_ij(my_natom)
212 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
213 : type(pawrad_type),intent(in) :: pawrad(ntypat)
214 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
215 : type(pawtab_type),intent(in) :: pawtab(ntypat)
216 :
217 : !Local variables-------------------------------
218 : !scalars
219 : integer :: cplex,iatom,ider,idir,ierr,ifft,ipert,irhoij,ispden,itypat,izero,iir,jjr,kkr
220 : integer :: jrhoij,klmn,klmn1,kmix,nfftot,nhatgrdim,nzlmopt,nk3xc,option,optxc
221 : logical :: nmxc,with_vxctau
222 : real(dp) :: alphaopt,compch_fft,compch_sph,doti,e1t10,e_ksnm1,e_xcdc_vxctau
223 : real(dp) :: eenth,fp0,gammp1,ro_dlt,ucvol_local,el_temp
224 : character(len=500) :: message
225 : type(xcdata_type) :: xcdata
226 : !arrays
227 : real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3)
228 : real(dp) :: gmet(3,3),gprimdlc(3,3),qpt(3),rmet(3,3),tsec(2)
229 0 : real(dp),allocatable :: nhatgr(:,:,:),rhoijtmp(:,:)
230 :
231 : ! *********************************************************************
232 :
233 : !DEBUG
234 : !write(std_out,*)' odamix : enter'
235 : !ENDDEBUG
236 :
237 0 : call timab(80,1,tsec)
238 :
239 : !To be adjusted for the call to rhotoxc
240 0 : add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
241 0 : nk3xc=1;nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
242 :
243 : !faire un test sur optres=1, usewvl=0, nspden=1,nhatgrdim
244 0 : if(optres/=1)then
245 0 : write(message,'(a,i0,a)')' optres=',optres,', not allowed in oda => stop '
246 0 : ABI_ERROR(message)
247 : end if
248 :
249 0 : if(dtset%usewvl/=0)then
250 0 : write(message,'(a,i0,a)')' usewvl=',dtset%usewvl,', not allowed in oda => stop '
251 0 : ABI_ERROR(message)
252 : end if
253 :
254 0 : if(dtset%nspden/=1)then
255 0 : write(message,'(a,i0,a)')' nspden=',dtset%nspden,', not allowed in oda => stop '
256 0 : ABI_ERROR(message)
257 : end if
258 :
259 0 : if (my_natom>0) then
260 0 : if(paw_ij(1)%has_dijhat==0)then
261 0 : message = ' dijhat variable must be allocated in odamix ! '
262 0 : ABI_ERROR(message)
263 : end if
264 0 : if(paw_ij(1)%cplex_dij==2.or.paw_ij(1)%qphase==2)then
265 0 : message = ' complex dij not allowed in odamix! '
266 0 : ABI_ERROR(message)
267 : end if
268 : end if
269 :
270 : !Test size of kinetic energy potential Vxctau
271 0 : with_vxctau = (present(vxctau).and.present(taur))
272 0 : if (with_vxctau) with_vxctau = (size(vxctau)>0.and.dtset%usekden/=0)
273 : if (with_vxctau) then
274 0 : if (size(vxctau)/=nfft*dtset%nspden*4) then
275 0 : ABI_BUG("Wrong size for vxctau!")
276 : end if
277 : end if
278 :
279 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
280 : !!!!!!!!!! calculation of f'(0)= Eband_new-EH_old-E_xcdc_old-Ek_old-E_loc_old-E_nonloc_old
281 : !!!!!!!!!! save previous energy E(rho_tild_n)
282 :
283 0 : fp0=energies%e_eigenvalues-energies%h0-two*energies%e_hartree-energies%e_xcdc
284 0 : if (usepaw==1) then
285 0 : do iatom=1,my_natom
286 0 : ABI_CHECK(pawrhoij(iatom)%qphase==1,'ODA mixing not allowed with a Q phase in PAW objects!')
287 0 : itypat=pawrhoij(iatom)%itypat
288 0 : do ispden=1,pawrhoij(iatom)%nspden
289 0 : jrhoij=1
290 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
291 0 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
292 0 : ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn)
293 0 : e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden))
294 0 : jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij
295 : end do
296 0 : klmn1=1
297 0 : do klmn=1,pawrhoij(iatom)%lmn2_size
298 0 : ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,ispden)*pawtab(itypat)%dltij(klmn)
299 0 : e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,ispden)-paw_ij(iatom)%dijhat(klmn,ispden))
300 0 : klmn1=klmn1+pawrhoij(iatom)%cplex_rhoij
301 : end do
302 : end do
303 0 : if (paw_ij(iatom)%ndij>=2.and.pawrhoij(iatom)%nspden==1) then
304 0 : jrhoij=1
305 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
306 0 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
307 0 : ro_dlt=pawrhoij(iatom)%rhoijp(jrhoij,1)*pawtab(itypat)%dltij(klmn)
308 0 : e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2))
309 0 : jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij
310 : end do
311 0 : klmn1=1
312 0 : do klmn=1,pawrhoij(iatom)%lmn2_size
313 0 : ro_dlt=-pawrhoij(iatom)%rhoijres(klmn1,1)*pawtab(itypat)%dltij(klmn)
314 0 : e1t10=e1t10+ro_dlt*(paw_ij(iatom)%dij(klmn,2)-paw_ij(iatom)%dijhat(klmn,2))
315 0 : klmn1=klmn1+pawrhoij(iatom)%cplex_rhoij
316 : end do
317 0 : e1t10=half*e1t10
318 : end if
319 : end do
320 0 : if (mpi_enreg%nproc_atom>1) then
321 0 : call xmpi_sum(e1t10,mpi_enreg%comm_atom,ierr)
322 : end if
323 0 : fp0=fp0-e1t10
324 : end if
325 0 : e_ksnm1=etotal
326 :
327 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
328 : !!!!! Calculation of quantities that do not depend on rho_n+1
329 :
330 : !PAW: eventually recompute compensation density (and gradients)
331 0 : nhatgrdim=0
332 0 : if (usepaw==1) then
333 0 : ider=-1;if (dtset%xclevel==2.or.usexcnhat==0) ider=0
334 0 : if (dtset%xclevel==2.and.usexcnhat==1) ider=ider+2
335 0 : if (ider>0) then
336 0 : nhatgrdim=1
337 0 : ABI_MALLOC(nhatgr,(nfft,dtset%nspden,3))
338 : end if
339 0 : if (ider>=0) then
340 0 : ider=0;izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero
341 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
342 : nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,&
343 : & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,&
344 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
345 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
346 0 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
347 : end if
348 : end if
349 :
350 : !------Compute Hartree and xc potentials----------------------------------
351 0 : nfftot=PRODUCT(ngfft(1:3))
352 :
353 : call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,&
354 0 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
355 :
356 0 : call xcdata_init(xcdata,dtset=dtset)
357 :
358 : !If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
359 0 : if (dtset%xc_tb09_c>99._dp) then
360 : call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, &
361 : & nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, &
362 : & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,usepaw, &
363 : & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
364 0 : & computation_type='all')
365 : end if
366 :
367 : !Compute xc potential (separate up and down if spin-polarized)
368 0 : optxc=1
369 : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
370 : & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,optxc,rhor,rprimd,&
371 : & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr,&
372 0 : & vxctau=vxctau,add_tfw=add_tfw_,strsxc=strsxc)
373 :
374 : !------Compute parts of total energy depending on potentials--------
375 :
376 0 : ucvol_local=ucvol
377 :
378 : !Compute Hartree energy energies%e_hartree
379 : call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
380 0 : & mpi_comm_sphgrid=mpi_enreg%comm_fft)
381 0 : energies%e_hartree=half*energies%e_hartree
382 :
383 : !Get electronic temperature from dtset
384 0 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
385 :
386 : !Compute local psp energy energies%e_localpsp
387 : call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,&
388 0 : & mpi_comm_sphgrid=mpi_enreg%comm_fft)
389 :
390 : !Compute double-counting XC energy energies%e_xcdc
391 : call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
392 0 : & mpi_comm_sphgrid=mpi_enreg%comm_fft)
393 0 : if (with_vxctau) then
394 : call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,&
395 0 : & vxctau(:,:,1),ucvol_local,mpi_comm_sphgrid=mpi_enreg%comm_fft)
396 0 : energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
397 : end if
398 :
399 0 : if (usepaw/=0) then
400 0 : nzlmopt=dtset%pawnzlm; option=2
401 0 : do iatom=1,my_natom
402 0 : itypat=paw_ij(iatom)%itypat
403 0 : ABI_MALLOC(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size))
404 0 : paw_ij(iatom)%has_dijhartree=1
405 : end do
406 : call pawdenpot(compch_sph,el_temp,gprimd,0,dtset%ixc,my_natom,dtset%natom,dtset%nspden,ntypat,&
407 : & dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,energies%paw,paw_ij,pawang,dtset%pawprtvol,&
408 : & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,&
409 : & dtset%xc_denpos,dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis,comm_atom=mpi_enreg%comm_atom,&
410 0 : & mpi_atmtab=mpi_enreg%my_atmtab)
411 0 : do iatom=1,my_natom
412 0 : ABI_FREE(paw_ij(iatom)%dijhartree)
413 0 : paw_ij(iatom)%has_dijhartree=0
414 : end do
415 : end if
416 :
417 0 : call entropy(dtset,energies)
418 :
419 : !Turn it into an electric enthalpy,refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]]
420 : ! the missing volume is added here
421 0 : energies%e_elecfield = zero
422 0 : if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then !!HONG
423 :
424 0 : energies%e_elecfield = -dot_product(dtset%red_efieldbar,red_ptot)
425 :
426 0 : call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
427 0 : eenth = zero
428 0 : do iir=1,3
429 0 : do jjr=1,3
430 0 : eenth= eenth+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 0 : eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth
434 0 : energies%e_elecfield = energies%e_elecfield + eenth
435 :
436 : end if
437 :
438 0 : energies%e_magfield = zero
439 : !if (dtset%berryopt == 5) then
440 : !emag = dot_product(mag_cart,dtset%bfield)
441 : !energies%e_magfield = emag
442 : !end if
443 :
444 : !HONG Turn it into an internal enthalpy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]],
445 : !but a little different: U=E_ks + (vol/8*pi) * g^{-1})_ij ebar_i ebar_j
446 0 : if (dtset%berryopt == 6 .or. dtset%berryopt == 16 ) then
447 0 : energies%e_elecfield=zero
448 0 : call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
449 0 : do iir=1,3
450 0 : do jjr=1,3
451 0 : 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
452 : end do
453 : end do
454 0 : energies%e_elecfield = ucvol_local/(8.0d0*pi)*energies%e_elecfield
455 : end if
456 :
457 : !HONG calculate internal energy and electric enthalpy for mixed BC case.
458 0 : if ( dtset%berryopt == 17 ) then
459 0 : energies%e_elecfield = zero
460 0 : call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
461 0 : A(:,:)=(4*pi/ucvol_local)*rmet(:,:)
462 0 : A1(:,:)=A(:,:)
463 0 : A_new(:,:)=A(:,:)
464 0 : efield_new(:)=dtset%red_efield(:)
465 : eenth = zero
466 :
467 0 : do kkr=1,3
468 0 : if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction
469 :
470 : ! step 1 add -ebar*p
471 0 : eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr)
472 :
473 : ! step 2 chang to e_new (change e to ebar)
474 0 : efield_new(kkr)=dtset%red_efieldbar(kkr)
475 :
476 : ! step 3 chang matrix A to A1
477 :
478 0 : do iir=1,3
479 0 : do jjr=1,3
480 0 : if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
481 0 : if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) &
482 0 : & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
483 0 : if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
484 : end do
485 : end do
486 :
487 0 : A(:,:)=A1(:,:)
488 0 : A_new(:,:)=A1(:,:)
489 : end if
490 :
491 : end do ! end fo kkr
492 :
493 :
494 0 : do iir=1,3
495 0 : do jjr=1,3
496 0 : eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
497 : end do
498 : end do
499 :
500 0 : energies%e_elecfield=energies%e_elecfield+eenth
501 :
502 : end if ! berryopt==17
503 :
504 : etotal = energies%e_kinetic+ energies%e_hartree + energies%e_xc + &
505 : & energies%e_localpsp + energies%e_nlpsp_vfock - energies%e_fock0 + energies%e_corepsp + &
506 : & energies%e_entropy + energies%e_elecfield + energies%e_magfield + &
507 0 : & energies%e_nucdip
508 : !etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - &
509 : !& energies%e_xcdc + energies%e_corepsp + &
510 : !& e_entropy + energies%e_elecfield
511 0 : etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
512 0 : if (usepaw==1) then
513 0 : etotal = etotal + energies%paw%epaw
514 : end if
515 :
516 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
517 : !!!!!!!!!!!!!! now, compute mixed densities
518 :
519 0 : gammp1=etotal-e_ksnm1-fp0
520 0 : if (fp0>0.d0) then
521 0 : write(std_out,*) "fp0 est positif"
522 : ! stop
523 : end if
524 0 : write(std_out,*) "fp0 ",fp0
525 0 : alphaopt=-fp0/two/gammp1
526 :
527 0 : if (alphaopt>one.or.alphaopt<0.d0) alphaopt=one
528 0 : if (abs(energies%h0)<=tol10) alphaopt=one
529 0 : write(std_out,*) " alphaopt",alphaopt
530 :
531 0 : energies%h0=(one-alphaopt)*energies%h0 + alphaopt*(energies%e_kinetic+energies%e_localpsp)
532 0 : energies%h0=energies%h0 + alphaopt*energies%e_nlpsp_vfock
533 :
534 0 : rhor= rhor+(alphaopt-one)*nvresid
535 0 : call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,0)
536 :
537 0 : if (usepaw==1) then
538 0 : if (my_natom>0) then
539 0 : ABI_MALLOC(rhoijtmp,(pawrhoij(1)%cplex_rhoij*pawrhoij(1)%lmn2_size,pawrhoij(1)%nspden))
540 : end if
541 0 : do iatom=1,my_natom
542 0 : rhoijtmp=zero
543 0 : if (pawrhoij(iatom)%cplex_rhoij==1) then
544 0 : if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
545 0 : do ispden=1,pawrhoij(iatom)%nspden
546 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
547 0 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
548 0 : rhoijtmp(klmn,ispden)=pawrhoij(iatom)%rhoijp(irhoij,ispden)
549 : end do
550 : end do
551 : end if
552 0 : do ispden=1,pawrhoij(iatom)%nspden
553 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
554 0 : klmn=pawrhoij(iatom)%kpawmix(kmix)
555 0 : rhoijtmp(klmn,ispden)=rhoijtmp(klmn,ispden)+(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn,ispden)
556 : end do
557 : end do
558 : else
559 0 : if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
560 0 : jrhoij=1
561 0 : do ispden=1,pawrhoij(iatom)%nspden
562 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
563 0 : klmn=2*pawrhoij(iatom)%rhoijselect(irhoij)-1
564 0 : rhoijtmp(klmn:klmn+1,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+1,ispden)
565 0 : jrhoij=jrhoij+2
566 : end do
567 : end do
568 : end if
569 0 : do ispden=1,pawrhoij(iatom)%nspden
570 0 : do kmix=1,pawrhoij(iatom)%lmnmix_sz
571 0 : klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
572 : rhoijtmp(klmn:klmn+1,ispden)=rhoijtmp(klmn:klmn+1,ispden) &
573 0 : & +(alphaopt-one)*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
574 : end do
575 : end do
576 : end if
577 : call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,pawrhoij(iatom)%nrhoijsel,&
578 : & pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,pawrhoij(iatom)%lmn2_size,&
579 0 : & pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
580 : end do ! iatom
581 0 : if (allocated(rhoijtmp)) then
582 0 : ABI_FREE(rhoijtmp)
583 : end if
584 : end if ! usepaw
585 :
586 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
587 : !!!!! Calcul des quantites qui dependent de rho_tilde_n+1 (rho apres mixing)
588 :
589 0 : if (usepaw==1) then
590 0 : if (ider>=0) then
591 0 : izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero
592 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,&
593 : & nfft,ngfft,nhatgrdim,dtset%nspden,ntypat,pawang,pawfgrtab,nhatgr,&
594 : & nhat,pawrhoij,pawrhoij,pawtab,qpt,rprimd,ucvol,dtset%usewvl,xred,&
595 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
596 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
597 0 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
598 : end if
599 : end if
600 :
601 : !------Compute Hartree and xc potentials----------------------------------
602 :
603 : call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,&
604 0 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
605 :
606 : !If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
607 0 : if (dtset%xc_tb09_c>99._dp) then
608 : call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, &
609 : & nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, &
610 : & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,usepaw, &
611 : & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
612 0 : & computation_type='all')
613 : end if
614 :
615 : !Compute xc potential (separate up and down if spin-polarized)
616 0 : optxc=1;if (nkxc>0) optxc=2
617 : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
618 : & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,optxc,rhor,rprimd,&
619 : & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr,&
620 0 : & vxctau=vxctau,add_tfw=add_tfw_,strsxc=strsxc)
621 :
622 0 : if (nhatgrdim>0) then
623 0 : ABI_FREE(nhatgr)
624 : end if
625 :
626 : !------Compute parts of total energy depending on potentials--------
627 :
628 0 : ucvol_local = ucvol
629 :
630 : !Compute Hartree energy energies%e_hartree
631 : call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
632 0 : & mpi_comm_sphgrid=mpi_enreg%comm_fft)
633 0 : energies%e_hartree=half*energies%e_hartree
634 :
635 : !Compute double-counting XC energy energies%e_xcdc
636 : call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
637 0 : & mpi_comm_sphgrid=mpi_enreg%comm_fft)
638 :
639 0 : if (usepaw==1) then
640 0 : do iatom=1,my_natom
641 0 : itypat=paw_ij(iatom)%itypat
642 0 : ABI_MALLOC(paw_ij(iatom)%dijhartree,(pawtab(itypat)%lmn2_size))
643 0 : paw_ij(iatom)%has_dijhartree=1
644 : end do
645 : call pawdenpot(compch_sph,el_temp,gprimd,0,dtset%ixc,my_natom,dtset%natom,dtset%nspden,&
646 : & ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,energies%paw,paw_ij,pawang,&
647 : & dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,&
648 : & dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis,&
649 0 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
650 0 : do iatom=1,my_natom
651 0 : ABI_FREE(paw_ij(iatom)%dijhartree)
652 0 : paw_ij(iatom)%has_dijhartree=0
653 : end do
654 : end if
655 :
656 0 : call entropy(dtset,energies)
657 :
658 : etotal=energies%h0+energies%e_hartree+energies%e_xc+energies%e_corepsp + &
659 0 : & energies%e_entropy + energies%e_elecfield + energies%e_magfield
660 0 : etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
661 0 : if (usepaw==1) then
662 0 : etotal = etotal + energies%paw%epaw
663 : end if
664 :
665 : !Compute energy residual
666 0 : deltae=etotal-elast
667 0 : elast=etotal
668 :
669 0 : do ispden=1,min(dtset%nspden,2)
670 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vpsp,vxc)
671 0 : do ifft=1,nfft
672 0 : vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
673 : end do
674 : end do
675 0 : if(dtset%nspden==4) vtrial(:,3:4)=vxc(:,3:4)
676 :
677 0 : call timab(80,2,tsec)
678 :
679 0 : end subroutine odamix
680 : !!***
681 :
682 : end module m_odamix
683 : !!***
|
