Branch data 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 : : !!***
|
