Line data Source code
1 : !!****m* ABINIT/m_forces
2 : !! NAME
3 : !! m_forces
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, FJ, MM, MT, SCE)
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 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_forces
26 :
27 : use defs_basis
28 : use defs_wvltypes
29 : use m_abicore
30 : use m_efield
31 : use m_errors
32 : use m_atomdata
33 : use m_dtset
34 :
35 : use defs_datatypes, only : pseudopotential_type
36 : use defs_abitypes, only : MPI_type
37 : use m_time, only : timab
38 : use m_geometry, only : gred2fcart, metric, xred2xcart
39 : use m_fock, only : fock_type
40 : use m_pawrad, only : pawrad_type
41 : use m_pawtab, only : pawtab_type
42 : use m_electronpositron, only : electronpositron_type,electronpositron_calctype
43 : use libxc_functionals, only : libxc_functionals_is_hybrid
44 : use m_fft, only : zerosym, fourdp
45 : use m_cgtools, only : mean_fftr
46 : use m_mpinfo, only : pre_gather, pre_scatter
47 : use m_atm2fft, only : atm2fft
48 : use m_mklocl, only : mklocl
49 : use m_predtk, only : prtxvf
50 : use m_xchybrid, only : xchybrid_ncpp_cc
51 : use m_mkcore, only : mkcore, mkcore_alt
52 : use m_mkcore_wvl, only : mkcore_wvl
53 :
54 : #if defined(HAVE_GPU_MARKERS)
55 : use m_nvtx_data
56 : #endif
57 :
58 : implicit none
59 :
60 : private
61 : !!***
62 :
63 : public :: forces
64 : public :: fresid
65 : !!***
66 :
67 : contains
68 : !!***
69 :
70 : !!****f* ABINIT/forces
71 : !! NAME
72 : !! forces
73 : !!
74 : !! FUNCTION
75 : !! Assemble gradients of various total energy terms with respect
76 : !! to reduced coordinates, including possible symmetrization,
77 : !! in order to produce forces.
78 : !!
79 : !! fcart(i,iat) = d(Etot)/(d(r(i,iat)))
80 : !!
81 : !! INPUTS
82 : !! atindx1(natom)=index table for atoms, inverse of atindx
83 : !! dtefield <type(efield_type)> = variables related to Berry phase
84 : !! dtset <type(dataset_type)>=all input variables in this dataset
85 : !! berryopt = 4/14: electric field is on -> add the contribution of the
86 : !! -ebar_i p_i - Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy
87 : !! = 6/16, or 7/17: electric displacement field is on -> add the contribution of the
88 : !! Omega/(8*pi) (g^{-1})_ij ebar_i ebar_j terms to the total energy
89 : !! | efield = cartesian coordinates of the electric field in atomic units
90 : !! | dfield = cartesian coordinates of the electric displacement field in atomic units
91 : !! | iatfix(3,natom)=1 for frozen atom along specified direction, 0 for unfrozen
92 : !! | ionmov=governs the movement of atoms (see help file)
93 : !! | densfor_pred=governs the mixed electronic-atomic part of the preconditioner
94 : !! | natom=number of atoms in cell
95 : !! | nconeq=number of atomic constraint equations
96 : !! | nspden=number of spin-density components
97 : !! | nsym=number of symmetries in space group
98 : !! | prtvol=integer controlling volume of printed output
99 : !! | typat(natom)=type integer for each atom in cell
100 : !! | wtatcon(3,natom,nconeq)=weights for atomic constraints
101 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
102 : !! grchempottn(3,natom)=d(E_chemical potential)/d(xred) (hartree)
103 : !! grcondft(3,natom)=d(E_constrainedDFT)/d(xred) (hartree)
104 : !! grewtn(3,natom)=d(Ewald)/d(xred) (hartree)
105 : !! grnl(3*natom)=gradients of Etot due to nonlocal contributions
106 : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree)
107 : !! gsqcut=cutoff value on G**2 for (large) sphere inside FFT box.
108 : !! gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)
109 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
110 : !! mgfft=maximum size of 1D FFTs
111 : !! mpi_enreg=information about MPI parallelization
112 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
113 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
114 : !! nattyp(ntypat)=number of atoms of each type
115 : !! nfft=(effective) number of FFT grid points (for this processor)
116 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
117 : !! ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
118 : !! ntypat=number of types of atoms
119 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
120 : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
121 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) array
122 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
123 : !! rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
124 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3
125 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
126 : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
127 : !! usefock=1 if fock operator is used; 0 otherwise.
128 : !! usevxctau=1 if if XC functional depends on kinetic energy density
129 : !! vresid(nfft,nspden)=potential residual (if non-collinear magn., only trace of it)
130 : !! vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space
131 : !! vxctau(nfft,nspden,4*usevxctau)=(only for meta-GGA): derivative of XC energy density
132 : !! wrt kinetic energy density (depsxcdtau)
133 : !! xred(3,natom)=reduced dimensionless atomic coordinates
134 : !! xred_old(3,natom)=previous reduced dimensionless atomic coordinates
135 : !!
136 : !! OUTPUT
137 : !! diffor=maximal absolute value of changes in the components of
138 : !! force between the input and the output.
139 : !! favg(3)=mean of the forces before correction for translational symmetry
140 : !! forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
141 : !! gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
142 : !! gresid(3,natom)=forces due to the residual of the density/potential
143 : !! grhf(3,natom)=Hellman-Feynman derivatives of the total energy
144 : !! grxc(9+3*natom)=d(Exc)/d(xred) if core charges are used
145 : !! maxfor=maximal absolute value of the output array force.
146 : !! synlgr(3,natom)=symmetrized d(enl)/d(xred)
147 : !!
148 : !! SIDE EFFECTS
149 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
150 : !! fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
151 : !! Note : unlike gred, this array has been corrected by enforcing
152 : !! the translational symmetry, namely that the sum of force
153 : !! on all atoms is zero.
154 : !!
155 : !! NOTES
156 : !! * Symmetrization of gradients with respect to reduced
157 : !! coordinates xred is conducted according to the expression
158 : !! [d(e)/d(t(n,a))]_symmetrized = (1/Nsym) Sum(S) symrec(n,m,S)*
159 : !! [d(e)/d(t(m,b))]_unsymmetrized
160 : !! where t(m,b)= (symrel^-1)(m,n)*(t(n,a)-tnons(n)) and tnons
161 : !! is a possible nonsymmorphic translation. The label "b" here
162 : !! refers to the atom which gets rotated into "a" under symmetry "S".
163 : !! symrel is the symmetry matrix in real space, which is the inverse
164 : !! transpose of symrec. symrec is the symmetry matrix in reciprocal
165 : !! space. sym_cartesian = R * symrel * R^-1 = G * symrec * G^-1
166 : !! where the columns of R and G are the dimensional primitive translations
167 : !! in real and reciprocal space respectively.
168 : !! * Note the use of "symrec" in the symmetrization expression above.
169 : !!
170 : !! SOURCE
171 :
172 23761 : subroutine forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,&
173 23761 : & forold,gred,grchempottn,grcondft,gresid,grewtn,&
174 23761 : & grhf,grnl,grvdw,grxc,gsqcut,indsym,&
175 : & maxfor,mgfft,mpi_enreg,n1xccc,n3xccc,&
176 23761 : & nattyp,nfft,ngfft,ngrvdw,ntypat,&
177 23761 : & pawrad,pawtab,ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,&
178 : & usefock,usevxctau,&
179 23761 : & vresid,vxc,vxctau,wvl,wvl_den,xred,&
180 : & electronpositron) ! optional argument
181 :
182 : !Arguments ------------------------------------
183 : !scalars
184 : integer,intent(in) :: mgfft,n1xccc,n3xccc,nfft,ngrvdw,ntypat,usefock,usevxctau
185 : real(dp),intent(in) :: gsqcut
186 : real(dp),intent(out) :: diffor,maxfor
187 : type(MPI_type),intent(in) :: mpi_enreg
188 : type(efield_type),intent(in) :: dtefield
189 : type(dataset_type),intent(in) :: dtset
190 : type(electronpositron_type),pointer,optional :: electronpositron
191 : type(pseudopotential_type),intent(in) :: psps
192 : type(wvl_internal_type), intent(in) :: wvl
193 : type(wvl_denspot_type), intent(inout) :: wvl_den
194 : type(fock_type),pointer, intent(inout) :: fock
195 : !arrays
196 : integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom)
197 : integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym)
198 : real(dp),intent(in) :: grchempottn(3,dtset%natom),grcondft(3,dtset%natom),grewtn(3,dtset%natom)
199 : real(dp),intent(in) :: grvdw(3,ngrvdw),grnl(3*dtset%natom)
200 : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
201 : real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden)
202 : real(dp),intent(in) :: vxc(nfft,dtset%nspden),vxctau(nfft,dtset%nspden,4*usevxctau)
203 : real(dp),intent(inout) :: fcart(3,dtset%natom),forold(3,dtset%natom)
204 : real(dp),intent(inout) :: vresid(nfft,dtset%nspden),xred(3,dtset%natom)
205 : real(dp),intent(out) :: favg(3),gred(3,dtset%natom),gresid(3,dtset%natom)
206 : real(dp),intent(out) :: grhf(3,dtset%natom),rprimd(3,3)
207 : real(dp),intent(inout) :: grxc(3,dtset%natom)
208 : real(dp),intent(out) :: synlgr(3,dtset%natom)
209 : type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw)
210 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
211 :
212 : !Local variables-------------------------------
213 : !scalars
214 : integer :: coredens_method,coretau_method,fdir,has_vxctau,iatom,idir,indx,ipositron,itypat,mu
215 : integer :: optatm,optdyfr,opteltfr,optgr,option,optn,optn2,optstr,optv,vloc_method
216 : real(dp) :: eei_dum1,eei_dum2,ucvol,ucvol_local,vol_element,entropy_dum1
217 : logical :: calc_epaw3_forces, efield_flag
218 : logical :: is_hybrid_ncpp
219 : !arrays
220 : integer :: qprtrb_dum(3)
221 : real(dp) :: dummy6(6),ep3(3),fioncart(3),gmet(3,3),gprimd(3,3)
222 : real(dp) :: rmet(3,3),strn_dummy6(6),strv_dummy6(6),tsec(2),vprtrb_dum(2)
223 23761 : real(dp),allocatable :: atmrho_dum(:),atmvloc_dum(:),dyfrlo_dum(:,:,:)
224 23761 : real(dp),allocatable :: dyfrn_dum(:,:,:),dyfrv_dum(:,:,:)
225 23761 : real(dp),allocatable :: dyfrx2_dum(:,:,:),eltfrn_dum(:,:),gauss_dum(:,:)
226 23761 : real(dp),allocatable :: epawf3red(:,:),fin(:,:),fionred(:,:),grl(:,:),grl_dum(:,:)
227 23761 : real(dp),allocatable :: grnl_tmp(:,:),grtn(:,:),grtn_indx(:,:),grxctau(:,:),v_dum(:),vxctotg(:,:)
228 23761 : real(dp),allocatable :: xccc3d_dum(:)
229 :
230 : ! *************************************************************************
231 :
232 23761 : call timab(69,1,tsec)
233 : ABI_NVTX_START_RANGE(NVTX_FORCES)
234 :
235 : !Save input value of forces
236 71283 : ABI_MALLOC(fin,(3,dtset%natom))
237 286261 : fin(:,:)=fcart(:,:)
238 :
239 : !Compute different geometric tensor, as well as ucvol, from rprimd
240 23761 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
241 :
242 : !Test size of vxctau
243 95044 : has_vxctau=size(vxctau)
244 23761 : if (has_vxctau/=nfft*dtset%nspden*4.and.has_vxctau/=0) then
245 0 : ABI_BUG('Wrong size for vxctau!')
246 : else if (has_vxctau/=0) then
247 23761 : has_vxctau=1
248 : end if
249 :
250 : !Check if we're in hybrid norm conserving pseudopotential
251 : is_hybrid_ncpp=(psps%usepaw==0 .and. &
252 23761 : & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
253 :
254 : !=======================================================================
255 : !========= Local pseudopotential and core charge contributions =========
256 : !=======================================================================
257 :
258 47522 : ABI_MALLOC(grl,(3,dtset%natom))
259 :
260 : !Determine by which method the local ionic potential and/or the pseudo core
261 : ! charge density contributions have to be computed
262 : !Local ionic potential:
263 : ! Method 1: PAW
264 : ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
265 23761 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
266 23761 : if (dtset%icoulomb>0) vloc_method=2
267 23761 : if (psps%usewvl==1) vloc_method=2
268 : !Pseudo core charge density:
269 : ! Method 1 construct through G space FT: PAW, nc_xccc_gspace==1
270 : ! Method 2 construct in real space : Norm-conserving PP default, wavelets
271 : coredens_method=1;if (psps%usepaw==0) coredens_method=2
272 23761 : if (psps%nc_xccc_gspace==1) coredens_method=1
273 23761 : if (psps%nc_xccc_gspace==0) coredens_method=2
274 23761 : if (psps%usewvl==1) coredens_method=2
275 : !Pseudo core charge kinetic energy density: same convention as coredens_method
276 23761 : coretau_method=0
277 23761 : if (dtset%usekden==1) then
278 60 : coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2
279 : end if
280 :
281 : !Local ionic potential and/or pseudo core charge by method 1
282 23761 : if (vloc_method==1.or.coredens_method==1.or.coretau_method==1) then
283 6697 : if (psps%nc_xccc_gspace==1.and.psps%usepaw==0.and.is_hybrid_ncpp) then
284 0 : ABI_BUG(' Not yet implemented !')
285 : end if
286 6697 : call timab(550,1,tsec)
287 : ! Allocate (unused) dummy variables, otherwise some compilers complain
288 6697 : ABI_MALLOC(gauss_dum,(0,0))
289 6697 : ABI_MALLOC(atmrho_dum,(0))
290 6697 : ABI_MALLOC(atmvloc_dum,(0))
291 6697 : ABI_MALLOC(dyfrn_dum,(0,0,0))
292 6697 : ABI_MALLOC(dyfrv_dum,(0,0,0))
293 6697 : ABI_MALLOC(eltfrn_dum,(0,0))
294 : ! Compute Vxc in reciprocal space
295 6697 : if (coredens_method==1.and.n3xccc>0) then
296 12549 : ABI_MALLOC(v_dum,(nfft))
297 12549 : ABI_MALLOC(vxctotg,(2,nfft))
298 157580223 : v_dum(:)=vxc(:,1);if (dtset%nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
299 4183 : call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
300 : call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
301 4183 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
302 4183 : ABI_FREE(v_dum)
303 : else
304 2514 : ABI_MALLOC(vxctotg,(0,0))
305 : end if
306 : ! Compute contribution to forces from Vloc and/or pseudo core density
307 6697 : optv=0;if (vloc_method==1) optv=1
308 6697 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
309 6697 : optatm=0;optdyfr=0;optgr=1;optstr=0;optn2=1;opteltfr=0
310 6697 : if (vloc_method==1.or.coredens_method==1) then
311 : call atm2fft(atindx1,atmrho_dum,atmvloc_dum,dyfrn_dum,dyfrv_dum,&
312 : & eltfrn_dum,gauss_dum,gmet,gprimd,&
313 : & grxc,grl,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
314 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb_dum,&
315 : & dtset%rcut,rhog,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,psps%vlspl,&
316 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
317 6697 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
318 : end if
319 43501 : if (n3xccc==0.and.coredens_method==1) grxc=zero
320 6697 : ABI_FREE(vxctotg)
321 6697 : if (dtset%usekden==1.and.usevxctau==1.and.coretau_method==1.and.n3xccc>0) then
322 : ! Compute contribution to forces from pseudo kinetic energy core density
323 29 : optv=0;optn=1;optn2=4
324 58 : ABI_MALLOC(grxctau,(3,dtset%natom))
325 29 : ABI_MALLOC(grl_dum,(0,0))
326 87 : ABI_MALLOC(v_dum,(nfft))
327 87 : ABI_MALLOC(vxctotg,(2,nfft))
328 1353053 : v_dum(:)=vxctau(:,1,1);if (dtset%nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxctau(:,2,1))
329 29 : call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0)
330 : call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
331 29 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
332 29 : ABI_FREE(v_dum)
333 : call atm2fft(atindx1,atmrho_dum,atmvloc_dum,dyfrn_dum,dyfrv_dum,&
334 : & eltfrn_dum,gauss_dum,gmet,gprimd,&
335 : & grxctau,grl_dum,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
336 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,qprtrb_dum,&
337 : & dtset%rcut,rhog,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxctotg,vxctotg,vxctotg,vprtrb_dum,psps%vlspl,&
338 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
339 29 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
340 261 : grxc(:,:)=grxc(:,:)+grxctau(:,:)
341 29 : ABI_FREE(grl_dum)
342 29 : ABI_FREE(grxctau)
343 29 : ABI_FREE(vxctotg)
344 : end if
345 : ! Deallocate temporary arrays
346 6697 : ABI_FREE(gauss_dum)
347 6697 : ABI_FREE(atmrho_dum)
348 6697 : ABI_FREE(atmvloc_dum)
349 6697 : ABI_FREE(dyfrn_dum)
350 6697 : ABI_FREE(dyfrv_dum)
351 6697 : ABI_FREE(eltfrn_dum)
352 6697 : call timab(550,2,tsec)
353 : end if
354 :
355 : !Local ionic potential by method 2
356 23761 : if (vloc_method==2) then
357 17067 : option=2
358 51201 : ABI_MALLOC(dyfrlo_dum,(3,3,dtset%natom))
359 34134 : ABI_MALLOC(grtn_indx,(3,dtset%natom))
360 51201 : ABI_MALLOC(v_dum,(nfft))
361 : call mklocl(dtset,dyfrlo_dum,eei_dum1,gmet,gprimd,grtn_indx,gsqcut,dummy6,mgfft,&
362 : & mpi_enreg,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,&
363 17067 : & qprtrb_dum,rhog,rhor,rprimd,ucvol,vprtrb_dum,v_dum,wvl,wvl_den,xred)
364 61247 : do iatom=1,dtset%natom
365 : ! Has to use the indexing array atindx1
366 193787 : grl(1:3,atindx1(iatom))=grtn_indx(1:3,iatom)
367 : end do
368 17067 : ABI_FREE(dyfrlo_dum)
369 17067 : ABI_FREE(grtn_indx)
370 17067 : ABI_FREE(v_dum)
371 : ! If gradients are computed in real space, we need to symmetrize the system before summing.
372 : ! Rshaltaf: I changed the following line to include surfaces BC
373 17067 : if (dtset%icoulomb == 1 .or. dtset%icoulomb == 2) then
374 0 : ABI_MALLOC(grnl_tmp,(3,dtset%natom))
375 0 : call sygrad(grnl_tmp,dtset%natom,grl,dtset%nsym,symrec,indsym)
376 0 : grl(:, :) = grnl_tmp(:, :)
377 0 : ABI_FREE(grnl_tmp)
378 : end if
379 : end if
380 :
381 : !Pseudo core electron density by method 2
382 23761 : if (coredens_method==2.or.coretau_method==2) then
383 17064 : if (n1xccc/=0) then
384 4093 : call timab(53,1,tsec)
385 4093 : option=2
386 12279 : ABI_MALLOC(dyfrx2_dum,(3,3,dtset%natom))
387 12279 : ABI_MALLOC(xccc3d_dum,(n3xccc))
388 4093 : if (coredens_method==2) then
389 4093 : if (is_hybrid_ncpp) then
390 : call xchybrid_ncpp_cc(dtset,eei_dum1,entropy_dum1,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
391 100 : & dummy6,eei_dum2,xccc3d_dum,grxc=grxc,xcccrc=psps%xcccrc,xccc1d=psps%xccc1d,xred=xred,n1xccc=n1xccc)
392 : else
393 3993 : if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then
394 : call mkcore(dummy6,dyfrx2_dum,grxc,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
395 : & ngfft(1),n1xccc, ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,vxc,&
396 3993 : & psps%xcccrc,psps%xccc1d,xccc3d_dum,xred)
397 0 : else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then
398 : call mkcore_alt(atindx1,dummy6,dyfrx2_dum,grxc,dtset%icoulomb,mpi_enreg,dtset%natom,nfft,&
399 : & dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
400 0 : & ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d_dum,xred,pawrad,pawtab,psps%usepaw)
401 : else if (psps%usewvl==1.and.psps%usepaw==1) then
402 : ucvol_local=ucvol
403 : #if defined HAVE_BIGDFT
404 : ! ucvol_local=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp)
405 : ! call mkcore_wvl_old(atindx1,dummy6,dyfrx2_dum,wvl%atoms%astruct%geocode,grxc,wvl%h,dtset%natom,&
406 : ! & nattyp,nfft,wvl_den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),dtset%nspden,ntypat,&
407 : ! & wvl%Glr%d%n1,wvl%Glr%d%n1i,wvl%Glr%d%n2,wvl%Glr%d%n2i,wvl%Glr%d%n3,wvl_den%denspot%dpbox%n3pi,&
408 : ! & n3xccc,option,pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d_dum,xred,&
409 : ! & mpi_comm_wvl=mpi_enreg%comm_wvl)
410 : call mkcore_wvl(atindx1,dummy6,grxc,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,&
411 : & n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d_dum,&
412 : & psps%xcccrc,xred,wvl_den,wvl,mpi_comm_wvl=mpi_enreg%comm_wvl)
413 : #endif
414 : end if
415 : end if
416 : end if
417 4093 : if (dtset%usekden==1.and.usevxctau==1.and.coretau_method==2) then
418 6 : ABI_MALLOC(grxctau,(3,dtset%natom))
419 : call mkcore_alt(atindx1,dummy6,dyfrx2_dum,grxctau,dtset%icoulomb,mpi_enreg,dtset%natom,nfft,&
420 : & dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
421 : & ucvol,vxctau(:,:,1),psps%xcccrc,psps%xccc1d,xccc3d_dum,xred,pawrad,pawtab,psps%usepaw,&
422 3 : & usekden=.true.)
423 23 : grxc(:,:)=grxc(:,:)+grxctau(:,:)
424 3 : ABI_FREE(grxctau)
425 : end if
426 4093 : ABI_FREE(xccc3d_dum)
427 4093 : ABI_FREE(dyfrx2_dum)
428 4093 : call timab(53,2,tsec)
429 : else
430 130507 : grxc(:,:)=zero
431 : end if
432 : end if
433 :
434 : !=======================================================================
435 : !===================== Nonlocal contributions ==========================
436 : !=======================================================================
437 :
438 : !Only has to apply symmetries
439 47522 : ABI_MALLOC(grnl_tmp,(3,dtset%natom))
440 89386 : do iatom=1,dtset%natom
441 286261 : indx=3*(iatom-1);grnl_tmp(1:3,atindx1(iatom))=grnl(indx+1:indx+3)
442 : end do
443 23761 : if (dtset%usewvl == 0) then
444 23761 : call sygrad(synlgr,dtset%natom,grnl_tmp,dtset%nsym,symrec,indsym)
445 : else
446 0 : synlgr = grnl_tmp
447 : end if
448 23761 : ABI_FREE(grnl_tmp)
449 :
450 : !=======================================================================
451 : !============ Density/potential residual contributions =================
452 : !=======================================================================
453 :
454 23761 : if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1.and.abs(dtset%densfor_pred)<=3) then
455 : call fresid(dtset,gresid,mpi_enreg,nfft,ngfft,ntypat,1,&
456 19370 : & pawtab,rhor,rprimd,ucvol,vresid,xred,xred,psps%znuclpsp)
457 4391 : else if (dtset%usewvl==0.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) then
458 : call fresidrsp(atindx1,dtset,gmet,gprimd,gresid,gsqcut,mgfft,&
459 : & mpi_enreg,psps%mqgrid_vl,nattyp,nfft,ngfft,ntypat,psps,pawtab,ph1d,&
460 4373 : & psps%qgrid_vl,rprimd,ucvol,psps%usepaw,vresid,psps%zionpsp,psps%znuclpsp)
461 : else
462 470 : gresid(:,:)=zero
463 : end if
464 :
465 : !=======================================================================
466 : !======================= Other contributions ===========================
467 : !=======================================================================
468 :
469 : !Ewald energy contribution to forces as already been computed in "ewald"
470 :
471 : !Potential residual contribution to forces as already been computed (forstr)
472 :
473 : !Add Berry phase contributions (berryopt == 4/6/7/14/16/17)
474 : !(compute the electric field force on the ion cores)
475 : efield_flag = (dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
476 23761 : & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)
477 25 : calc_epaw3_forces = (efield_flag .and. dtset%optforces /= 0 .and. psps%usepaw == 1)
478 : if ( efield_flag ) then
479 50 : ABI_MALLOC(fionred,(3,dtset%natom))
480 449 : fionred(:,:)=zero
481 131 : do iatom=1,dtset%natom
482 106 : itypat=dtset%typat(iatom)
483 : ! force on ion due to electric field, cartesian representation
484 424 : fioncart(:)=psps%ziontypat(itypat)*dtset%efield(:)
485 : ! form fionred = rprimd^T * fioncart, note that forces transform
486 : ! oppositely to coordinates, because they are derivative with respect to
487 : ! coordinates
488 131 : call dgemv('T',3,3,one,rprimd,3,fioncart,1,zero,fionred(1:3,iatom),1)
489 : ! do mu=1,3
490 : ! fionred(mu,iatom)=rprimd(1,mu)*fioncart(1) &
491 : !& +rprimd(2,mu)*fioncart(2) &
492 : !& +rprimd(3,mu)*fioncart(3)
493 : ! end do
494 : end do
495 : end if
496 :
497 : !(compute additional F3-type force due to projectors for electric field with PAW)
498 25 : if ( efield_flag .and. calc_epaw3_forces ) then
499 18 : ABI_MALLOC(epawf3red,(3,dtset%natom))
500 : ! dtefield%epawf3(iatom,idir,fdir) contains
501 81 : epawf3red(:,:)=zero
502 27 : do iatom=1,dtset%natom
503 81 : do fdir = 1, 3
504 216 : do idir = 1, 3
505 : ! vol_element is volume/pt for integration of epawf3. volume is BZ volume
506 : ! so 1/ucvol, and number of kpts is nstr(idir)*nkstr(idir)
507 162 : vol_element=one/(ucvol*dtefield%nstr(idir)*dtefield%nkstr(idir))
508 216 : ep3(idir) = vol_element*dtefield%epawf3(iatom,idir,fdir)
509 : end do
510 234 : epawf3red(fdir,iatom) = -ucvol*dot_product(dtset%red_efieldbar(1:3),ep3(1:3))
511 : end do
512 : end do ! end loop over iatom
513 : end if
514 :
515 : !This was incorrect coding. Bug found by Jiawang Hong
516 : !if (dtset%berryopt==4) then
517 : !allocate(fionred(3,dtset%natom));fionred(:,:)=zero
518 : !iatom = 0
519 : !do itypat=1,ntypat
520 : !do iattyp=1,nattyp(itypat)
521 : !iatom=iatom+1
522 : !fioncart(:)=psps%ziontypat(itypat)*dtset%efield(:)
523 : !do mu=1,3
524 : !fionred(mu,iatom)=rprimd(1,mu)*fioncart(1) &
525 : !& +rprimd(2,mu)*fioncart(2) &
526 : !& +rprimd(3,mu)*fioncart(3)
527 : !end do
528 : !end do
529 : !end do
530 : !end if
531 :
532 : !=======================================================================
533 : !======= Assemble the various contributions to the forces ==============
534 : !=======================================================================
535 :
536 : !Collect grads of etot wrt reduced coordinates
537 : !This gives non-symmetrized Hellman-Feynman reduced gradients
538 47522 : ABI_MALLOC(grtn,(3,dtset%natom))
539 286261 : grtn(:,:)=grl(:,:)+grchempottn(:,:)+grcondft(:,:)+grewtn(:,:)+synlgr(:,:)+grxc(:,:)
540 :
541 23761 : if (usefock==1 .and. associated(fock)) then
542 320 : if (fock%fock_common%optfor) then
543 2070 : grtn(:,:)=grtn(:,:)+fock%fock_common%forces(:,:)
544 : end if
545 : end if
546 :
547 23833 : if (ngrvdw==dtset%natom) grtn(:,:)=grtn(:,:)+grvdw(:,:)
548 : ! note that fionred is subtracted, because it really is a force and we need to
549 : ! turn it back into a gradient. The gred2fcart routine below includes the minus
550 : ! sign to convert gradients back to forces
551 24185 : if ( efield_flag ) grtn(:,:)=grtn(:,:)-fionred(:,:)
552 : ! epawf3red is added, because it actually is a gradient, not a force
553 97 : if ( efield_flag .and. calc_epaw3_forces ) grtn(:,:) = grtn(:,:) + epawf3red(:,:)
554 :
555 : !Symmetrize explicitly for given space group and store in grhf :
556 23761 : call sygrad(grhf,dtset%natom,grtn,dtset%nsym,symrec,indsym)
557 :
558 : !If residual forces are too large, there must be a problem: cancel them !
559 23761 : if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>0.and.abs(dtset%densfor_pred)/=5) then
560 89255 : do iatom=1,dtset%natom
561 285791 : do mu=1,3
562 262048 : if (abs(gresid(mu,iatom))>10000._dp*abs(grtn(mu,iatom))) gresid(mu,iatom)=zero
563 : end do
564 : end do
565 : end if
566 :
567 : !Add residual potential correction
568 286261 : grtn(:,:)=grtn(:,:)+gresid(:,:)
569 :
570 : !Additional stuff for electron-positron
571 23761 : ipositron=0
572 23761 : if (present(electronpositron)) then
573 23761 : if (associated(electronpositron)) then
574 289 : if (allocated(electronpositron%gred_ep)) ipositron=electronpositron_calctype(electronpositron)
575 : end if
576 : end if
577 289 : if (abs(ipositron)==1) then
578 5440 : grtn(:,:)=grtn(:,:)-grxc(:,:)-grchempottn(:,:)-grcondft(:,:)-grewtn(:,:)-gresid(:,:)-two*grl(:,:)
579 : ! grtn(:,:)=grtn(:,:)-grxc(:,:)-grewtn(:,:)-gresid(:,:)-two*grl(:,:)
580 16120 : grl(:,:)=-grl(:,:);grxc(:,:)=zero;gresid(:,:)=zero
581 100 : if (ngrvdw==dtset%natom) grtn(:,:)=grtn(:,:)-grvdw(:,:)
582 : if ( dtset%berryopt== 4 .or. dtset%berryopt== 6 .or. dtset%berryopt== 7 .or. &
583 : & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) then
584 0 : grtn(:,:)=grtn(:,:)+fionred(:,:)
585 0 : fionred(:,:)=zero
586 : end if
587 : end if
588 39009 : if (ipositron>0) grtn(:,:)=grtn(:,:)+electronpositron%gred_ep(:,:)
589 :
590 : !Symmetrize all grads explicitly for given space group:
591 23761 : if (dtset%usewvl == 0) then
592 23761 : call sygrad(gred,dtset%natom,grtn,dtset%nsym,symrec,indsym)
593 : else
594 0 : gred = grtn
595 : end if
596 :
597 : !Conversion to cartesian coordinates (bohr) AND
598 : !Subtract off average force from each force component
599 : !to avoid spurious drifting of atoms across cell.
600 : ! notice that gred2fcart multiplies gred by -1 to convert it
601 : ! from a gradient (input) to a force (output)
602 :
603 23953 : call gred2fcart(favg,(dtset%jellslab==0 .and. dtset%nzchempot==0),fcart,gred,gprimd,dtset%natom)
604 :
605 : !Compute maximal force and maximal difference
606 23761 : maxfor=zero;diffor=zero
607 89386 : do iatom=1,dtset%natom
608 286261 : do mu=1,3
609 262500 : if (dtset%iatfix(mu,iatom) /= 1) then
610 165847 : maxfor=max(maxfor,abs(fcart(mu,iatom)))
611 165847 : diffor=max(diffor,abs(fcart(mu,iatom)-fin(mu,iatom)))
612 31028 : else if (dtset%ionmov==4 .or. dtset%ionmov==5) then
613 : ! Make the force vanish on fixed atoms when ionmov=4 or 5
614 : ! This is because fixing of atom cannot be imposed at the
615 : ! level of a routine similar to brdmin or moldyn for these options.
616 672 : fcart(mu,iatom)=zero
617 : end if
618 : end do
619 : end do
620 :
621 : !Apply any generalized constraints to the forces
622 23761 : if (dtset%nconeq>0) call constrf(diffor,fcart,forold,gred,dtset%iatfix,dtset%ionmov,maxfor,&
623 740 : & dtset%natom,dtset%nconeq,dtset%prtvol,rprimd,dtset%wtatcon,xred)
624 :
625 : !=======================================================================
626 : !Memory deallocations
627 23761 : ABI_FREE(grl)
628 23761 : ABI_FREE(grtn)
629 23761 : ABI_FREE(fin)
630 23761 : if ( efield_flag ) then
631 25 : ABI_FREE(fionred)
632 25 : if ( calc_epaw3_forces ) then
633 9 : ABI_FREE(epawf3red)
634 : end if
635 : end if
636 :
637 : ABI_NVTX_END_RANGE()
638 23761 : call timab(69,2,tsec)
639 :
640 23761 : end subroutine forces
641 : !!***
642 :
643 : !!****f* ABINIT/sygrad
644 : !!
645 : !! NAME
646 : !! sygrad
647 : !!
648 : !! FUNCTION
649 : !! Symmetrize derivatives of energy with respect to coordinates.
650 : !! Unsymmetrized gradients are input as dedt; symmetrized grads are then placed in gred.
651 : !! If nsym=1 simply copy dedt into gred (only symmetry is identity).
652 : !!
653 : !! INPUTS
654 : !! natom=number of atoms in cell
655 : !! dedt(3,natom)=unsymmetrized gradients wrt dimensionless tn (hartree)
656 : !! nsym=number of symmetry operators in group
657 : !! symrec(3,3,nsym)=symmetries of group in terms of operations on
658 : !! reciprocal space primitive translations--see comments below
659 : !! indsym(4,nsym,natom)=label given by subroutine symatm, indicating atom
660 : !! label which gets rotated into given atom by given symmetry
661 : !! (first three elements are related primitive translation--
662 : !! see symatm where this is computed)
663 : !!
664 : !! OUTPUT
665 : !! gred(3,3,natom)=symmetrized gradients wrt reduced coordinates (hartree)
666 : !!
667 : !! NOTES
668 : !! symmetrization of gradients with respect to reduced
669 : !! coordinates tn is conducted according to the expression
670 : !! $[d(e)/d(t(n,a))]_{symmetrized} = (1/Nsym)*Sum(S)*symrec(n,m,S)*
671 : !! [d(e)/d(t(m,b))]_{unsymmetrized}$
672 : !! where $t(m,b)= (symrel^{-1})(m,n)*(t(n,a)-tnons(n))$ and tnons
673 : !! is a possible nonsymmorphic translation. The label "b" here
674 : !! refers to the atom which gets rotated into "a" under symmetry "S".
675 : !! symrel is the symmetry matrix in real space, which is the inverse
676 : !! transpose of symrec. symrec is the symmetry matrix in reciprocal
677 : !! space. $sym_{cartesian} = R * symrel * R^{-1} = G * symrec * G^{-1}$
678 : !! where the columns of R and G are the dimensional primitive translations
679 : !! in real and reciprocal space respectively.
680 : !! Note the use of "symrec" in the symmetrization expression above.
681 : !!
682 : !! SOURCE
683 :
684 71283 : subroutine sygrad(gred,natom,dedt,nsym,symrec,indsym)
685 :
686 : !Arguments ------------------------------------
687 : !scalars
688 : integer,intent(in) :: natom,nsym
689 : !arrays
690 : integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym)
691 : real(dp),intent(in) :: dedt(3,natom)
692 : real(dp),intent(out) :: gred(3,natom)
693 :
694 : !Local variables-------------------------------
695 : !scalars
696 : integer :: ia,ind,isym,mu
697 : real(dp),parameter :: tol=1.0d-30
698 : real(dp) :: summ
699 :
700 : ! *************************************************************************
701 : !
702 71283 : if (nsym==1) then
703 : ! only symmetry is identity so simply copy
704 332601 : gred(:,:)=dedt(:,:)
705 : else
706 : ! actually conduct symmetrization
707 168972 : do ia=1,natom
708 526182 : do mu=1,3
709 : summ=0._dp
710 4325130 : do isym=1,nsym
711 3967920 : ind=indsym(4,isym,ia)
712 : summ=summ+dble(symrec(mu,1,isym))*dedt(1,ind)+&
713 : & dble(symrec(mu,2,isym))*dedt(2,ind)+&
714 4325130 : & dble(symrec(mu,3,isym))*dedt(3,ind)
715 : end do
716 357210 : gred(mu,ia)=summ/dble(nsym)
717 476280 : if(abs(gred(mu,ia))<tol)gred(mu,ia)=0.0_dp
718 : end do
719 : end do
720 : end if
721 :
722 71283 : end subroutine sygrad
723 : !!***
724 :
725 : !!****f* ABINIT/fresidrsp
726 : !!
727 : !! NAME
728 : !! fresidrsp
729 : !!
730 : !! FUNCTION
731 : !! Compute the forces due to the residual of the potential (or density)
732 : !! in RECIPROCAL SPACE, using
733 : !! - the atomic density read in psp file (PAW or NC with nctval_spl e.g. psp8 format)
734 : !! - a gaussian atomic density (norm-conserving psps if nctval_spl is not available)
735 : !!
736 : !! INPUTS
737 : !! atindx1(natom)=index table for atoms, inverse of atindx
738 : !! dtset <type(dataset_type)>=all input variables in this dataset
739 : !! | densty(ntypat,4)=parameters for initialisation of the density of each atom type
740 : !! | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver
741 : !! | ixc= choice of exchange-correlation scheme
742 : !! | natom=number of atoms in cell.
743 : !! | nspden=number of spin-density components
744 : !! | typat(natom)=integer type for each atom in cell
745 : !! gmet(3,3)=reciprocal space metric
746 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
747 : !! gsqcut=cutoff value on G**2 for sphere inside fft box
748 : !! mgfft=maximum size of 1D FFTs
749 : !! mpi_enreg=information about MPI parallelization
750 : !! mqgrid=number of grid pts in q array for atomic density spline n^AT(q)
751 : !! nattyp(ntypat)=number of atoms of each type in cell
752 : !! nfft=(effective) number of FFT grid points (for this processor)
753 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
754 : !! ntypat=number of types of atoms in cell.
755 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
756 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
757 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase information for given atom coordinates.
758 : !! qgrid(mqgrid)=q grid for spline atomic valence density n^AT(q) from 0 to qmax
759 : !! ucvol=unit cell volume (bohr**3).
760 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
761 : !! vresid(nfft,nspden)=potential residual - (non-collinear magn. : only V11 and V22 are used)
762 : !! zion(ntypat)=charge on each type of atom (real number)
763 : !! znucl(ntypat)=atomic number, for each type of atom
764 : !!
765 : !! OUTPUT
766 : !! gresid(3,natom)=forces due to the residual of the potential
767 : !!
768 : !! SOURCE
769 :
770 4373 : subroutine fresidrsp(atindx1,dtset,gmet,gprimd,gresid,gsqcut,mgfft,mpi_enreg,mqgrid,nattyp,nfft,&
771 4373 : & ngfft,ntypat,psps,pawtab,ph1d,qgrid,rprimd,ucvol,usepaw,vresid,zion,znucl)
772 :
773 : !Arguments ------------------------------------
774 : !scalars
775 : integer,intent(in) :: mgfft,mqgrid,nfft,ntypat,usepaw
776 : real(dp),intent(in) :: gsqcut,ucvol
777 : type(pseudopotential_type),intent(in) :: psps
778 : type(MPI_type),intent(in) :: mpi_enreg
779 : type(dataset_type),intent(in) :: dtset
780 : !arrays
781 : integer,intent(in) :: atindx1(dtset%natom),nattyp(ntypat),ngfft(18)
782 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*dtset%natom)
783 : real(dp),intent(in) :: qgrid(mqgrid),vresid(nfft,dtset%nspden),zion(ntypat)
784 : real(dp),intent(inout) :: rprimd(3,3)
785 : real(dp),intent(in) :: znucl(ntypat)
786 : real(dp),intent(out) :: gresid(3,dtset%natom)
787 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
788 :
789 : !Local variables-------------------------------
790 : !scalars
791 : integer :: itypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv
792 : logical :: usegauss
793 : !arrays
794 : integer :: dummy3(3)
795 : real(dp) :: dummy2(2)
796 : real(dp) :: dummy_in1(0),dummy_in2(0)
797 : real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
798 : real(dp) :: strn_dummy6(6),strv_dummy6(6)
799 4373 : real(dp),allocatable :: gauss(:,:),vresg(:,:),work(:)
800 :
801 : ! *************************************************************************
802 :
803 : !Inits
804 4373 : optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=0;optn=1
805 13119 : ABI_MALLOC(vresg,(2,nfft))
806 :
807 : !Transfer potential residual to reciprocal space
808 : !Use only Vres=Vres11+Vres22=Vres_up+Vres_dn
809 13119 : ABI_MALLOC(work,(nfft))
810 231692637 : work(:)=vresid(:,1)
811 1192481 : if (dtset%nspden>=2) work(:)=work(:)+vresid(:,2)
812 4373 : call fourdp(1,vresg,work,-1,mpi_enreg,nfft,1,ngfft,0)
813 4373 : ABI_FREE(work)
814 :
815 : !Determine whether a gaussan atomic density has to be used or not
816 4373 : usegauss=.true.
817 5071 : if (usepaw==0) usegauss = any(.not.psps%nctab(1:ntypat)%has_tvale)
818 9988 : if (usepaw==1) usegauss=(minval(pawtab(1:ntypat)%has_tvale)==0)
819 4373 : if (usegauss) then
820 1837 : optn2=3
821 5511 : ABI_MALLOC(gauss,(2,ntypat))
822 3686 : do itypat=1,ntypat
823 1849 : gauss(1,itypat)=zion(itypat)
824 3686 : gauss(2,itypat) = atom_length(dtset%densty(itypat,1),zion(itypat),znucl(itypat))
825 : end do
826 1837 : call wrtout(std_out," Computing residual forces using gaussian functions as atomic densities", "COLL")
827 : else
828 2536 : optn2=2
829 2536 : ABI_MALLOC(gauss,(2,0))
830 2536 : call wrtout(std_out," Computing residual forces using atomic densities taken from pseudos", "COLL")
831 : end if
832 :
833 : !Compute forces due to residual
834 : call atm2fft(atindx1,dummy_out1,dummy_out2,dummy_out3,dummy_out4,&
835 : & dummy_out5,gauss,gmet,gprimd,gresid,dummy_out6,gsqcut,mgfft,&
836 : & mqgrid,dtset%natom,nattyp,nfft,ngfft,ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
837 : & psps,pawtab,ph1d,qgrid,dummy3,dtset%rcut,dummy_in1,rprimd,strn_dummy6,strv_dummy6,ucvol,usepaw,&
838 : & vresg,vresg,vresg,dummy2,dummy_in2,comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
839 4373 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
840 :
841 : !In case of nspden>=2, has to apply 1/2 factor
842 4701 : if (dtset%nspden>=2) gresid=gresid*half
843 :
844 4373 : ABI_FREE(gauss)
845 4373 : ABI_FREE(vresg)
846 :
847 4373 : end subroutine fresidrsp
848 : !!***
849 :
850 : !!****f* ABINIT/fresid
851 : !! NAME
852 : !! fresid
853 : !!
854 : !! FUNCTION
855 : !! If option=1, compute the forces due to the residual of the potential
856 : !! If option=2, generate approximate new density from old one,
857 : !! old atomic positions, and new atomic positions
858 : !!
859 : !! INPUTS
860 : !! dtset <type(dataset_type)>=all input variables in this dataset
861 : !! | icoulomb=0 periodic treatment of Hartree potential, 1 use of Poisson solver
862 : !! | natom=number of atoms in cell.
863 : !! | nspden=number of spin-density components
864 : !! | typat(natom)=integer type for each atom in cell
865 : !! | usepaw= 0 for non paw calculation; =1 for paw calculation
866 : !! | xclevel= level of the XC functional
867 : !! mpi_enreg=information about MPI parallelization
868 : !! nfft=(effective) number of FFT grid points (for this processor)
869 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
870 : !! ntypat=number of types of atoms in cell.
871 : !! option=see below
872 : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
873 : !! rhor(nfft,nspden)=electron density in electrons/bohr**3 (slices of it if FTT parallelism).
874 : !! rprimd(3,3)=dimensional primitive translation vectors (bohr)
875 : !! ucvol=unit cell volume (bohr**3).
876 : !! xred_new(3,natom)=new reduced coordinates for atoms in unit cell
877 : !! xred_old(3,natom)=old reduced coordinates for atoms in unit cell
878 : !! znucl(ntypat)=real(dp), atomic number of atom type
879 : !!
880 : !! OUTPUT
881 : !! gresid(3,natom)=forces due to the residual of the potential
882 : !!
883 : !! SIDE EFFECTS
884 : !! work(nfft,nspden)=functions on the fft grid (slices of it if FTT parallelism):
885 : !! if option==1, the POTENTIAL residual is input
886 : !! if option==2, the interpolated density is output
887 : !!
888 : !! NOTES
889 : !! FFT parallelism:
890 : !! At the beginning of this routine, the plane-waves are ditributed over all the processors.
891 : !! In the main part, all the processors perform the same calculations over the whole FFT grid.
892 : !! At the end, each processor gets its part of the whole FFT grid.
893 : !! These modifications are not efficient when large FFT grids are used.
894 : !! So they have to be considered as a first step before a comprehensive parallelization of this routine.
895 : !!
896 : !! SOURCE
897 :
898 21387 : subroutine fresid(dtset,gresid,mpi_enreg,nfft,ngfft,ntypat,option,&
899 21387 : & pawtab,rhor,rprimd,ucvol,work,xred_new,xred_old,znucl)
900 :
901 : !Arguments ------------------------------------
902 : !scalars
903 : integer,intent(in) :: nfft,ntypat,option
904 : real(dp),intent(in) :: ucvol
905 : type(MPI_type),intent(in) :: mpi_enreg
906 : type(dataset_type),intent(in) :: dtset
907 : !arrays
908 : integer,intent(in) :: ngfft(18)
909 : real(dp),intent(in) :: rhor(nfft,dtset%nspden),rprimd(3,3)
910 : real(dp),intent(in) :: xred_new(3,dtset%natom),xred_old(3,dtset%natom)
911 : real(dp),intent(in) :: znucl(ntypat)
912 : real(dp),intent(inout) :: work(nfft,dtset%nspden)
913 : real(dp),intent(out) :: gresid(3,dtset%natom)
914 : type(pawtab_type),intent(in) :: pawtab(ntypat*dtset%usepaw)
915 :
916 : !Local variables-------------------------------
917 : !real(dp), parameter :: app_remain=0.001_dp
918 : !scalars
919 : integer,parameter :: natnum=110
920 : integer :: atmove,i1,i1_new,i1m,i1p,i2,i2_new,i2m,i2p,i3,i3_new,i3m,i3p
921 : integer :: iatom,ifft,ifft_new,iloop,ind2m,ind2m3m,ind2p,ind2p3p,ind3m,ind3p
922 : integer :: index,index_new,ishift,ishift1,ishift2,ishift3,ispden,ixp,mshift,mu
923 : integer :: n1,n2,n3,n4,nfft_tmp,nfftot,nu,quit
924 : real(dp),parameter :: app_remain=0.01_dp
925 : real(dp) :: diff_rem1,diff_rem2,diff_rem3,difmag,difmag2
926 : real(dp) :: difmag2_fact,difmag2_part,drho1,drho11,drho12,drho13,drho14
927 : real(dp) :: drho1dn,drho1mx,drho1my,drho1mz,drho1tot,drho1up,drho2,drho21
928 : real(dp) :: drho22,drho23,drho24,drho2dn,drho2mx,drho2my,drho2mz,drho2tot
929 : real(dp) :: drho2up,drho3,drho31,drho32,drho33,drho34,drho3dn,drho3mx,drho3my
930 : real(dp) :: drho3mz,drho3tot,drho3up,drhox00,drhox01,drhox10,drhox11,drhoxy0
931 : real(dp) :: drhoxy1,drhoxyz,fact,range,range2,rcov,rcov2,rcovm1,rdiff1
932 : real(dp) :: rdiff2,rdiff3,vresid1,vresid2,vresid3,vresid4,xx
933 : type(atomdata_t) :: atom
934 : !arrays
935 : integer :: diff_igrid(3),igrid(3),irange(3)
936 21387 : integer,allocatable :: ii(:,:)
937 : real(dp) :: diff_grid(3),diff_rem(3),diff_tau(3),diff_xred(3),lencp(3)
938 : real(dp) :: rho_tot(4),rhosum(4),rmet(3,3),scale(3),tau(3)
939 21387 : real(dp),allocatable :: approp(:),atmrho(:,:),rhor_tot(:,:),rrdiff(:,:)
940 21387 : real(dp),allocatable :: work_tot(:,:)
941 21387 : logical,allocatable :: my_sphere(:)
942 :
943 : ! *************************************************************************
944 :
945 : !Compute lengths of cross products for pairs of primitive
946 : !translation vectors (used in setting index search range below)
947 : lencp(1)=cross_fr(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
948 21387 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
949 : lencp(2)=cross_fr(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
950 21387 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
951 : lencp(3)=cross_fr(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
952 21387 : & rprimd(1,2),rprimd(2,2),rprimd(3,2))
953 :
954 : !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
955 : !(recall ucvol=R1.(R2xR3))
956 85548 : scale(:)=ucvol/lencp(:)
957 :
958 : !initialize diff_igrid, otherwise valgrind complains
959 21387 : diff_igrid=0
960 :
961 : !Compute metric tensor in real space rmet
962 85548 : do nu=1,3
963 278031 : rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
964 : end do
965 :
966 : !FFT parallelization: Starting from now, calculations are performed on the whole FFT grid
967 : !and no more on slices. The nfft variable becomes nfft_tmp until the end
968 21387 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
969 21387 : n4=n3/mpi_enreg%nproc_fft
970 85548 : nfftot=PRODUCT(ngfft(1:3));nfft_tmp=nfftot
971 21387 : if(mpi_enreg%paral_kgb==1) then
972 0 : ABI_MALLOC(rhor_tot,(nfftot,dtset%nspden))
973 0 : ABI_MALLOC(work_tot,(nfftot,dtset%nspden))
974 0 : do ispden=1,dtset%nspden
975 0 : call pre_gather(rhor(:,ispden),rhor_tot(:,ispden),n1,n2,n3,n4,mpi_enreg)
976 0 : call pre_gather(work(:,ispden),work_tot(:,ispden),n1,n2,n3,n4,mpi_enreg)
977 : end do
978 : end if
979 :
980 246991 : gresid(1:3,1:dtset%natom)=0.0_dp
981 21387 : quit=0
982 :
983 : !Initialize appropriation function
984 64161 : ABI_MALLOC(approp,(nfft_tmp))
985 85548 : ABI_MALLOC(atmrho,(nfft_tmp,dtset%nspden))
986 64161 : ABI_MALLOC(my_sphere,(nfft_tmp))
987 :
988 382952698 : approp(:)=app_remain
989 : !First loop over atoms in unit cell : build appropriation function
990 : !Second loop : compute forces
991 64161 : do iloop=1,2
992 :
993 : ! Take into account the remaining density
994 42774 : if(option==2 .and. iloop==2)then
995 2017 : if(mpi_enreg%paral_kgb==1) then
996 : ! FFT parallelization: All the processors perform the same calculation.
997 : ! We divided by nproc_fft in order to "remove" the xmpi_sum made in mean_fftr
998 0 : do ispden=1,dtset%nspden
999 0 : do ifft=1,nfft_tmp
1000 0 : work_tot(ifft,ispden)=rhor_tot(ifft,ispden)*approp(ifft)*app_remain
1001 : end do
1002 : end do
1003 0 : call mean_fftr(work_tot,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1004 0 : rhosum(1:dtset%nspden)=rhosum(1:dtset%nspden)/mpi_enreg%nproc_fft
1005 : else
1006 4209 : do ispden=1,dtset%nspden
1007 38938549 : do ifft=1,nfft_tmp
1008 38936532 : work(ifft,ispden)=rhor(ifft,ispden)*approp(ifft)*app_remain
1009 : end do
1010 : end do
1011 2017 : call mean_fftr(work,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1012 : end if
1013 :
1014 : ! This will be used to restore proper normalization of density
1015 4209 : rho_tot(1:dtset%nspden)=rhosum(1:dtset%nspden)*nfftot
1016 : end if
1017 :
1018 155576 : do iatom=1,dtset%natom
1019 :
1020 : ! Get the covalent radius
1021 112802 : call atomdata_from_znucl(atom,znucl(dtset%typat(iatom)))
1022 112802 : rcov = atom%rcov
1023 : ! PAW choose PAW radius instead...
1024 112802 : if (dtset%usepaw==1) rcov=max(rcov,pawtab(dtset%typat(iatom))%rpaw)
1025 :
1026 : ! Set search range
1027 112802 : rcov2=rcov**2
1028 112802 : range=2._dp*rcov
1029 112802 : range2=range**2
1030 112802 : rcovm1=1.0_dp/rcov
1031 :
1032 : ! Use range to compute an index range along R(1:3)
1033 : ! (add 1 to make sure it covers full range)
1034 112802 : irange(1)=1+nint((range/scale(1))*dble(n1))
1035 112802 : irange(2)=1+nint((range/scale(2))*dble(n2))
1036 112802 : irange(3)=1+nint((range/scale(3))*dble(n3))
1037 :
1038 : ! Allocate ii and rrdiff
1039 451208 : mshift=2*maxval(irange(1:3))+1
1040 338406 : ABI_MALLOC(ii,(mshift,3))
1041 338406 : ABI_MALLOC(rrdiff,(mshift,3))
1042 :
1043 : ! Consider each component in turn
1044 451208 : do mu=1,3
1045 :
1046 : ! Convert reduced coord of given atom to [0,1)
1047 338406 : tau(mu)=mod(xred_old(mu,iatom)+1._dp-aint(xred_old(mu,iatom)),1._dp)
1048 :
1049 : ! Use tau to find nearest grid point along R(mu)
1050 : ! (igrid=0 is the origin; shift by 1 to agree with usual index)
1051 338406 : igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
1052 :
1053 : ! Set up a counter that explore the relevant range
1054 : ! of points around the atom
1055 338406 : ishift=0
1056 7114704 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
1057 6776298 : ishift=ishift+1
1058 6776298 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
1059 7114704 : rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
1060 : end do
1061 :
1062 : ! If option 2, set up quantities related with the change of atomic coordinates
1063 451208 : if(option==2 .and. iloop==2)then
1064 17055 : diff_xred(mu)=xred_new(mu,iatom)-xred_old(mu,iatom)
1065 : ! Convert to [0,1)
1066 17055 : diff_tau(mu)=mod(diff_xred(mu)+1._dp-aint(diff_xred(mu)),1._dp)
1067 : ! Convert to [0,ngfft)
1068 17055 : diff_grid(mu)=diff_tau(mu)*dble(ngfft(mu))
1069 : ! Integer part
1070 17055 : diff_igrid(mu)=int(diff_grid(mu))
1071 : ! Compute remainder
1072 17055 : diff_rem(mu)=diff_grid(mu)-diff_igrid(mu)
1073 :
1074 : ! DEBUG
1075 : ! write(std_out,*)' mu,diff',mu,diff_igrid(mu),diff_rem(mu)
1076 : ! ENDDEBUG
1077 :
1078 : end if
1079 :
1080 : ! End loop on mu
1081 : end do
1082 :
1083 : ! May be the atom is fixed
1084 112802 : atmove=1
1085 112802 : if(option==2 .and. iloop==2)then
1086 5685 : if(diff_xred(1)**2+diff_xred(2)**2+diff_xred(3)**2 < 1.0d-24)then
1087 : atmove=0
1088 : else
1089 4233 : diff_rem1=diff_rem(1)
1090 4233 : diff_rem2=diff_rem(2)
1091 4233 : diff_rem3=diff_rem(3)
1092 : end if
1093 : end if
1094 :
1095 : ! If second loop, initialize atomic density, and the variable
1096 : ! that says whether a fft point belongs to the sphere of the atom
1097 107117 : if(iloop==2) then
1098 1708975250 : atmrho(:,:)=0.0_dp
1099 1596391643 : my_sphere(:)=.false.
1100 : end if
1101 :
1102 : ! Conduct triple loop over restricted range of grid points for iatom
1103 :
1104 2356636 : do ishift3=1,1+2*irange(3)
1105 : ! map back to [1,ngfft(3)] for usual fortran index in unit cell
1106 2243834 : i3=ii(ishift3,3)
1107 2243834 : i3m=i3-1 ; if(i3==1)i3m=n3
1108 2243834 : i3p=i3+1 ; if(i3==n3)i3p=1
1109 :
1110 : ! find vector from atom location to grid point (reduced)
1111 2243834 : rdiff3=rrdiff(ishift3,3)
1112 :
1113 61458442 : do ishift2=1,1+2*irange(2)
1114 59101806 : i2=ii(ishift2,2)
1115 59101806 : i2m=i2-1 ; if(i2==1)i2m=n2
1116 59101806 : i2p=i2+1 ; if(i2==n2)i2p=1
1117 59101806 : index=n1*(i2-1+n2*(i3-1))
1118 59101806 : ind3m=n1*(i2-1+n2*(i3m-1))
1119 59101806 : ind3p=n1*(i2-1+n2*(i3p-1))
1120 59101806 : ind2m=n1*(i2m-1+n2*(i3-1))
1121 59101806 : ind2p=n1*(i2p-1+n2*(i3-1))
1122 59101806 : ind2p3p=n1*(i2p-1+n2*(i3p-1))
1123 :
1124 59101806 : rdiff2=rrdiff(ishift2,2)
1125 : ! Prepare the computation of difmag2
1126 : difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
1127 59101806 : & +2.0_dp*rmet(3,2)*rdiff3*rdiff2
1128 59101806 : difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
1129 :
1130 2069358086 : do ishift1=1,1+2*irange(1)
1131 2008012446 : rdiff1=rrdiff(ishift1,1)
1132 :
1133 : ! Compute (rgrid-tau-Rprim)**2
1134 2008012446 : difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
1135 :
1136 : ! Only accept contribution inside defined range
1137 : ! This condition means that x, calculated below, cannot exceed 2.0_dp
1138 2067114252 : if (difmag2<range2) then
1139 :
1140 : ! Will compute contribution to appropriation function based on
1141 : ! rcov2, range2 and difmag2
1142 677755748 : i1=ii(ishift1,1)
1143 677755748 : ifft=i1+index
1144 :
1145 677755748 : if(iloop==1)then
1146 :
1147 : ! Build appropriation function
1148 338877874 : if (difmag2<rcov2)then
1149 42373024 : approp(ifft)=approp(ifft)+1.0_dp
1150 : else
1151 296504850 : difmag=sqrt(difmag2)
1152 296504850 : xx=difmag*rcovm1
1153 : ! The following function is 1. at xx=1, 0. at xx=2, with vanishing
1154 : ! derivatives at these points.
1155 296504850 : approp(ifft)=approp(ifft)+((2.0_dp*xx-9.0_dp)*xx+12.0_dp)*xx-4.0_dp
1156 : end if
1157 :
1158 : else
1159 :
1160 338877874 : if (difmag2<rcov2) then
1161 : fact=one
1162 : else
1163 296504850 : difmag=sqrt(difmag2)
1164 296504850 : xx=difmag*rcovm1
1165 296504850 : fact=((2.0_dp*xx-9.0_dp)*xx+12.0_dp)*xx-4.0_dp
1166 : end if
1167 :
1168 : ! Build atomic density
1169 338877874 : if(mpi_enreg%paral_kgb==1) then
1170 : atmrho(ifft,1:dtset%nspden)=atmrho(ifft,1:dtset%nspden) &
1171 0 : & +rhor_tot(ifft,1:dtset%nspden)*fact*approp(ifft)
1172 : else
1173 : atmrho(ifft,1:dtset%nspden)=atmrho(ifft,1:dtset%nspden) &
1174 731581060 : & +rhor(ifft,1:dtset%nspden)*fact*approp(ifft)
1175 : end if
1176 :
1177 : ! Compute the sphere of the atom : it is different for
1178 : ! option 1 and for option 2
1179 338877874 : i1p=i1+1 ; if(i1==n1)i1p=1
1180 338877874 : if(option==1)then
1181 307499042 : i1m=i1-1 ; if(i1==1)i1m=n1
1182 307499042 : my_sphere(ifft)=.true.
1183 307499042 : my_sphere(i1p+index)=.true. ; my_sphere(i1m+index)=.true.
1184 307499042 : my_sphere(i1+ind2p)=.true. ; my_sphere(i1+ind2m)=.true.
1185 307499042 : my_sphere(i1+ind3p)=.true. ; my_sphere(i1+ind3m)=.true.
1186 : else
1187 31378832 : my_sphere(ifft)=.true. ; my_sphere(i1p+index)=.true.
1188 31378832 : my_sphere(i1+ind2p)=.true. ; my_sphere(i1p+ind2p)=.true.
1189 31378832 : my_sphere(i1+ind3p)=.true. ; my_sphere(i1p+ind3p)=.true.
1190 31378832 : my_sphere(i1+ind2p3p)=.true. ; my_sphere(i1p+ind2p3p)=.true.
1191 : end if
1192 :
1193 : end if
1194 :
1195 : ! End of condition on the range
1196 : end if
1197 :
1198 : ! End loop on ishift1
1199 : end do
1200 :
1201 : ! End loop on ishift2
1202 : end do
1203 :
1204 : ! End loop on ishift3
1205 : end do
1206 : ! At the end of the second loop for each atom, compute the force
1207 : ! from the atomic densities, or translate density.
1208 : ! In the first case, use a two-point finite-difference approximation,
1209 : ! since this calculation serves only to decrease the error,
1210 : ! and should not be very accurate, but fast.
1211 : ! In the second case, using a crude trilinear interpolation scheme
1212 : ! for the same reason.
1213 : !
1214 : ! The section is skipped if option==2 and the atom is fixed
1215 112802 : if(iloop==2 .and. (option==1 .or. atmove==1) )then
1216 :
1217 1371341 : do i3=1,n3
1218 1316392 : i3m=i3-1 ; if(i3==1)i3m=n3
1219 1316392 : i3p=i3+1 ; if(i3==n3)i3p=1
1220 : ! note: diff_igrid is only set if(option==2 .and. iloop==2)
1221 1316392 : i3_new=i3+diff_igrid(3) ; if(i3_new > n3)i3_new=i3_new-n3
1222 41297387 : do i2=1,n2
1223 39926046 : i2m=i2-1 ; if(i2==1)i2m=n2
1224 39926046 : i2p=i2+1 ; if(i2==n2)i2p=1
1225 39926046 : i2_new=i2+diff_igrid(2) ; if(i2_new > n2)i2_new=i2_new-n2
1226 39926046 : index=n1*(i2-1+n2*(i3-1))
1227 39926046 : index_new=n1*(i2_new-1+n2*(i3_new-1))
1228 39926046 : ind3m=n1*(i2-1+n2*(i3m-1))
1229 39926046 : ind3p=n1*(i2-1+n2*(i3p-1))
1230 39926046 : ind2m=n1*(i2m-1+n2*(i3-1))
1231 39926046 : ind2p=n1*(i2p-1+n2*(i3-1))
1232 39926046 : ind2m3m=n1*(i2m-1+n2*(i3m-1))
1233 1627916196 : do i1=1,n1
1234 1586673758 : ifft=i1+index
1235 1626599804 : if(my_sphere(ifft))then
1236 :
1237 309607368 : i1m=i1-1 ; if(i1==1)i1m=n1
1238 :
1239 309607368 : if(option==1)then
1240 : ! Treat option 1 : computation of residual forces
1241 287723562 : i1p=i1+1 ; if(i1==n1)i1p=1
1242 : ! Distinguish spin-unpolarized and spin-polarized
1243 287723562 : if(dtset%nspden==1)then ! Non magnetic
1244 : ! Note that the factor needed to obtain a true finite difference
1245 : ! estimation of the derivative will be applied afterwards, for speed
1246 255981919 : drho1=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1247 255981919 : drho2=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1248 255981919 : drho3=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1249 255981919 : if(mpi_enreg%paral_kgb==1) then
1250 0 : vresid1=work_tot(ifft,1)
1251 : else
1252 255981919 : vresid1=work(ifft,1)
1253 : end if
1254 255981919 : gresid(1,iatom)=gresid(1,iatom)+drho1*vresid1
1255 255981919 : gresid(2,iatom)=gresid(2,iatom)+drho2*vresid1
1256 255981919 : gresid(3,iatom)=gresid(3,iatom)+drho3*vresid1
1257 31741643 : else if(dtset%nspden==2) then ! Collinear magnetism
1258 28783850 : drho1tot=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1259 28783850 : drho2tot=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1260 28783850 : drho3tot=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1261 28783850 : drho1up=atmrho(i1p+index,2)-atmrho(i1m+index,2)
1262 28783850 : drho2up=atmrho(i1+ind2p,2) -atmrho(i1+ind2m,2)
1263 28783850 : drho3up=atmrho(i1+ind3p,2) -atmrho(i1+ind3m,2)
1264 28783850 : drho1dn=drho1tot-drho1up
1265 28783850 : drho2dn=drho2tot-drho2up
1266 28783850 : drho3dn=drho3tot-drho3up
1267 28783850 : if(mpi_enreg%paral_kgb==1) then
1268 0 : vresid1=work_tot(ifft,1)
1269 0 : vresid2=work_tot(ifft,2)
1270 : else
1271 28783850 : vresid1=work(ifft,1)
1272 28783850 : vresid2=work(ifft,2)
1273 : end if
1274 28783850 : gresid(1,iatom)=gresid(1,iatom)+drho1up*vresid1+drho1dn*vresid2
1275 28783850 : gresid(2,iatom)=gresid(2,iatom)+drho2up*vresid1+drho2dn*vresid2
1276 28783850 : gresid(3,iatom)=gresid(3,iatom)+drho3up*vresid1+drho3dn*vresid2
1277 : else ! Non-collinear magnetism
1278 2957793 : drho1tot=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1279 2957793 : drho1mx =atmrho(i1p+index,2)-atmrho(i1m+index,2)
1280 2957793 : drho1my =atmrho(i1p+index,3)-atmrho(i1m+index,3)
1281 2957793 : drho1mz =atmrho(i1p+index,4)-atmrho(i1m+index,4)
1282 2957793 : drho2tot=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1283 2957793 : drho2mx =atmrho(i1+ind2p,2) -atmrho(i1+ind2m,2)
1284 2957793 : drho2my =atmrho(i1+ind2p,3) -atmrho(i1+ind2m,3)
1285 2957793 : drho2mz =atmrho(i1+ind2p,4) -atmrho(i1+ind2m,4)
1286 2957793 : drho3tot=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1287 2957793 : drho3mx =atmrho(i1+ind3p,2) -atmrho(i1+ind3m,2)
1288 2957793 : drho3my =atmrho(i1+ind3p,3) -atmrho(i1+ind3m,3)
1289 2957793 : drho3mz =atmrho(i1+ind3p,4) -atmrho(i1+ind3m,4)
1290 2957793 : drho11=half*(drho1tot+drho1mz)
1291 2957793 : drho12=half*(drho1tot-drho1mz)
1292 2957793 : drho13= half*drho1mx
1293 2957793 : drho14=-half*drho1my
1294 2957793 : drho21=half*(drho2tot+drho2mz)
1295 2957793 : drho22=half*(drho2tot-drho2mz)
1296 2957793 : drho23= half*drho2mx
1297 2957793 : drho24=-half*drho2my
1298 2957793 : drho31=half*(drho3tot+drho3mz)
1299 2957793 : drho32=half*(drho3tot-drho3mz)
1300 2957793 : drho33= half*drho3mx
1301 2957793 : drho34=-half*drho3my
1302 2957793 : if(mpi_enreg%paral_kgb==1) then
1303 0 : vresid1=work_tot(ifft,1)
1304 0 : vresid2=work_tot(ifft,2)
1305 0 : vresid3=work_tot(ifft,3)
1306 0 : vresid4=work_tot(ifft,4)
1307 : else
1308 2957793 : vresid1=work(ifft,1)
1309 2957793 : vresid2=work(ifft,2)
1310 2957793 : vresid3=work(ifft,3)
1311 2957793 : vresid4=work(ifft,4)
1312 : end if
1313 2957793 : gresid(1,iatom)=gresid(1,iatom)+drho11*vresid1+drho12*vresid2+two*(drho13*vresid3+drho14*vresid4)
1314 2957793 : gresid(2,iatom)=gresid(2,iatom)+drho21*vresid1+drho22*vresid2+two*(drho23*vresid3+drho24*vresid4)
1315 2957793 : gresid(3,iatom)=gresid(3,iatom)+drho31*vresid1+drho32*vresid2+two*(drho33*vresid3+drho34*vresid4)
1316 : end if
1317 : ! Treat the case option==2 now : trilinear interpolation of the density
1318 : else
1319 21883806 : i1_new=i1+diff_igrid(1) ; if(i1_new > n1)i1_new=i1_new-n1
1320 21883806 : ifft_new=i1_new+index_new
1321 46555152 : do ispden=1,dtset%nspden
1322 : drhox00=(atmrho(i1m+index,ispden)-atmrho(i1+index,ispden))*diff_rem1 &
1323 24671346 : & +atmrho(i1+index,ispden)
1324 : drhox10=(atmrho(i1m+ind2m,ispden)-atmrho(i1+ind2m,ispden))*diff_rem1 &
1325 24671346 : & +atmrho(i1+ind2m,ispden)
1326 : drhox01=(atmrho(i1m+ind3m,ispden)-atmrho(i1+ind3m,ispden))*diff_rem1 &
1327 24671346 : & +atmrho(i1+ind3m,ispden)
1328 : drhox11=(atmrho(i1m+ind2m3m,ispden)-atmrho(i1+ind2m3m,ispden))*diff_rem1 &
1329 24671346 : & +atmrho(i1+ind2m3m,ispden)
1330 24671346 : drhoxy0=(drhox10-drhox00)*diff_rem2+drhox00
1331 24671346 : drhoxy1=(drhox11-drhox01)*diff_rem2+drhox01
1332 24671346 : drhoxyz=(drhoxy1-drhoxy0)*diff_rem3+drhoxy0
1333 24671346 : if(mpi_enreg%paral_kgb==1) then
1334 0 : work_tot(ifft_new,ispden)=work_tot(ifft_new,ispden)+drhoxyz
1335 : else
1336 24671346 : work(ifft_new,ispden)=work(ifft_new,ispden)+drhoxyz
1337 : end if
1338 46555152 : rho_tot(ispden)=rho_tot(ispden)+drhoxyz
1339 : end do
1340 : end if
1341 :
1342 : ! End condition of belonging to the sphere of influence of the atom
1343 : end if
1344 : end do
1345 : end do
1346 : end do
1347 : ! The finite-difference factor applied here also take
1348 : ! into account diverse factors
1349 54949 : fact=-ucvol/dble(nfftot)
1350 54949 : gresid(1,iatom)=gresid(1,iatom)*dble(n1)*.5_dp*fact
1351 54949 : gresid(2,iatom)=gresid(2,iatom)*dble(n2)*.5_dp*fact
1352 54949 : gresid(3,iatom)=gresid(3,iatom)*dble(n3)*.5_dp*fact
1353 : end if
1354 :
1355 : ! Update work if the atom is fixed.
1356 56401 : if(iloop==2 .and. option==2 .and. atmove==0)then
1357 1452 : if(mpi_enreg%paral_kgb==1) then
1358 : ! FFT parallelization: All the processors perform the same calculation.
1359 : ! We divided by nproc_fft in order to "remove" the xmpi_sum made in mean_fftr
1360 0 : do ispden=1,dtset%nspden
1361 0 : do ifft=1,nfft_tmp
1362 0 : work_tot(ifft,ispden)=work_tot(ifft,ispden)+atmrho(ifft,ispden)
1363 : end do
1364 : end do
1365 0 : call mean_fftr(atmrho,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1366 0 : rhosum(1:dtset%nspden)=rhosum(1:dtset%nspden)/mpi_enreg%nproc_fft
1367 : else
1368 2940 : do ispden=1,dtset%nspden
1369 10819304 : do ifft=1,nfft_tmp
1370 10817852 : work(ifft,ispden)=work(ifft,ispden)+atmrho(ifft,ispden)
1371 : end do
1372 : end do
1373 1452 : call mean_fftr(atmrho,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1374 : end if
1375 :
1376 2940 : rho_tot(1:dtset%nspden)=rho_tot(1:dtset%nspden)+rhosum(1:dtset%nspden)*nfftot
1377 : end if
1378 :
1379 112802 : ABI_FREE(ii)
1380 268378 : ABI_FREE(rrdiff)
1381 :
1382 : ! End loop on atoms
1383 : end do
1384 :
1385 : ! DEBUG
1386 : ! if(option==2)then
1387 : ! if(iloop==1)then
1388 : ! write(std_out,*)' fresid : rhor, approp'
1389 : ! do ifft=1,n1
1390 : ! write(std_out,*)ifft,rhor(ifft,1),approp(ifft)
1391 : ! end do
1392 : ! end if
1393 : ! if(iloop==2)then
1394 : ! write(std_out,*)' fresid : rhor, approp, work(:,:)'
1395 : ! do ifft=1,n1
1396 : ! write(std_out,'(i4,3es18.8)' )ifft,rhor(ifft,1),approp(ifft),work(ifft,1)
1397 : ! end do
1398 : ! do ifft=1,nfft_tmp
1399 : ! if(work(ifft,1)<0.0_dp)then
1400 : ! write(std_out,*)' f_fft negative value',work(ifft,1),' for ifft=',ifft
1401 : ! end if
1402 : ! if(rhor(ifft,1)<0.0_dp)then
1403 : ! write(std_out,*)' rhor negative value',rhor(ifft,1),' for ifft=',ifft
1404 : ! end if
1405 : ! end do
1406 : ! end if
1407 : ! end if
1408 : ! ENDDEBUG
1409 :
1410 : if(quit==1)exit
1411 :
1412 : ! At the end of the first loop, where the appropriation function is generated,
1413 : ! invert it, to save cpu time later.
1414 382995472 : if(iloop==1)approp(:)=1.0_dp/approp(:)
1415 :
1416 : ! End first or second pass through atoms
1417 : end do
1418 :
1419 : !Restore proper normalisation of density
1420 : !(Non-collinear magnetism: n, mx,my,mz integral conservation)
1421 21387 : if(option==2)then
1422 2017 : if(mpi_enreg%paral_kgb==1) then
1423 : ! FFT parallelization: All the processors perform the same calculation.
1424 : ! We divided by nproc_fft in order to "remove" the xmpi_sum made in mean_fftr
1425 : ! Trangel: mpicomm now is optional in mean_fftr, no need to divide over nproc_fft
1426 0 : call mean_fftr(rhor_tot,rhosum,nfft_tmp,nfftot,dtset%nspden)
1427 : ! rhosum(1:dtset%nspden)=rhosum(1:dtset%nspden)/mpi_enreg%nproc_fft
1428 : else
1429 2017 : call mean_fftr(rhor,rhosum,nfft_tmp,nfftot,dtset%nspden,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1430 : end if
1431 : ! "!OCL NOPREEX" to avoid zero division after optimization (-Of) by MM
1432 : ! (Even if nspden=1, "1.0/rho_tot" will appear on vpp fujitsu
1433 : ! OCL NOPREEX
1434 2017 : if(mpi_enreg%paral_kgb==1) then
1435 0 : do ispden=1,dtset%nspden
1436 0 : fact=rhosum(ispden)*dble(nfftot)/rho_tot(ispden)
1437 0 : work_tot(:,ispden)=fact*work_tot(:,ispden)
1438 0 : call pre_scatter(work(:,ispden),work_tot(:,ispden),n1,n2,n3,n4,mpi_enreg)
1439 : end do
1440 : else
1441 4209 : do ispden=1,dtset%nspden
1442 2192 : fact=rhosum(ispden)*dble(nfftot)/rho_tot(ispden)
1443 38938549 : work(:,ispden)=fact*work(:,ispden)
1444 : end do
1445 : end if
1446 : ! DEBUG
1447 : ! Here, zero all the hard work, for checking purposes !
1448 : ! work(:,:)=rhor(:,:)
1449 : ! ENDDEBUG
1450 : end if
1451 :
1452 21387 : ABI_FREE(approp)
1453 21387 : ABI_FREE(atmrho)
1454 21387 : ABI_FREE(my_sphere)
1455 21387 : if(mpi_enreg%paral_kgb==1) then
1456 0 : ABI_FREE(rhor_tot)
1457 0 : ABI_FREE(work_tot)
1458 : end if
1459 :
1460 : !DEBUG
1461 : !write(std_out,*)' fresid : exit '
1462 : !do iatom=1,dtset%natom
1463 : !write(std_out,*)iatom,gresid(1:3,iatom)
1464 : !enddo
1465 : !ENDDEBUG
1466 :
1467 : contains
1468 :
1469 64161 : function cross_fr(xx,yy,zz,aa,bb,cc)
1470 : !Define magnitude of cross product of two vectors
1471 : real(dp) :: cross_fr
1472 : real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1473 64161 : cross_fr=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1474 64161 : end function cross_fr
1475 :
1476 : end subroutine fresid
1477 : !!***
1478 :
1479 : !!****f* ABINIT/constrf
1480 : !! NAME
1481 : !! constrf
1482 : !!
1483 : !! FUNCTION
1484 : !! Computes projected forces, fpcart, which satisfy a set of
1485 : !! constraint equations of the form
1486 : !! Sum[mu,iatom]: wtatcon(mu,iatom,iconeq)*fpcart(mu,iatom) = 0 (iconeq=1,nconeq).
1487 : !! These projected forces are returned in fcart and thus replace
1488 : !! the original forces.
1489 : !!
1490 : !! INPUTS
1491 : !! iatfix(3,natom)=1 for frozen atom along each direction, 0 for unfrozen
1492 : !! natom=number of atoms in cell
1493 : !! nconeq=number of atomic constraint equations
1494 : !! prtvol=control print volume and debugging
1495 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1496 : !! wtatcon(3,natom,nconeq)=weights for atomic constraints
1497 : !! xred(3,natom)=reduced dimensionless atomic coordinates
1498 : !!
1499 : !! OUTPUT
1500 : !! diffor=maximum absolute change in component of projected fcart between present
1501 : !! and previous SCF cycle
1502 : !! gred(3,natom)=grads of Etot wrt reduced coordinates (hartree)
1503 : !! maxfor=maximum absolute value of fcart
1504 : !!
1505 : !! SIDE EFFECTS
1506 : !! fcart(3,natom)=cartesian forces (hartree/bohr) on input, projected forces on output
1507 : !! forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
1508 : !!
1509 : !! TODO
1510 : !!
1511 : !! SOURCE
1512 :
1513 740 : subroutine constrf(diffor,fcart,forold,gred,iatfix,ionmov,maxfor,natom,&
1514 740 : & nconeq,prtvol,rprimd,wtatcon,xred)
1515 :
1516 : use m_linalg_interfaces
1517 :
1518 : !Arguments ------------------------------------
1519 : !scalars
1520 : integer,intent(in) :: ionmov,natom,nconeq,prtvol
1521 : real(dp),intent(out) :: diffor,maxfor
1522 : !arrays
1523 : integer,intent(in) :: iatfix(3,natom)
1524 : real(dp),intent(in) :: rprimd(3,3),wtatcon(3,natom,nconeq)
1525 : real(dp),intent(inout) :: fcart(3,natom),forold(3,natom),xred(3,natom)
1526 : real(dp),intent(inout) :: gred(3,natom) !vz_i
1527 :
1528 : !Local variables -------------------------
1529 : !scalars
1530 : integer :: iatom,iconeq,iconeq1,iconeq2,index,info,mu,prtvel
1531 : character(len=500) :: message
1532 : !arrays
1533 740 : real(dp),allocatable :: fcartvec(:),fpcart(:,:),fvector(:),vel_dummy(:,:)
1534 740 : real(dp),allocatable :: wmatrix(:,:),wtatconvec(:,:),wtcoeffs(:),xcart(:,:)
1535 :
1536 : !************************************************************************
1537 :
1538 : !Allocate temporary variables
1539 2220 : ABI_MALLOC(fpcart,(3,natom))
1540 2960 : ABI_MALLOC(fcartvec,(3*natom))
1541 2220 : ABI_MALLOC(fvector,(nconeq))
1542 1480 : ABI_MALLOC(vel_dummy,(3,natom))
1543 2960 : ABI_MALLOC(wmatrix,(nconeq,nconeq))
1544 2960 : ABI_MALLOC(wtatconvec,(3*natom,nconeq))
1545 1480 : ABI_MALLOC(wtcoeffs,(nconeq))
1546 1480 : ABI_MALLOC(xcart,(3,natom))
1547 :
1548 : !If prtvol>10, output coordinates and forces prior to projecting
1549 740 : if(prtvol>=10)then
1550 132 : write(message,'(a)')' constrf - coordinates and forces prior to constraint projections:'
1551 132 : call wrtout(std_out,message,'COLL')
1552 132 : call xred2xcart(natom,rprimd,xcart,xred)
1553 132 : prtvel=0
1554 132 : call prtxvf(fcart,gred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred)
1555 : end if
1556 :
1557 : !Transfer fcart and wtatcon to flat vectors
1558 : index=0
1559 4040 : do iatom=1,natom
1560 13940 : do mu=1,3
1561 9900 : index=index+1
1562 9900 : fcartvec(index)=fcart(mu,iatom)
1563 23100 : wtatconvec(index,:)=wtatcon(mu,iatom,:)
1564 : end do
1565 : end do
1566 :
1567 : !Compute a matrix (wmatrix) and vector (fvector) such that solving
1568 : !the linear equations wmatrix*wcoeffs=fvector gives the coefficients
1569 : !of wtatcon (wcoeffs) needed to compute the projected forces
1570 1480 : do iconeq2=1,nconeq
1571 740 : fvector(iconeq2)=ddot(3*natom,fcartvec,1,wtatconvec(1,iconeq2),1)
1572 2220 : do iconeq1=1,nconeq
1573 1480 : wmatrix(iconeq1,iconeq2)=ddot(3*natom,wtatconvec(1,iconeq1),1,wtatconvec(1,iconeq2),1)
1574 : end do
1575 : end do
1576 :
1577 : !Solve the system of linear equations, wmatrix*wcoeffs=fvector
1578 740 : call dposv('U',nconeq,1,wmatrix,nconeq,fvector,nconeq,info)
1579 :
1580 740 : if (info/=0) then
1581 : write(message, '(a,a,a,a,a)' )&
1582 0 : & 'Constraint matrix is not positive definite,',ch10,&
1583 0 : & 'probably because constraints are linearly dependent.',ch10,&
1584 0 : & 'Action: Check for linear dependence of constraints.'
1585 0 : ABI_ERROR(message)
1586 : end if
1587 :
1588 : !The solution vector is returned in fvector, so copy it to a more sensible location
1589 1480 : wtcoeffs(:)=fvector(:)
1590 :
1591 : !Compute the projected forces, which now satisfy all the constraint equations
1592 13940 : fpcart(:,:)=fcart(:,:)
1593 1480 : do iconeq=1,nconeq
1594 14680 : fpcart(:,:)=fpcart(:,:)-wtcoeffs(iconeq)*wtatcon(:,:,iconeq)
1595 : end do
1596 :
1597 : !Reconvert constrained forces back from fpcart to gred
1598 4040 : do iatom=1,natom
1599 13940 : do mu=1,3
1600 : gred(mu,iatom)= - (rprimd(1,mu)*fpcart(1,iatom)+&
1601 : & rprimd(2,mu)*fpcart(2,iatom)+&
1602 13200 : & rprimd(3,mu)*fpcart(3,iatom))
1603 : end do
1604 : end do
1605 :
1606 : !If prtvol>=10, output coordinates and forces after projecting
1607 740 : if(prtvol>=10)then
1608 132 : write(message,'(a)')' constrf - coordinates and forces after constraint projections:'
1609 132 : call wrtout(std_out,message,'COLL')
1610 132 : prtvel=0
1611 132 : call prtxvf(fpcart,gred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred)
1612 : end if
1613 :
1614 : !Copy the constrained forces, fpcart, back to fcart
1615 13940 : fcart(:,:)=fpcart(:,:)
1616 :
1617 : !Compute maximal force and maximal difference of the projected forces,
1618 : !overriding the values already computed in forces
1619 740 : maxfor=0.0_dp
1620 740 : diffor=0.0_dp
1621 4040 : do iatom=1,natom
1622 13940 : do mu=1,3
1623 13200 : if (iatfix(mu,iatom) /= 1) then
1624 7263 : maxfor=max(maxfor,abs(fcart(mu,iatom)))
1625 7263 : diffor=max(diffor,abs(fcart(mu,iatom)-forold(mu,iatom)))
1626 2637 : else if (ionmov==4 .or. ionmov==5) then
1627 : ! Make the force vanish on fixed atoms when ionmov=4 or 5
1628 : ! This is because fixing of atom cannot be imposed at the
1629 : ! level of a routine similar to brdmin or moldyn for these options.
1630 0 : fcart(mu,iatom)=0.0_dp
1631 : end if
1632 : end do
1633 : end do
1634 13940 : forold(:,:)=fcart(:,:)
1635 :
1636 : !Dellocate temporary variables
1637 740 : ABI_FREE(fpcart)
1638 740 : ABI_FREE(fcartvec)
1639 740 : ABI_FREE(fvector)
1640 740 : ABI_FREE(vel_dummy)
1641 740 : ABI_FREE(wmatrix)
1642 740 : ABI_FREE(wtatconvec)
1643 740 : ABI_FREE(wtcoeffs)
1644 740 : ABI_FREE(xcart)
1645 :
1646 740 : end subroutine constrf
1647 : !!***
1648 :
1649 : end module m_forces
1650 : !!***
|
