Branch data Line data Source code
1 : : !!****m* ABINIT/m_setvtr
2 : : !! NAME
3 : : !! m_setvtr
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (XG, GMR, FJ, MT, EB, SPr)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_setvtr
22 : :
23 : : use defs_basis
24 : : use defs_wvltypes
25 : : use m_abicore
26 : : use m_errors
27 : : use m_abi2big
28 : : use m_xmpi
29 : : use m_xcdata
30 : : use m_dtset
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 : xred2xcart
36 : : use m_cgtools, only : dotprod_vn
37 : : use m_ewald, only : ewald
38 : : use m_energies, only : energies_type
39 : : use m_electronpositron, only : electronpositron_type, electronpositron_calctype, rhohxcpositron
40 : : use libxc_functionals, only : libxc_functionals_is_hybrid
41 : : use m_pawang, only : pawang_type
42 : : use m_pawrad, only : pawrad_type
43 : : use m_pawrhoij, only : pawrhoij_type
44 : : use m_pawtab, only : pawtab_type
45 : : use m_jellium, only : jellium
46 : : use m_spacepar, only : hartre
47 : : use m_dens, only : constrained_dft_t,constrained_dft_ini,constrained_dft_free,mag_penalty
48 : : use m_vdw_dftd2, only : vdw_dftd2
49 : : use m_vdw_dftd3, only : vdw_dftd3
50 : : use m_atm2fft, only : atm2fft
51 : : use m_rhotoxc, only : rhotoxc
52 : : use m_mklocl, only : mklocl
53 : : use m_xchybrid, only : xchybrid_ncpp_cc
54 : : use m_mkcore, only : mkcore, mkcore_alt
55 : : use m_psolver, only : psolver_rhohxc
56 : : use m_wvl_psi, only : wvl_psitohpsi
57 : : use m_mkcore_wvl, only : mkcore_wvl
58 : : use m_xc_tb09, only : xc_tb09_update_c
59 : :
60 : : #if defined HAVE_BIGDFT
61 : : use BigDFT_API, only: denspot_set_history
62 : : #endif
63 : :
64 : : implicit none
65 : :
66 : : private
67 : : !!***
68 : :
69 : : public :: setvtr
70 : : !!***
71 : :
72 : : contains
73 : : !!***
74 : :
75 : : !!****f* m_setvtr/setvtr
76 : : !! NAME
77 : : !! setvtr
78 : : !!
79 : : !! FUNCTION
80 : : !! Set up the trial potential and some energy terms
81 : : !!
82 : : !! INPUTS
83 : : !! [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
84 : : !! atindx1(dtset%natom)=index table for atoms, inverse of atindx
85 : : !! dtset <type(dataset_type)>=all input variables in this dataset
86 : : !! | if =0,1 no xc kernel, =2 spin-averaged (LDA) kernel
87 : : !! | densfor_pred=govern the choice of preconditioner for the SCF cycle
88 : : !! | iscf=determines the way the SCF cycle is handled
89 : : !! | natom=number of atoms in cell.
90 : : !! | nspden=number of spin-density components
91 : : !! | qprtrb(3)= integer wavevector of possible perturbing potential
92 : : !! | in basis of reciprocal lattice translations
93 : : !! | typat(natom)=type integer for each atom in cell
94 : : !! | vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
95 : : !! | perturbing potential is added of the form
96 : : !! | V(G)=(vprtrb(1)+I*vprtrb(2))/2 at the values G=qprtrb and
97 : : !! | (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb (integers)
98 : : !! | for each type of atom, from psp (used in norm-conserving only)
99 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
100 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
101 : : !! gsqcut=cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential)
102 : : !! istep=step number in the main loop of scfcv
103 : : !! mgfft=maximum size of 1D FFTs
104 : : !! moved_rhor=1 if the density was moved just before
105 : : !! mpi_enreg=information about MPI parallelization
106 : : !! nattyp(ntypat)=number of atoms of each type in cell.
107 : : !! nfft=(effective) number of FFT grid points (for this processor)
108 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
109 : : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
110 : : !! nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
111 : : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
112 : : !! nkxc=second dimension of the array kxc
113 : : !! ntypat=number of types of atoms in unit cell.
114 : : !! n1xccc=dimension of xccc1d; 0 if no XC core correction is used
115 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
116 : : !! optene=>0 if some additional energies have to be computed
117 : : !! pawang <type(pawang_type)> =paw angular mesh and related data
118 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
119 : : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
120 : : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
121 : : !! ph1d(2,3*(2*mgfft+1)*natom)=phase (structure factor) information.
122 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
123 : : !! rhog(2,nfft)=Fourier transform of electron density
124 : : !! rhor(nfft,nspden)=electron density in electrons/bohr**3.
125 : : !! | definition for spin components:
126 : : !! | case of nspden = 2
127 : : !! | rhor(:,1) => rho_up + rho_dwn
128 : : !! | rhor(:,2) => rho_up
129 : : !! | case of nspden = 4
130 : : !! | rhor(:,1) => rho_upup + rho_dwndwn
131 : : !! | rhor(:,2:4) => {m_x,m_y,m_z}
132 : : !! rmet(3,3)=real space metric (bohr**2)
133 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
134 : : !! ucvol = unit cell volume (bohr^3)
135 : : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
136 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
137 : : !! [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density
138 : : !!
139 : : !! OUTPUT
140 : : !! energies <type(energies_type)>=all part of total energy.
141 : : !! | e_xc=exchange-correlation energy (hartree)
142 : : !! | In case of hybrid compensation algorithm:
143 : : !! | e_hybcomp_E0=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density
144 : : !! | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density
145 : : !! | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent density
146 : : !! ==== if optene==2 or 4
147 : : !! | e_localpsp=local psp energy (hartree)
148 : : !! ==== if dtset%icoulomb == 0
149 : : !! | e_ewald=Ewald energy (hartree)
150 : : !! ==== if optene>=1
151 : : !! | e_hartree=Hartree part of total energy (hartree)
152 : : !! ==== if optene==3 or 4
153 : : !! | e_xcdc=exchange-correlation double-counting energy (hartree)
154 : : !! ==== if dtset%vdw_xc == 5 or 6 or 7
155 : : !! | e_vdw_dftd=Dispersion energy from DFT-D Van der Waals correction (hartree)
156 : : !! grchempottn(3,natom)=grads of spatially-varying chemical energy (hartree)
157 : : !! grewtn(3,natom)=grads of Ewald energy (hartree)
158 : : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree)
159 : : !! kxc(nfft,nkxc)=exchange-correlation kernel, will be computed if nkxc/=0 .
160 : : !! see routine rhotoxc for a more complete description
161 : : !! strsxc(6)=xc contribution to stress tensor (hartree/bohr^3)
162 : : !! vxcavg=mean of the vxc potential
163 : : !!
164 : : !! SIDE EFFECTS
165 : : !! moved_atm_inside=1 if the atomic positions were moved inside the SCF loop.
166 : : !! vhartr(nfft)=Hartree potential (Hartree)
167 : : !! vpsp(nfft)=local psp (Hartree)
168 : : !! vtrial(nfft,nspden)= trial potential (Hartree)
169 : : !! vxc(nfft,nspden)= xc potential (Hartree)
170 : : !! [electronpositron <type(electronpositron_type)>]=quantities for the electron-positron annihilation (optional argument)
171 : : !! [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output
172 : : !! i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density
173 : : !! [vxctau(nfftf,dtset%nspden,4*usevxctau)]=derivative of XC energy density with respect to
174 : : !! kinetic energy density (metaGGA cases) (optional output)
175 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
176 : : !! [xcctau3d(n3xccc*usekden)]=3D core electron kinetic energy density for XC core correction, bohr^-3
177 : : !!
178 : : !! NOTES
179 : : !! In case of PAW calculations:
180 : : !! All computations are done on the fine FFT grid.
181 : : !! All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
182 : : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
183 : : !! Developers have to be careful when introducing others arrays: they have to be stored on the fine FFT grid.
184 : : !!
185 : : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
186 : : !!
187 : : !! SOURCE
188 : :
189 : 18730 : subroutine setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,&
190 : 9365 : & istep,kxc,mgfft,moved_atm_inside,moved_rhor,mpi_enreg,&
191 : 28095 : & nattyp,nfft,ngfft,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,ntypat,n1xccc,n3xccc,&
192 [ + + ]: 9365 : & optene,pawang,pawrad,pawrhoij,pawtab,ph1d,psps,rhog,rhor,rmet,rprimd,strsxc,&
193 [ - + ]: 9365 : & ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,xccc3d,xred,&
194 [ + + + - : 36760 : & electronpositron,taur,vxc_hybcomp,vxctau,add_tfw,xcctau3d,calc_ewald) ! optionals arguments
- + + + +
- - + ]
195 : :
196 : : !Arguments ------------------------------------
197 : : !scalars
198 : : integer,intent(in) :: istep,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nhatgrdim,nkxc,ntypat
199 : : integer,intent(in) :: optene,usexcnhat
200 : : integer,intent(inout) :: moved_atm_inside,moved_rhor
201 : : logical,intent(in),optional :: add_tfw
202 : : logical,intent(in),optional :: calc_ewald
203 : : real(dp),intent(in) :: gsqcut,ucvol
204 : : real(dp),intent(out) :: vxcavg
205 : : type(MPI_type),intent(in) :: mpi_enreg
206 : : type(dataset_type),intent(inout) :: dtset
207 : : type(electronpositron_type),pointer,optional :: electronpositron
208 : : type(energies_type),intent(inout) :: energies
209 : : type(pawang_type),intent(in) :: pawang
210 : : type(pseudopotential_type),intent(in) :: psps
211 : : type(wvl_data), intent(inout) :: wvl
212 : : !arrays
213 : : integer, intent(in) :: atindx1(dtset%natom),nattyp(ntypat),ngfft(18)
214 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
215 : : real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw)
216 : : real(dp),intent(in) :: nhatgr(:,:,:) !(nfft,dtset%nspden,3*nhatgrdim)
217 : : real(dp),intent(in) :: rhog(2,nfft)
218 : : real(dp),intent(inout) :: rmet(3,3),rprimd(3,3)
219 : : real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
220 : : real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft)
221 : : real(dp),intent(inout),optional :: taur(nfft,dtset%nspden*dtset%usekden)
222 : : real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden)
223 : : real(dp),intent(out),optional :: vxctau(:,:,:) !vxctau(nfft,dtset%nspden,4*usevxctau)
224 : : real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden)
225 : : real(dp),intent(inout) :: xccc3d(n3xccc)
226 : : real(dp),intent(inout),optional ::xcctau3d(n3xccc*dtset%usekden)
227 : : real(dp),intent(in) :: xred(3,dtset%natom)
228 : : real(dp),intent(out) :: grchempottn(3,dtset%natom)
229 : : real(dp),intent(out) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc),strsxc(6)
230 : : type(pawrhoij_type),intent(in) :: pawrhoij(:)
231 : : type(pawrad_type),intent(in) :: pawrad(ntypat*dtset%usepaw)
232 : : type(pawtab_type),intent(in) :: pawtab(ntypat*dtset%usepaw)
233 : :
234 : : !Local variables-------------------------------
235 : : !scalars
236 : : integer :: coredens_method,coretau_method,mpi_comm_sphgrid,nk3xc
237 : : integer :: iatom,ifft,ipositron,ispden,nfftot
238 : : integer :: optatm,optdyfr,opteltfr,optgr,option,option_eff,optn,optn2,optstr,optv,vloc_method
239 : : real(dp) :: doti,e_xcdc_vxctau,ebb,ebn,evxc,ucvol_local,rpnrm,dum
240 : : logical :: add_tfw_,is_hybrid_ncpp,non_magnetic_xc,with_vxctau,wvlbigdft,lewald
241 : : character(len=500) :: message
242 : 9365 : type(constrained_dft_t) :: constrained_dft
243 : : type(xcdata_type) :: xcdata,xcdatahyb
244 : : !arrays
245 : 9365 : real(dp),allocatable :: xcart(:,:)
246 : : real(dp),parameter :: identity(1:4)=(/1._dp,1._dp,0._dp,0._dp/)
247 : : real(dp) :: dummy6(6),tsec(2)
248 : : real(dp) :: grewtn_fake(3,1)
249 : : real(dp) :: dummy_in(0)
250 : : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
251 : : real(dp) :: strn_dummy6(6), strv_dummy6(6)
252 : : real(dp) :: vhspinfield(4)
253 : 9365 : real(dp),allocatable :: grtn(:,:),dyfr_dum(:,:,:),gr_dum(:,:)
254 : 9365 : real(dp),allocatable :: rhojellg(:,:),rhojellr(:),rhowk(:,:),vjell(:)
255 : 9365 : real(dp),allocatable :: v_constr_dft_r(:,:),rhog_dum(:,:)
256 : :
257 : : ! *********************************************************************
258 : :
259 : 9365 : call timab(91,1,tsec)
260 : :
261 : : !Test size of kinetic energy potential Vxctau
262 : 9365 : with_vxctau = (present(vxctau).and.present(taur))
263 [ + + + + : 36077 : if (with_vxctau) with_vxctau = (size(vxctau)>0.and.dtset%usekden/=0)
+ + - + ]
264 : : if (with_vxctau) then
265 [ + + - + ]: 204 : if (size(vxctau)/=nfft*dtset%nspden*4) then
266 : 0 : ABI_BUG("Wrong size for vxctau!")
267 : : end if
268 : : end if
269 : :
270 : : !Check if we're in hybrid norm conserving pseudopotential with a core correction
271 : : is_hybrid_ncpp=(dtset%usepaw==0 .and. n3xccc/=0 .and. &
272 [ + + + + : 9365 : & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
+ + + + ]
273 : :
274 : : !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
275 [ - + - - ]: 18730 : wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
276 : :
277 : : !Get size of FFT grid
278 [ + + ]: 37460 : nfftot=PRODUCT(ngfft(1:3))
279 : :
280 : : !mpi
281 : 9365 : mpi_comm_sphgrid=mpi_enreg%comm_fft
282 [ - + ]: 9365 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
283 : :
284 : : !Test electron-positron case
285 [ + + ]: 9365 : ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron)
286 : :
287 : : !Test addition of Weiszacker gradient correction to Thomas-Fermi kin energy
288 [ + + ]: 9365 : add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
289 : :
290 : : !Get Ewald energy and Ewald forces, as well as vdW-DFTD energy and forces, and chemical potential energy and forces.
291 : : !-------------------------------------------------------------------------------------------------------------------
292 : 9365 : call timab(5,1,tsec)
293 [ + + ]: 9365 : if (ipositron/=1) then
294 : 9322 : lewald = .true.
295 [ + + ]: 9322 : if (present(calc_ewald)) lewald=calc_ewald
296 [ + + ]: 1755 : if (lewald) then
297 [ - + - - : 7620 : if (dtset%icoulomb == 0 .or. (dtset%usewvl == 0 .and. dtset%icoulomb == 2)) then
- - ]
298 : : ! Periodic system, need to compute energy and forces due to replica and
299 : : ! to correct the shift in potential calculation.
300 : : call ewald(energies%e_ewald,gmet,grewtn,gsqcut,dtset%icutcoul,dtset%natom,ngfft,dtset%nkpt,ntypat,&
301 : 7620 : &dtset%rcut,rmet,rprimd,dtset%typat,ucvol,dtset%vcutgeo,xred,psps%ziontypat)
302 : : ! For a periodic system bearing a finite charge, the monopole correction to the
303 : : ! energy is relevant.
304 : : ! See Leslie and Gillan, JOURNAL OF PHYSICS C-SOLID STATE PHYSICS 18, 973 (1985)
305 [ + + ]: 7620 : if(abs(dtset%cellcharge(1))>tol8) then
306 : : call ewald(energies%e_monopole,gmet,grewtn_fake,gsqcut,dtset%icutcoul,1,ngfft,dtset%nkpt,1,&
307 [ + + ]: 104 : &dtset%rcut,rmet,rprimd,(/1/),ucvol,dtset%vcutgeo,(/0.0_dp,0.0_dp,0.0_dp/),(/dtset%cellcharge(1)/))
308 : 52 : energies%e_monopole=-energies%e_monopole
309 : : end if
310 [ # # ]: 0 : else if (dtset%icoulomb == 1) then
311 : : ! In a non periodic system (real space computation), the G=0 divergence
312 : : ! doesn't occur and ewald is not needed. Only the ion/ion interaction
313 : : ! energy is relevant and used as Ewald energy and gradient.
314 : 0 : call ionion_realSpace(dtset, energies%e_ewald, grewtn, rprimd, xred, psps%ziontypat)
315 [ # # ]: 0 : else if (dtset%icoulomb == 2) then
316 : : call ionion_surface(dtset, energies%e_ewald, grewtn, mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, rprimd, &
317 : 0 : & wvl%descr, wvl%den, xred)
318 : : end if
319 : : end if
320 [ + + ]: 9322 : if (dtset%nzchempot>0) then
321 : 23 : call spatialchempot(energies%e_chempot,dtset%chempot,grchempottn,dtset%natom,ntypat,dtset%nzchempot,dtset%typat,xred)
322 : : end if
323 [ + + + - ]: 9322 : if (dtset%vdw_xc==5.and.ngrvdw==dtset%natom) then
324 : : call vdw_dftd2(energies%e_vdw_dftd,dtset%ixc,dtset%natom,ntypat,1,dtset%typat,rprimd,&
325 : 10 : & dtset%vdw_tol,xred,psps%znucltypat,gred_vdw_dftd2=grvdw)
326 : : end if
327 [ + + + - ]: 9322 : if ((dtset%vdw_xc==6.or.dtset%vdw_xc==7).and.ngrvdw==dtset%natom) then
328 : : call vdw_dftd3(energies%e_vdw_dftd,dtset%ixc,dtset%natom,&
329 : : & ntypat,1,dtset%typat,rprimd,dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,&
330 : 8 : & xred,psps%znucltypat,gred_vdw_dftd3=grvdw)
331 : : end if
332 : : else
333 : 43 : energies%e_ewald=zero
334 : 43 : energies%e_chempot=zero
335 [ + + + + ]: 1975 : grchempottn=zero
336 [ + + + + ]: 1975 : grewtn=zero
337 : 43 : energies%e_vdw_dftd=zero
338 [ - + - - : 43 : if (ngrvdw>0) grvdw=zero
- - ]
339 : : end if
340 : 9365 : call timab(5,2,tsec)
341 : :
342 : : !Compute parts of total energy depending on potentials
343 : : !--------------------------------------------------------------
344 [ + - ]: 9365 : if (dtset%usewvl == 0) then
345 : 9365 : ucvol_local = ucvol
346 : : #if defined HAVE_BIGDFT
347 : : else
348 : : ! We need to tune the volume when wavelets are used because, not all FFT points are used.
349 : : ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
350 : : ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp)
351 : : #endif
352 : : end if
353 : :
354 : : !Determine by which method the local ionic potential and/or the pseudo core charge density
355 : : ! have to be computed
356 : : !Local ionic potential:
357 : : ! Method 1: PAW
358 : : ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
359 [ + + ]: 9365 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
360 [ - + ]: 9365 : if (dtset%icoulomb>0) vloc_method=2
361 [ - + ]: 9365 : if (psps%usewvl==1) vloc_method=2
362 : : !Pseudo core charge density:
363 : : ! Method 1: PAW, nc_xccc_gspace
364 : : ! Method 2: Norm-conserving PP, wavelets
365 : : coredens_method=1;if (psps%usepaw==0) coredens_method=2
366 [ + + ]: 9365 : if (psps%nc_xccc_gspace==1) coredens_method=1
367 [ + + ]: 9365 : if (psps%nc_xccc_gspace==0) coredens_method=2
368 [ - + ]: 9365 : if (psps%usewvl==1) coredens_method=2
369 : 9365 : coretau_method=0
370 [ + + ]: 9365 : if (dtset%usekden==1) then
371 [ + + ]: 74 : coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2
372 : : end if
373 : : !In some specific cases, XC has to be handled as non-magnetic
374 [ + + + + ]: 9365 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
375 : :
376 : : !Local ionic potential and/or pseudo core charge by method 1
377 [ + + ]: 9365 : if (vloc_method==1.or.coredens_method==1) then
378 : 3162 : call timab(552,1,tsec)
379 [ + + ]: 3162 : optv=0;if (vloc_method==1) optv=1
380 [ + - ]: 3162 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
381 : 3162 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
382 : : call atm2fft(atindx1,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,&
383 : : & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
384 : : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
385 : : & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
386 : : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
387 : 3162 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
388 : 3162 : call timab(552,2,tsec)
389 : : end if
390 [ + + ]: 9365 : if (coretau_method==1) then
391 : 34 : call timab(552,1,tsec)
392 : 34 : optv=0;optn=1
393 : 34 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=4
394 : : call atm2fft(atindx1,xcctau3d,dummy_out6,dummy_out1,dummy_out2,dummy_out3,dummy_in,&
395 : : & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
396 : : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
397 : : & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
398 : : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
399 : 34 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
400 : 34 : call timab(552,2,tsec)
401 : : end if
402 : :
403 : : !Local ionic potential by method 2
404 [ + + ]: 9365 : if (vloc_method==2) then
405 : 6206 : option=1
406 [ - + + - : 18618 : ABI_MALLOC(gr_dum,(3,dtset%natom))
+ - + - -
+ - + ]
407 [ - + + - : 18618 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
+ - + - -
+ - + ]
408 [ - + + - : 18618 : ABI_MALLOC(rhog_dum,(2,nfft))
+ - + - -
+ - + ]
409 : : call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,&
410 : : & gr_dum,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,&
411 : : & nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,&
412 : 6206 : & dtset%qprtrb,rhog_dum,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
413 [ - + ]: 6206 : ABI_FREE(gr_dum)
414 [ - + ]: 6206 : ABI_FREE(dyfr_dum)
415 [ - + ]: 6206 : ABI_FREE(rhog_dum)
416 : : end if
417 : :
418 : : !3D pseudo core electron density xccc3d by method 2
419 [ + + + + ]: 9365 : if (coredens_method==2.and.n1xccc/=0) then
420 : 2867 : call timab(91,2,tsec)
421 : 2867 : call timab(92,1,tsec)
422 : 2867 : option=1
423 [ - + + - : 8601 : ABI_MALLOC(gr_dum,(3,dtset%natom))
+ - + - -
+ - + ]
424 [ - + + - : 8601 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
+ - + - -
+ - + ]
425 [ + - + - ]: 2867 : if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then
426 : : call mkcore(dummy6,dyfr_dum,gr_dum,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
427 : : & ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,&
428 : 2867 : & vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred)
429 [ # # # # : 0 : else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then
# # ]
430 : : call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,&
431 : : & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
432 : 0 : & ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred,pawrad,pawtab,psps%usepaw)
433 : : else if (psps%usewvl==1.and.psps%usepaw==1) then
434 : : #if defined HAVE_BIGDFT
435 : : ! call mkcore_wvl_old(atindx1,dummy6,dyfr_dum,wvl%descr%atoms%astruct%geocode,gr_dum,wvl%descr%h,&
436 : : ! & dtset%natom,nattyp,nfft,wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),&
437 : : ! & dtset%nspden,ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,&
438 : : ! & wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,wvl%den%denspot%dpbox%n3pi,n3xccc,option,&
439 : : ! & pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d,xred,&
440 : : ! & mpi_comm_wvl=mpi_enreg%comm_wvl)
441 : : call mkcore_wvl(atindx1,dummy6,gr_dum,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,&
442 : : & n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d,&
443 : : & psps%xcccrc,xred,wvl%den,wvl%descr,mpi_comm_wvl=mpi_enreg%comm_wvl)
444 : : #endif
445 : : end if
446 [ - + ]: 2867 : ABI_FREE(gr_dum)
447 [ - + ]: 2867 : ABI_FREE(dyfr_dum)
448 : 2867 : call timab(92,2,tsec)
449 : 2867 : call timab(91,1,tsec)
450 : : end if
451 [ + + ]: 9365 : if (coretau_method==2) then
452 : 40 : call timab(91,2,tsec)
453 : 40 : call timab(92,1,tsec)
454 : 40 : option=1
455 [ - + + - : 120 : ABI_MALLOC(gr_dum,(3,dtset%natom))
+ - + - -
+ - + ]
456 [ - + + - : 120 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
+ - + - -
+ - + ]
457 : : call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,&
458 : : & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
459 : : & ucvol,vxc,psps%xcccrc,psps%xccc1d,xcctau3d,xred,pawrad,pawtab,psps%usepaw,&
460 : 40 : & usekden=.true.)
461 [ - + ]: 40 : ABI_FREE(gr_dum)
462 [ - + ]: 40 : ABI_FREE(dyfr_dum)
463 : 40 : call timab(92,2,tsec)
464 : 40 : call timab(91,1,tsec)
465 : : end if
466 : :
467 : : !Adds the jellium potential to the local part of ionic potential
468 [ + + ]: 9365 : if (dtset%jellslab/=0) then
469 [ - + + - : 18 : ABI_MALLOC(vjell,(nfft))
+ - + - -
+ - + ]
470 [ - + + - : 18 : ABI_MALLOC(rhojellg,(2,nfft))
+ - + - -
+ - + ]
471 [ - + + - : 12 : ABI_MALLOC(rhojellr,(nfft))
- + - + ]
472 : 6 : option=1
473 : : call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,dtset%nspden,option,&
474 : 6 : & dtset%slabwsrad,rhojellg,rhojellr,rprimd,vjell,dtset%slabzbeg,dtset%slabzend)
475 : : ! Compute background-background energy
476 : 6 : call dotprod_vn(1,rhojellr,ebb,doti,nfft,nfftot,1,1,vjell,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
477 : 6 : ebb=half*ebb
478 : : ! Compute electrostatic energy between background and nuclei before adding vjell to vpsp
479 : 6 : call dotprod_vn(1,rhojellr,ebn,doti,nfft,nfftot,1,1,vpsp,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
480 : : ! Update e_ewald with ebb and ebn
481 : 6 : energies%e_ewald=energies%e_ewald+ebb+ebn
482 : : ! Compute gradient of ebn wrt tn
483 : : ! This is not yet coded for usewvl or icoulomb=1
484 [ + - - + ]: 6 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
485 : : write(message,'(3a)')&
486 : 0 : & 'The computation of forces due to jellium background',ch10,&
487 : 0 : & 'has to be verified in the PAW formalism.'
488 : 0 : ABI_WARNING(message)
489 : :
490 [ # # # # : 0 : ABI_MALLOC(grtn,(3,dtset%natom))
# # # # #
# # # ]
491 : 0 : optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=1;optn=0;optn2=1
492 : : call atm2fft(atindx1,dummy_out1,vpsp,dummy_out2,dummy_out3,dummy_out4,dummy_in,&
493 : : & gmet,gprimd,dummy_out5,grtn,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
494 : : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
495 : : & dtset%rcut,rhojellg,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
496 : : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
497 : 0 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
498 : :
499 : : ! Update grewtn with gradient of ebn wrt tn
500 [ # # ]: 0 : do iatom=1,dtset%natom
501 [ # # ]: 0 : grewtn(1:3,iatom)=grewtn(1:3,iatom)+grtn(1:3,iatom)
502 : : end do
503 [ # # ]: 0 : ABI_FREE(grtn)
504 : : else ! of usepaw==1
505 : 6 : option=2
506 [ - + + - : 18 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
+ - + - -
+ - + ]
507 [ - + + - : 18 : ABI_MALLOC(grtn,(3,dtset%natom))
+ - + - -
+ - + ]
508 : : call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,&
509 : : & grtn,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,&
510 : : & nfft,ngfft,1,ntypat,option,pawtab,ph1d,psps,dtset%qprtrb,rhojellg,&
511 : 6 : & rhojellr,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
512 : : ! Update grewtn with gradient of ebn wrt tn (reestablish order of atoms)
513 [ + + ]: 12 : do iatom=1,dtset%natom
514 [ + + ]: 30 : grewtn(1:3,atindx1(iatom))=grewtn(1:3,atindx1(iatom))+grtn(1:3,iatom)
515 : : end do
516 [ - + ]: 6 : ABI_FREE(dyfr_dum)
517 [ - + ]: 6 : ABI_FREE(grtn)
518 : : end if ! of usepaw==1
519 [ + + ]: 104886 : vpsp(:)=vpsp(:)+vjell(:)
520 [ - + ]: 6 : ABI_FREE(vjell)
521 [ - + ]: 6 : ABI_FREE(rhojellg)
522 [ - + ]: 6 : ABI_FREE(rhojellr)
523 : : end if
524 : :
525 : : !Additional stuff for electron-positron calculation
526 : : !Compute the electronic/positronic local (Hartree) potential
527 [ + + + + ]: 4917236 : if (ipositron==1) vpsp=-vpsp
528 : :
529 : : !If we are at the initialisation, or
530 : : !if the atom positions has changed and the non-linear core correction
531 : : !is included, or the rhor has changed, one needs to compute the xc stuff.
532 : : !One needs also to compute the Hartree stuff if the density changed,
533 : : !or at initialisation.
534 : : !--------------------------------------------------------------
535 : :
536 [ + + + + : 9365 : if(istep==1 .or. n1xccc/=0 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) then
+ + + + +
- ]
537 : :
538 : 9365 : option=0
539 [ + + + + : 9365 : if(istep==1 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) option=1
+ + + + ]
540 [ + + ]: 9365 : if (nkxc>0) option=2
541 [ + + ]: 9365 : if (dtset%iscf==-1) option=-2
542 : :
543 [ + + ]: 9365 : if (ipositron/=1) then
544 [ + - + - ]: 9322 : if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then
545 : :
546 : : ! >>>> Hartree potential
547 [ + + ]: 9322 : if(option/=0 .and. option/=10)then
548 : : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,&
549 : 9248 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
550 : : end if
551 : :
552 : : ! >>>> Exchange-correlation potential
553 : 9322 : call xcdata_init(xcdata,dtset=dtset)
554 [ + + ]: 9322 : if(mod(dtset%fockoptmix,100)==11)then
555 : 18 : xcdatahyb=xcdata
556 : : ! Setup the auxiliary xc functional information
557 : 18 : call xcdata_init(xcdata,dtset=dtset,auxc_ixc=0,ixc=dtset%auxc_ixc)
558 : : end if
559 : : ! Not yet able to deal fully with the full XC kernel in case of GGA + spin
560 [ + + + + : 9322 : option_eff=option;if (option==2.and.xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12
+ + ]
561 : 9322 : nk3xc=1
562 : :
563 : : ! If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
564 [ + + ]: 9322 : if (dtset%xc_tb09_c>99._dp) then
565 : : call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, &
566 : : & nfft,ngfft,nhat,psps%usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, &
567 : : & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,psps%usepaw, &
568 : : & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
569 [ - + - - : 9321 : & computation_type='all')
- - - - -
- - + ]
570 : : end if
571 : :
572 [ + + ]: 9322 : if (ipositron==0) then
573 : :
574 : : ! Compute energies%e_xc and associated quantities
575 [ + + - + ]: 9286 : if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then
576 : : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
577 : : & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
578 : : & option_eff,rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,strsxc=strsxc,&
579 [ - + - - : 10410 : & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,xcctau3d=xcctau3d,grho1_over_rho1=rpnrm)
- - - - -
- + + + +
+ + + + +
+ - + ]
580 : : else
581 : : ! Only when is_hybrid_ncpp, and moreover, the xc functional is not the auxiliary xc functional, then call xchybrid_ncpp_cc
582 : : call xchybrid_ncpp_cc(dtset,energies%e_xc,energies%entropy_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
583 : 33 : & strsxc,vxcavg,xccc3d,vxc=vxc)
584 : : end if
585 : :
586 : : ! Possibly compute energies%e_hybcomp_E0
587 [ + + ]: 9286 : if(mod(dtset%fockoptmix,100)==11)then
588 : : ! This call to rhotoxc uses the hybrid xc functional
589 [ + - ]: 18 : if(.not.is_hybrid_ncpp)then
590 : : call rhotoxc(energies%e_hybcomp_E0,dum,kxc,mpi_enreg,nfft,ngfft,&
591 : : & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
592 : : & option_eff,rhor,rprimd,usexcnhat,vxc_hybcomp,vxcavg,xccc3d,xcdatahyb,&
593 [ - + - - : 18 : & strsxc=strsxc,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
- - - - -
- - + - -
+ - - + -
+ - + - +
- - - - ]
594 : : else
595 : : call xchybrid_ncpp_cc(dtset,energies%e_hybcomp_E0,dum,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
596 [ # # # # : 0 : & strsxc,vxcavg,xccc3d,vxc=vxc_hybcomp)
# # # # #
# # # #
# ]
597 : : end if
598 : :
599 : : ! Combine hybrid and auxiliary quantities
600 : 18 : energies%e_xc=energies%e_xc*dtset%auxc_scal
601 : 18 : energies%entropy_xc=energies%entropy_xc*dtset%auxc_scal
602 : 18 : energies%e_hybcomp_E0=energies%e_hybcomp_E0-energies%e_xc
603 [ + + + + ]: 73764 : vxc(:,:)=vxc(:,:)*dtset%auxc_scal
604 [ + + + + ]: 73764 : vxc_hybcomp(:,:)=vxc_hybcomp(:,:)-vxc(:,:)
605 : : end if
606 : :
607 [ + - ]: 36 : else if (ipositron==2) then
608 : : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
609 : : & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
610 : : & option_eff,rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,strsxc=strsxc,&
611 : : & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,&
612 [ - + - - : 36 : & electronpositron=electronpositron)
- - - - -
- + - - +
- + - + ]
613 : : end if
614 : :
615 [ # # ]: 0 : elseif(.not. wvlbigdft) then
616 : : ! Use the free boundary solver
617 : : call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
618 : : & dtset%icoulomb, dtset%ixc, &
619 : : & mpi_enreg, nfft, ngfft,&
620 : : & nhat,psps%usepaw,&
621 : : & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
622 : : & usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, &
623 : : & vxcavg,wvl%descr,wvl%den,wvl%e,&
624 : 0 : & xccc3d,dtset%xclevel,dtset%xc_denpos)
625 : : end if
626 : : else
627 : 43 : energies%e_xc=zero
628 : : call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,&
629 : 43 : & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
630 : : end if
631 [ + + ]: 9365 : if (ipositron/=0) then
632 [ + - ]: 79 : if (optene>=1) then
633 : : call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,&
634 : 79 : & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
635 : : end if
636 [ + + ]: 9924334 : vhartr=vhartr+electronpositron%vha_ep
637 : : end if
638 : : end if
639 : :
640 : : !Compute the trial potential
641 : : !-------------------------------------------------------------
642 [ + - ]: 9365 : if (.not. wvlbigdft) then
643 : : ! Now, compute trial Hxc potential. Local psp potential will be added later.
644 [ + + - + ]: 9365 : if(moved_atm_inside==0 .or.dtset%iscf>=10) then
645 : :
646 : : ! Compute starting Hxc potential.
647 : : ! Multiply by identity, should not change anything if nspden /= 4
648 [ + + ]: 19672 : do ispden=1,dtset%nspden
649 [ + + ]: 235048507 : vtrial(:,ispden)=vhartr(:)*identity(ispden)+vxc(:,ispden)
650 : : end do
651 : :
652 : : else
653 : :
654 : : ! One should be here only when moved_atm_inside==1
655 : : ! The (H)xc now added corrects the previous one.
656 [ + + ]: 106 : if(dtset%densfor_pred==1)then
657 : : ! xc was substracted off. This should be rationalized later
658 [ + + ]: 86 : do ispden=1,dtset%nspden
659 [ + + ]: 461910 : vtrial(:,ispden)=vtrial(:,ispden)+vxc(:,ispden)
660 : : end do
661 [ + + + - ]: 71 : else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then
662 : : ! Hxc was substracted off. This should be rationalized later
663 [ + + ]: 84 : do ispden=1,dtset%nspden
664 [ + + ]: 1573623 : vtrial(:,ispden)=vtrial(:,ispden)+vhartr(:)*identity(ispden)+vxc(:,ispden)
665 : : end do
666 : : end if
667 : : end if
668 : :
669 : : ! Adds the local part of the potential
670 [ + + + + ]: 9365 : if ((moved_atm_inside==0).or.(dtset%densfor_pred/=3)) then
671 [ + + ]: 19526 : do ispden=1,min(2,dtset%nspden)
672 [ + + ]: 232131234 : do ifft=1,nfft
673 : 232121869 : vtrial(ifft,ispden)=vtrial(ifft,ispden)+vpsp(ifft)
674 : : end do
675 : : end do
676 : : end if
677 : :
678 : : ! Adds the compensating vxc for hybrids
679 [ + + ]: 9365 : if(mod(dtset%fockoptmix,100)==11)then
680 [ + + + + ]: 73764 : vtrial(:,:)=vtrial(:,:)+vxc_hybcomp(:,:)
681 : : end if
682 : :
683 [ - + ]: 9365 : if(dtset%usewvl==1) then
684 : 0 : call wvl_vtrial_abi2big(1,vtrial,wvl%den)
685 : : end if
686 : :
687 : : else
688 : :
689 : : ! Compute with covering comms the different part of the potential.
690 : : #if defined HAVE_BIGDFT
691 : : ! Copy e_ewald.
692 : : wvl%e%energs%eion = energies%e_ewald
693 : : ! Setup the mixing, if necessary
694 : : call denspot_set_history(wvl%den%denspot,dtset%iscf,dtset%nsppol, &
695 : : & wvl%den%denspot%dpbox%ndims(1),wvl%den%denspot%dpbox%ndims(2))
696 : : #endif
697 [ # # # # : 0 : ABI_MALLOC(xcart,(3, dtset%natom))
# # # # #
# # # ]
698 : 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
699 : : call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
700 : : & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, &
701 : : & istep, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl, dtset%nspden, &
702 : : & rpnrm, .true.,evxc, wvl,.true., xcart, strsxc,&
703 : 0 : & vtrial, vxc)
704 [ # # ]: 0 : if (optene==3.or.optene==4) energies%e_xcdc=evxc
705 [ # # ]: 0 : ABI_FREE(xcart)
706 : :
707 : : end if
708 : :
709 : : !Add the hspinfield to vtrial
710 [ + + + + : 37442 : if (any(abs(dtset%hspinfield(:))>tol8)) then
+ + ]
711 : 12 : vhspinfield(:) = zero ! vhspinfield_ij = -1/2*sigma_ij^alpha*B_alpha
712 [ + + ]: 12 : if(dtset%nspden==2)then
713 : 2 : vhspinfield(1) = -half*dtset%hspinfield(3) ! v_dwndwn = -1/2*B_z
714 : 2 : vhspinfield(2) = half*dtset%hspinfield(3) ! v_upup = 1/2*B_z
715 [ + + ]: 39368 : do ifft=1,nfft
716 : 39366 : vtrial(ifft,1) = vtrial(ifft,1) + vhspinfield(1) !SPr: added 1st component
717 : 39368 : vtrial(ifft,2) = vtrial(ifft,2) + vhspinfield(2)
718 : : end do !ifft
719 : : end if
720 [ + + ]: 12 : if(dtset%nspden==4)then
721 : 10 : vhspinfield(1)=-half*dtset%hspinfield(3) ! v_dwndwn => v_11
722 : 10 : vhspinfield(2)= half*dtset%hspinfield(3) ! v_upup => v_22
723 : 10 : vhspinfield(3)=-half*dtset%hspinfield(1) ! Re(v_dwnup) = Re(v_updwn) => Re(v_12)
724 : 10 : vhspinfield(4)= half*dtset%hspinfield(2) ! Im(v_dwnup) =-Im(v_dwnup) => Im(v_12)
725 [ + + ]: 50 : do ispden=1,dtset%nspden
726 [ + + ]: 787370 : do ifft=1,nfft
727 : 787360 : vtrial(ifft,ispden) = vtrial(ifft,ispden) + vhspinfield(ispden)
728 : : end do
729 : : end do
730 : : end if
731 : : end if
732 : :
733 : : !Compute the constrained potential for the magnetic moments
734 [ + + ]: 9365 : if (dtset%magconon==1.or.dtset%magconon==2) then
735 : : ! Initialize the datastructure constrained_dft, for penalty function constrained magnetization
736 : : call constrained_dft_ini(dtset%chrgat,constrained_dft,dtset%constraint_kind,dtset%magconon,dtset%magcon_lambda,&
737 : : & mpi_enreg,dtset%natom,dtset%nfft,dtset%ngfft,dtset%nspden,dtset%ntypat,dtset%ratsm,&
738 : 5 : & dtset%ratsph,rprimd,dtset%spinat,dtset%typat,xred,dtset%ziontypat,dtset%znucl,dtset%qgbt,dtset%use_gbt,dtset%spinaxis)
739 [ - + + - : 20 : ABI_MALLOC(v_constr_dft_r, (nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
740 [ + + + + ]: 65557 : v_constr_dft_r = zero
741 : 5 : call mag_penalty(constrained_dft,mpi_enreg,rhor,v_constr_dft_r,xred,dtset%qgbt,dtset%use_gbt)
742 [ + + ]: 5 : if(dtset%nspden==4)then
743 [ + + ]: 15 : do ispden=1,dtset%nspden ! (SPr: both components should be used? EB: Yes it should be the case, corrected now)
744 [ + + ]: 49167 : do ifft=1,nfft
745 : 49164 : vtrial(ifft,ispden) = vtrial(ifft,ispden) + v_constr_dft_r(ifft,ispden)
746 : : end do !ifft
747 : : end do !ispden
748 [ + - ]: 2 : else if(dtset%nspden==2)then
749 [ + + ]: 8194 : do ifft=1,nfft
750 : : ! TODO : MJV: check that magnetic constraint works also for nspden 2 or add input variable condition
751 : : ! EB: ispden=2 is rho_up only: to be tested
752 : : ! SPr: for ispden=2, both components should be used (e.g. see definition for vhspinfield)?
753 : 8192 : vtrial(ifft,1) = vtrial(ifft,1) + v_constr_dft_r(ifft,1) !SPr: added the first component here
754 : 8194 : vtrial(ifft,2) = vtrial(ifft,2) + v_constr_dft_r(ifft,2)
755 : : end do !ifft
756 : : end if
757 [ - + ]: 5 : ABI_FREE(v_constr_dft_r)
758 : 5 : call constrained_dft_free(constrained_dft)
759 : : end if
760 : :
761 : : !Compute parts of total energy depending on potentials
762 : : !--------------------------------------------------------------
763 : :
764 : : !For icoulomb==0 or usewvl Ehartree is calculated in psolver_rhohxc().
765 : : !For PAW we recalculate this since nhat was not taken into account
766 : : !in psolver_rhohxc: E_H= int v_H (n+nhat) dr
767 : :
768 [ + + + - : 9365 : if (optene>=1 .and. .not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then
- + - - ]
769 : : ! Compute Hartree energy ehart
770 : : ! Already available in the Psolver case through psolver_rhohxc().
771 [ + + ]: 3750 : if (ipositron/=1) then
772 : : call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
773 : 3707 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
774 [ + + ]: 3707 : if (ipositron==0) energies%e_hartree = half * energies%e_hartree
775 [ + + ]: 3707 : if (ipositron==2) energies%e_hartree = half * (energies%e_hartree-electronpositron%e_hartree)
776 : : else
777 : 43 : energies%e_hartree=zero
778 : : end if
779 : : end if
780 : :
781 [ + + ]: 9365 : if(mod(dtset%fockoptmix,100)==11)then
782 [ + - ]: 18 : if (.not. wvlbigdft) then
783 : : call dotprod_vn(1,rhor,energies%e_hybcomp_v0,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol_local,&
784 [ - + - - : 18 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
- - - - -
+ ]
785 : 18 : energies%e_hybcomp_v=energies%e_hybcomp_v0
786 : : end if
787 : : end if
788 : :
789 [ + - + + : 9365 : if (optene==2.or.optene==4 .and. .not. wvlbigdft) then
+ - ]
790 : : ! Compute local psp energy eei
791 : : call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,&
792 : 3750 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
793 : : end if
794 : :
795 [ + - + + : 9365 : if (optene==3.or.optene==4 .and. .not. wvlbigdft) then
+ - ]
796 : : ! Compute double-counting XC energy enxcdc
797 [ + + ]: 3750 : if (ipositron/=1) then
798 [ + + + + ]: 3707 : if (dtset%usepaw==0.or.usexcnhat/=0) then
799 : : call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
800 : 1323 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
801 : : else
802 [ - + + - : 9536 : ABI_MALLOC(rhowk,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
803 [ + - + - : 56135544 : rhowk=rhor-nhat
- + - - +
+ + + ]
804 : : call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
805 : 2384 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
806 [ - + ]: 2384 : ABI_FREE(rhowk)
807 : : end if
808 [ + + ]: 3707 : if (with_vxctau) then
809 : : call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),&
810 [ - + - - : 34 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
- - - - -
+ ]
811 : 34 : energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
812 : : end if
813 [ + + ]: 3707 : if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
814 : : else
815 : 43 : energies%e_xcdc=zero
816 : : end if
817 : : end if
818 : :
819 : : !--------------------------------------------------------------
820 : :
821 : : !The initialisation for the new atomic positions has been done
822 : 9365 : moved_atm_inside=0
823 : :
824 : 9365 : call timab(91,2,tsec)
825 : :
826 [ + + + + : 28066 : end subroutine setvtr
- + - + -
+ - + - +
- + - + -
+ ]
827 : : !!***
828 : :
829 : : !!****m* m_setvtr/spatialchempot
830 : : !! NAME
831 : : !! spatialchempot
832 : : !! FUNCTION
833 : : !! Treat spatially varying chemical potential.
834 : : !! Compute energy and derivative with respect to dimensionless reduced atom coordinates of the
835 : : !! spatially varying chemical potential. No contribution to stresses.
836 : : !!
837 : : !! INPUTS
838 : : !! chempot(3,nzchempot,ntypat)=input array with information about the chemical potential (see input variable description)
839 : : !! natom=number of atoms in unit cell
840 : : !! ntypat=number of type of atoms
841 : : !! nzchempot=number of limiting planes for chemical potential
842 : : !! typat(natom)=integer label of each type of atom (1,2,...)
843 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
844 : : !!
845 : : !! OUTPUT
846 : : !! e_chempot=chemical potential energy in hartrees
847 : : !! grchempottn(3,natom)=grads of e_chempot wrt xred(3,natom), hartrees.
848 : : !!
849 : : !! SOURCE
850 : :
851 : 23 : subroutine spatialchempot(e_chempot,chempot,grchempottn,natom,ntypat,nzchempot,typat,xred)
852 : :
853 : : !Arguments ------------------------------------
854 : : !scalars
855 : : integer,intent(in) :: natom,ntypat,nzchempot
856 : : real(dp),intent(out) :: e_chempot
857 : : !arrays
858 : : integer,intent(in) :: typat(natom)
859 : : real(dp),intent(in) :: chempot(3,nzchempot,ntypat),xred(3,natom)
860 : : real(dp),intent(out) :: grchempottn(3,natom)
861 : :
862 : : !Local variables-------------------------------
863 : : !scalars
864 : : integer :: iatom,itypat,iz
865 : : real(dp) :: a_2,a_3,cp0,cp1,dcp0,dcp1,ddz,deltaz,deltaziz
866 : : real(dp) :: dqz,dz1,qz,zred,z0
867 : : !character(len=500) :: message
868 : :
869 : : ! *************************************************************************
870 : :
871 : : !DEBUG
872 : : !write(std_out,'(a)')' spatialchempot : enter '
873 : : !write(std_out,'(a,2i6)')' nzchempot,ntypat',nzchempot,ntypat
874 : : !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,1)=',chempot(1:3,1:2,1)
875 : : !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,2)=',chempot(1:3,1:2,2)
876 : : !ENDDEBUG
877 : :
878 : 23 : e_chempot=zero
879 [ + + + + ]: 207 : grchempottn(:,:)=zero
880 : :
881 : : !Loop on the different atoms
882 [ + + ]: 69 : do iatom=1,natom
883 : :
884 : 46 : itypat=typat(iatom)
885 : 46 : zred=xred(3,iatom)
886 : :
887 : : ! Determine the delimiting plane just lower to zred
888 : : ! First compute the relative zred with respect to the first delimiting plane
889 : : ! Take into account a tolerance :
890 : 46 : deltaz=zred-chempot(1,1,itypat)
891 [ + - - + ]: 46 : deltaz=modulo(deltaz+tol12,1.0d0)-tol12
892 : : ! deltaz is positive (or higher than -tol12), and lower than one-tol12.
893 [ + - ]: 56 : do iz=2,nzchempot+1
894 [ + + ]: 56 : if(iz/=nzchempot+1)then
895 : 46 : deltaziz=chempot(1,iz,itypat)-chempot(1,1,itypat)
896 : : else
897 : : deltaziz=one
898 : : end if
899 [ + + ]: 56 : if(deltaziz>deltaz)exit
900 : : end do
901 : :
902 : : ! Defines coordinates and values inside the delimiting interval,
903 : : ! with respect to the lower delimiting plane
904 : 46 : z0=chempot(1,iz-1,itypat)-chempot(1,1,itypat) ; cp0=chempot(2,iz-1,itypat) ; dcp0=chempot(3,iz-1,itypat)
905 [ + + ]: 46 : if(iz/=nzchempot+1)then
906 : 36 : dz1=chempot(1,iz,itypat)-chempot(1,iz-1,itypat) ; cp1=chempot(2,iz,itypat) ; dcp1=chempot(3,iz,itypat)
907 : : else
908 : 10 : dz1=(chempot(1,1,itypat)+one)-chempot(1,nzchempot,itypat) ; cp1=chempot(2,1,itypat) ; dcp1=chempot(3,1,itypat)
909 : : end if
910 : 46 : ddz=deltaz-z0
911 : :
912 : : !DEBUG
913 : : ! write(std_out,'(a,2i5)')' Delimiting planes, iz-1 and iz=',iz-1,iz
914 : : ! write(std_out,'(a,2es13.3)')' z0, dz1= :',z0,dz1
915 : : ! write(std_out,'(a,2es13.3)')' cp0, cp1= :',cp0,cp1
916 : : ! write(std_out,'(a,2es13.3)')' dcp0, dcp1= :',dcp0,dcp1
917 : : ! write(std_out,'(a,2es13.3)')' deltaz,ddz=',deltaz,ddz
918 : : !ENDDEBUG
919 : :
920 : : ! Determine the coefficient of the third-order polynomial taking z0 as origin
921 : : ! P(dz=z-z0)= a_3*dz**3 + a_2*dz**2 + a_1*dz + a_0 ; obviously a_0=cp0 and a_1=dcp0
922 : : ! Define qz=a_3*dz + a_2 and dqz=3*a_3*dz + 2*a_2
923 : 46 : qz=((cp1-cp0)-dcp0*dz1)/dz1**2
924 : 46 : dqz=(dcp1-dcp0)/dz1
925 : 46 : a_3=(dqz-two*qz)/dz1
926 : 46 : a_2=three*qz-dqz
927 : :
928 : : ! Compute value and gradient of the chemical potential, at ddz wrt to lower delimiting plane
929 : 46 : e_chempot=e_chempot+(a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0)
930 : 69 : grchempottn(3,iatom)=three*a_3*ddz**2 + two*a_2*ddz + dcp0
931 : :
932 : : !DEBUG
933 : : ! write(std_out,'(a,4es16.6)')' qz,dqz=',qz,dqz
934 : : ! write(std_out,'(a,4es16.6)')' cp0,dcp0,a_2,a_3=',cp0,dcp0,a_2,a_3
935 : : ! write(std_out,'(a,2es13.3)')' dcp0*ddz + cp0=',dcp0*ddz + cp0
936 : : ! write(std_out,'(a,2es13.3)')' a_2*ddz**2=',a_2*ddz**2
937 : : ! write(std_out,'(a,2es13.3)')' a_3*ddz**3=',a_3*ddz**3
938 : : ! write(std_out,'(a,2es13.3)')' contrib=',a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0
939 : : ! write(std_out,'(a,2es13.3)')' e_chempot=',e_chempot
940 : : ! write(std_out,'(a,3es20.10)')' grchempottn=',grchempottn(:,iatom)
941 : : !ENDDEBUG
942 : :
943 : : end do
944 : :
945 : : !DEBUG
946 : : !write(std_out,'(a)')' spatialchempot : exit '
947 : : !write(std_out,'(a,es16.6)') ' e_chempot=',e_chempot
948 : : !ENDDEBUG
949 : :
950 : 23 : end subroutine spatialchempot
951 : : !!***
952 : :
953 : : !!****f* ABINIT/ionion_realspace
954 : : !!
955 : : !! NAME
956 : : !! ionion_realspace
957 : : !!
958 : : !! FUNCTION
959 : : !! Compute the ion/ion interaction energies and forces in real space
960 : : !! case. Use ewald() instead if computations are done in reciprocal
961 : : !! space since it also includes the correction for the shift done in
962 : : !! potentials calculations and includes replica interactions.
963 : : !!
964 : : !! INPUTS
965 : : !! dtset <type(dataset_type)>=all input variables in this dataset
966 : : !! rmet(3,3)=metric tensor in real space (bohr^2)
967 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
968 : : !! zion(ntypat)=charge on each type of atom (real number)
969 : : !!
970 : : !! OUTPUT
971 : : !! eew=final ion/ion energy in hartrees
972 : : !! grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.
973 : : !!
974 : : !! SOURCE
975 : :
976 : 0 : subroutine ionion_realSpace(dtset, eew, grewtn, rprimd, xred, zion)
977 : :
978 : : !Arguments ------------------------------------
979 : : !scalars
980 : : real(dp),intent(out) :: eew
981 : : type(dataset_type),intent(in) :: dtset
982 : : !arrays
983 : : real(dp),intent(in) :: rprimd(3,3),zion(dtset%ntypat)
984 : : real(dp),intent(in) :: xred(3,dtset%natom)
985 : : real(dp),intent(out) :: grewtn(3,dtset%natom)
986 : :
987 : : !Local variables-------------------------------
988 : : !scalars
989 : : integer :: ia1,ia2,iatom,igeo
990 : : real(dp) :: r
991 : : !arrays
992 : 0 : real(dp),allocatable :: grew_cart(:,:),xcart(:,:)
993 : :
994 : : ! *************************************************************************
995 : :
996 : : !Store xcart for each atom
997 [ # # # # : 0 : ABI_MALLOC(xcart,(3, dtset%natom))
# # # # #
# # # ]
998 : 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
999 : :
1000 : : !Summing the interaction between ions.
1001 : 0 : eew = 0._dp
1002 [ # # ]: 0 : do ia1 = 1, dtset%natom, 1
1003 [ # # ]: 0 : do ia2 = ia1 + 1, dtset%natom, 1
1004 : : r = sqrt((xcart(1, ia1) - xcart(1, ia2)) ** 2 + &
1005 : : & (xcart(2, ia1) - xcart(2, ia2)) ** 2 + &
1006 : 0 : & (xcart(3, ia1) - xcart(3, ia2)) ** 2)
1007 : 0 : eew = eew + zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / r
1008 : : end do
1009 : : end do
1010 : :
1011 : : !Allocate temporary array to store cartesian gradients.
1012 [ # # # # : 0 : ABI_MALLOC(grew_cart,(3, dtset%natom))
# # # # ]
1013 : :
1014 : : !Summing the forces for each atom
1015 [ # # ]: 0 : do ia1 = 1, dtset%natom, 1
1016 [ # # ]: 0 : grew_cart(:, ia1) = 0._dp
1017 [ # # ]: 0 : do ia2 = 1, dtset%natom, 1
1018 [ # # ]: 0 : if (ia1 /= ia2) then
1019 : : r = (xcart(1, ia1) - xcart(1, ia2)) ** 2 + &
1020 : : & (xcart(2, ia1) - xcart(2, ia2)) ** 2 + &
1021 : 0 : & (xcart(3, ia1) - xcart(3, ia2)) ** 2
1022 [ # # ]: 0 : do igeo = 1, 3, 1
1023 : : grew_cart(igeo, ia1) = grew_cart(igeo, ia1) - (xcart(igeo, ia1) - xcart(igeo, ia2)) * &
1024 : 0 : & zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / (r ** 1.5_dp)
1025 : : end do
1026 : : end if
1027 : : end do
1028 : : end do
1029 : :
1030 [ # # ]: 0 : ABI_FREE(xcart)
1031 : :
1032 : : !Transform cartesian gradients to reduced gradients.
1033 [ # # ]: 0 : do iatom = 1, dtset%natom, 1
1034 [ # # ]: 0 : do igeo = 1, 3, 1
1035 : : grewtn(igeo, iatom) = rprimd(1, igeo) * grew_cart(1, iatom) + &
1036 : : & rprimd(2, igeo) * grew_cart(2, iatom) + &
1037 : 0 : & rprimd(3, igeo) * grew_cart(3, iatom)
1038 : : end do
1039 : : end do
1040 [ # # ]: 0 : ABI_FREE(grew_cart)
1041 : :
1042 : 0 : end subroutine ionion_realSpace
1043 : : !!***
1044 : :
1045 : : !!****f* ABINIT/ionion_surface
1046 : : !!
1047 : : !! NAME
1048 : : !! ionion_surface
1049 : : !!
1050 : : !! FUNCTION
1051 : : !! Compute the ion/ion interaction energies and forces in real space
1052 : : !! case. Use ewald() instead if computations are done in reciprocal
1053 : : !! space since it also includes the correction for the shift done in
1054 : : !! potentials calculations and includes replica interactions.
1055 : : !!
1056 : : !! INPUTS
1057 : : !! dtset <type(dataset_type)>=all input variables in this dataset
1058 : : !! rmet(3,3)=metric tensor in real space (bohr^2)
1059 : : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
1060 : : !! zion(ntypat)=charge on each type of atom (real number)
1061 : : !!
1062 : : !! OUTPUT
1063 : : !! eew=final ion/ion energy in hartrees
1064 : : !! grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.
1065 : : !!
1066 : : !! SOURCE
1067 : :
1068 : 0 : subroutine ionion_surface(dtset, eew, grewtn, me, nproc, rprimd, wvl, wvl_den, xred)
1069 : :
1070 : : #if defined HAVE_BIGDFT
1071 : : use BigDFT_API, only: IonicEnergyandForces
1072 : : #endif
1073 : :
1074 : : !Arguments ------------------------------------
1075 : : !scalars
1076 : : integer, intent(in) :: me, nproc
1077 : : real(dp),intent(out) :: eew
1078 : : type(dataset_type),intent(in) :: dtset
1079 : : type(wvl_internal_type), intent(in) :: wvl
1080 : : type(wvl_denspot_type), intent(inout) :: wvl_den
1081 : : !arrays
1082 : : real(dp),intent(in) :: rprimd(3,3)
1083 : : real(dp),intent(in) :: xred(3,dtset%natom)
1084 : : real(dp),intent(out) :: grewtn(3,dtset%natom)
1085 : :
1086 : : !Local variables-------------------------------
1087 : : !scalars
1088 : : integer :: dispersion, iatom, igeo
1089 : : real(dp) :: psoffset
1090 : : !arrays
1091 : 0 : real(dp),allocatable :: xcart(:,:)
1092 : 0 : real(dp),pointer :: grew_cart(:,:),fdisp(:,:)
1093 : : #if defined HAVE_BIGDFT
1094 : : real(dp) :: edisp
1095 : : real(dp) :: ewaldstr(6)
1096 : : #endif
1097 : :
1098 : : ! *************************************************************************
1099 : :
1100 : : !Store xcart for each atom
1101 [ # # # # : 0 : ABI_MALLOC(xcart,(3, dtset%natom))
# # # # #
# ]
1102 : 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
1103 : :
1104 : 0 : nullify(fdisp)
1105 : 0 : nullify(grew_cart)
1106 : 0 : dispersion = 0
1107 : 0 : psoffset = 0._dp
1108 : : #if defined HAVE_BIGDFT
1109 : : call IonicEnergyandForces(me, nproc, wvl_den%denspot%dpbox,&
1110 : : & wvl%atoms, dtset%efield, xcart, &
1111 : : & eew, grew_cart, dispersion, edisp, fdisp,&
1112 : : & ewaldstr,wvl%Glr%d%n1,wvl%Glr%d%n2,wvl%Glr%d%n3,&
1113 : : & wvl_den%denspot%V_ext, wvl_den%denspot%pkernel,psoffset)
1114 : :
1115 : : if (associated(fdisp)) then
1116 : : ABI_FREE(fdisp)
1117 : : end if
1118 : : #endif
1119 : :
1120 [ # # ]: 0 : ABI_FREE(xcart)
1121 : :
1122 : : !Transform cartesian gradients to reduced gradients.
1123 [ # # ]: 0 : do iatom = 1, dtset%natom, 1
1124 [ # # ]: 0 : do igeo = 1, 3, 1
1125 : : grewtn(igeo, iatom) = -rprimd(1, igeo) * grew_cart(1, iatom) - &
1126 : : & rprimd(2, igeo) * grew_cart(2, iatom) - &
1127 : 0 : & rprimd(3, igeo) * grew_cart(3, iatom)
1128 : : end do
1129 : : end do
1130 : : if (associated(grew_cart)) then
1131 : : ABI_FREE(grew_cart)
1132 : : end if
1133 : :
1134 : : #if !defined HAVE_BIGDFT
1135 : : if (.false.) write(std_out,*) me,nproc,wvl%h(1),wvl_den%symObj
1136 : : #endif
1137 : :
1138 [ # # ]: 0 : end subroutine ionion_surface
1139 : : !!***
1140 : :
1141 : : end module m_setvtr
1142 : : !!***
|
