Branch data 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 : 24016 : subroutine forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,&
173 : 24016 : & forold,gred,grchempottn,grcondft,gresid,grewtn,&
174 : 24016 : & grhf,grnl,grvdw,grxc,gsqcut,indsym,&
175 : : & maxfor,mgfft,mpi_enreg,n1xccc,n3xccc,&
176 : 24016 : & nattyp,nfft,ngfft,ngrvdw,ntypat,&
177 : 24016 : & pawrad,pawtab,ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,&
178 : : & usefock,usevxctau,&
179 : 24016 : & 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 : 24016 : real(dp),allocatable :: atmrho_dum(:),atmvloc_dum(:),dyfrlo_dum(:,:,:)
224 : 24016 : real(dp),allocatable :: dyfrn_dum(:,:,:),dyfrv_dum(:,:,:)
225 : 24016 : real(dp),allocatable :: dyfrx2_dum(:,:,:),eltfrn_dum(:,:),gauss_dum(:,:)
226 : 24016 : real(dp),allocatable :: epawf3red(:,:),fin(:,:),fionred(:,:),grl(:,:),grl_dum(:,:)
227 : 24016 : real(dp),allocatable :: grnl_tmp(:,:),grtn(:,:),grtn_indx(:,:),grxctau(:,:),v_dum(:),vxctotg(:,:)
228 : 24016 : real(dp),allocatable :: xccc3d_dum(:)
229 : :
230 : : ! *************************************************************************
231 : :
232 : 24016 : call timab(69,1,tsec)
233 : : ABI_NVTX_START_RANGE(NVTX_FORCES)
234 : :
235 : : !Save input value of forces
236 [ - + + - : 72048 : ABI_MALLOC(fin,(3,dtset%natom))
+ - + - -
+ - + ]
237 [ + + + + ]: 288540 : fin(:,:)=fcart(:,:)
238 : :
239 : : !Compute different geometric tensor, as well as ucvol, from rprimd
240 : 24016 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
241 : :
242 : : !Test size of vxctau
243 [ + + ]: 96064 : has_vxctau=size(vxctau)
244 [ + + - + ]: 24016 : 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 : 24016 : 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 [ + + + + : 24016 : & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
+ + ]
253 : :
254 : : !=======================================================================
255 : : !========= Local pseudopotential and core charge contributions =========
256 : : !=======================================================================
257 : :
258 [ - + + - : 48032 : 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 [ + + ]: 24016 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
266 [ - + ]: 24016 : if (dtset%icoulomb>0) vloc_method=2
267 [ - + ]: 24016 : 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 [ + + ]: 24016 : if (psps%nc_xccc_gspace==1) coredens_method=1
273 [ + + ]: 24016 : if (psps%nc_xccc_gspace==0) coredens_method=2
274 [ - + ]: 24016 : if (psps%usewvl==1) coredens_method=2
275 : : !Pseudo core charge kinetic energy density: same convention as coredens_method
276 : 24016 : coretau_method=0
277 [ + + ]: 24016 : 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 [ + + - + ]: 24016 : if (vloc_method==1.or.coredens_method==1.or.coretau_method==1) then
283 [ + - + + : 6944 : 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 : 6944 : call timab(550,1,tsec)
287 : : ! Allocate (unused) dummy variables, otherwise some compilers complain
288 [ - + ]: 6944 : ABI_MALLOC(gauss_dum,(0,0))
289 [ - + ]: 6944 : ABI_MALLOC(atmrho_dum,(0))
290 [ - + ]: 6944 : ABI_MALLOC(atmvloc_dum,(0))
291 [ - + ]: 6944 : ABI_MALLOC(dyfrn_dum,(0,0,0))
292 [ - + ]: 6944 : ABI_MALLOC(dyfrv_dum,(0,0,0))
293 [ - + ]: 6944 : ABI_MALLOC(eltfrn_dum,(0,0))
294 : : ! Compute Vxc in reciprocal space
295 [ + - + + ]: 6944 : if (coredens_method==1.and.n3xccc>0) then
296 [ - + + - : 13290 : ABI_MALLOC(v_dum,(nfft))
+ - + - -
+ - + ]
297 [ - + + - : 13290 : ABI_MALLOC(vxctotg,(2,nfft))
+ - + - -
+ - + ]
298 [ + + + + : 161215534 : v_dum(:)=vxc(:,1);if (dtset%nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
+ + ]
299 : 4430 : 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 : 4430 : & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
302 [ - + ]: 4430 : 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 [ + + ]: 6944 : optv=0;if (vloc_method==1) optv=1
308 [ + - ]: 6944 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
309 : 6944 : optatm=0;optdyfr=0;optgr=1;optstr=0;optn2=1;opteltfr=0
310 [ + - ]: 6944 : 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 : 6944 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
318 : : end if
319 [ + + + - : 43748 : if (n3xccc==0.and.coredens_method==1) grxc=zero
+ + + + ]
320 [ - + ]: 6944 : ABI_FREE(vxctotg)
321 [ + + + - : 6944 : 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 [ - + ]: 6944 : ABI_FREE(gauss_dum)
347 [ - + ]: 6944 : ABI_FREE(atmrho_dum)
348 [ - + ]: 6944 : ABI_FREE(atmvloc_dum)
349 [ - + ]: 6944 : ABI_FREE(dyfrn_dum)
350 [ - + ]: 6944 : ABI_FREE(dyfrv_dum)
351 [ - + ]: 6944 : ABI_FREE(eltfrn_dum)
352 : 6944 : call timab(550,2,tsec)
353 : : end if
354 : :
355 : : !Local ionic potential by method 2
356 [ + + ]: 24016 : if (vloc_method==2) then
357 : 17075 : option=2
358 [ - + + - : 51225 : ABI_MALLOC(dyfrlo_dum,(3,3,dtset%natom))
+ - + - -
+ - + ]
359 [ - + + - : 34150 : ABI_MALLOC(grtn_indx,(3,dtset%natom))
- + - + ]
360 [ - + + - : 51225 : 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 : 17075 : & qprtrb_dum,rhog,rhor,rprimd,ucvol,vprtrb_dum,v_dum,wvl,wvl_den,xred)
364 [ + + ]: 61265 : do iatom=1,dtset%natom
365 : : ! Has to use the indexing array atindx1
366 [ + + ]: 193835 : grl(1:3,atindx1(iatom))=grtn_indx(1:3,iatom)
367 : : end do
368 [ - + ]: 17075 : ABI_FREE(dyfrlo_dum)
369 [ - + ]: 17075 : ABI_FREE(grtn_indx)
370 [ - + ]: 17075 : 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 [ - + ]: 17075 : 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 [ + + ]: 24016 : if (coredens_method==2.or.coretau_method==2) then
383 [ + + ]: 17072 : if (n1xccc/=0) then
384 : 4101 : call timab(53,1,tsec)
385 : 4101 : option=2
386 [ - + + - : 12303 : ABI_MALLOC(dyfrx2_dum,(3,3,dtset%natom))
+ - + - -
+ - + ]
387 [ - + + - : 12303 : ABI_MALLOC(xccc3d_dum,(n3xccc))
+ - + - -
+ - + ]
388 [ + - ]: 4101 : if (coredens_method==2) then
389 [ + + ]: 4101 : 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 [ + - + - ]: 4001 : 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 : 4001 : & 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 [ + + + - : 4101 : 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 [ - + ]: 4101 : ABI_FREE(xccc3d_dum)
427 [ - + ]: 4101 : ABI_FREE(dyfrx2_dum)
428 : 4101 : 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 [ - + + - : 48032 : ABI_MALLOC(grnl_tmp,(3,dtset%natom))
- + - + ]
440 [ + + ]: 90147 : do iatom=1,dtset%natom
441 [ + + ]: 288540 : indx=3*(iatom-1);grnl_tmp(1:3,atindx1(iatom))=grnl(indx+1:indx+3)
442 : : end do
443 [ + - ]: 24016 : if (dtset%usewvl == 0) then
444 : 24016 : call sygrad(synlgr,dtset%natom,grnl_tmp,dtset%nsym,symrec,indsym)
445 : : else
446 [ # # # # ]: 0 : synlgr = grnl_tmp
447 : : end if
448 [ - + ]: 24016 : ABI_FREE(grnl_tmp)
449 : :
450 : : !=======================================================================
451 : : !============ Density/potential residual contributions =================
452 : : !=======================================================================
453 : :
454 [ + - + + ]: 24016 : 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 : 19625 : & 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 [ + + ]: 24016 : & 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 [ - + + - : 48032 : ABI_MALLOC(grtn,(3,dtset%natom))
- + - + ]
539 [ + + + + ]: 288540 : grtn(:,:)=grl(:,:)+grchempottn(:,:)+grcondft(:,:)+grewtn(:,:)+synlgr(:,:)+grxc(:,:)
540 : :
541 [ + + + - ]: 24016 : 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 [ + + + + : 24088 : 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 [ + + + + : 24440 : 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 : 24016 : 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 [ + - + + ]: 24016 : if (dtset%usewvl==0.and.abs(dtset%densfor_pred)>0.and.abs(dtset%densfor_pred)/=5) then
560 [ + + ]: 90016 : do iatom=1,dtset%natom
561 [ + + ]: 288070 : do mu=1,3
562 [ + + ]: 264072 : 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 [ + + + + ]: 288540 : grtn(:,:)=grtn(:,:)+gresid(:,:)
569 : :
570 : : !Additional stuff for electron-positron
571 : 24016 : ipositron=0
572 [ + - ]: 24016 : if (present(electronpositron)) then
573 [ + + ]: 24016 : 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 [ + + + + : 39264 : if (ipositron>0) grtn(:,:)=grtn(:,:)+electronpositron%gred_ep(:,:)
+ + ]
589 : :
590 : : !Symmetrize all grads explicitly for given space group:
591 [ + - ]: 24016 : if (dtset%usewvl == 0) then
592 : 24016 : 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 [ + + + + ]: 24208 : call gred2fcart(favg,(dtset%jellslab==0 .and. dtset%nzchempot==0),fcart,gred,gprimd,dtset%natom)
604 : :
605 : : !Compute maximal force and maximal difference
606 : 24016 : maxfor=zero;diffor=zero
607 [ + + ]: 90147 : do iatom=1,dtset%natom
608 [ + + ]: 288540 : do mu=1,3
609 [ + + ]: 264524 : if (dtset%iatfix(mu,iatom) /= 1) then
610 : 167365 : maxfor=max(maxfor,abs(fcart(mu,iatom)))
611 : 167365 : 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 [ + + ]: 24016 : 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 [ - + ]: 24016 : ABI_FREE(grl)
628 [ - + ]: 24016 : ABI_FREE(grtn)
629 [ - + ]: 24016 : ABI_FREE(fin)
630 [ + + ]: 24016 : 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 : 24016 : call timab(69,2,tsec)
639 : :
640 [ - + - + ]: 24016 : 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 : 72048 : 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 [ + + ]: 72048 : if (nsym==1) then
703 : : ! only symmetry is identity so simply copy
704 [ + + + + ]: 338490 : gred(:,:)=dedt(:,:)
705 : : else
706 : : ! actually conduct symmetrization
707 [ + + ]: 169254 : do ia=1,natom
708 [ + + ]: 527130 : do mu=1,3
709 : : summ=0._dp
710 [ + + ]: 4347846 : do isym=1,nsym
711 : 3989970 : 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 : 4347846 : & dble(symrec(mu,3,isym))*dedt(3,ind)
715 : : end do
716 : 357876 : gred(mu,ia)=summ/dble(nsym)
717 [ + + ]: 477168 : if(abs(gred(mu,ia))<tol)gred(mu,ia)=0.0_dp
718 : : end do
719 : : end do
720 : : end if
721 : :
722 : 72048 : 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 : 21642 : subroutine fresid(dtset,gresid,mpi_enreg,nfft,ngfft,ntypat,option,&
899 : 21642 : & 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 : 21642 : 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 : 21642 : real(dp),allocatable :: approp(:),atmrho(:,:),rhor_tot(:,:),rrdiff(:,:)
940 : 21642 : real(dp),allocatable :: work_tot(:,:)
941 : 21642 : 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 : 21642 : & rprimd(1,3),rprimd(2,3),rprimd(3,3))
949 : : lencp(2)=cross_fr(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
950 : 21642 : & rprimd(1,1),rprimd(2,1),rprimd(3,1))
951 : : lencp(3)=cross_fr(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
952 : 21642 : & 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 [ + + ]: 86568 : scale(:)=ucvol/lencp(:)
957 : :
958 : : !initialize diff_igrid, otherwise valgrind complains
959 : 21642 : diff_igrid=0
960 : :
961 : : !Compute metric tensor in real space rmet
962 [ + + ]: 86568 : do nu=1,3
963 [ + + ]: 281346 : 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 : 21642 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
969 : 21642 : n4=n3/mpi_enreg%nproc_fft
970 [ + + ]: 86568 : nfftot=PRODUCT(ngfft(1:3));nfft_tmp=nfftot
971 [ - + ]: 21642 : 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 [ + + + + ]: 249270 : gresid(1:3,1:dtset%natom)=0.0_dp
981 : 21642 : quit=0
982 : :
983 : : !Initialize appropriation function
984 [ - + + - : 64926 : ABI_MALLOC(approp,(nfft_tmp))
+ - + - -
+ - + ]
985 [ - + + - : 86568 : ABI_MALLOC(atmrho,(nfft_tmp,dtset%nspden))
- + + - +
- + - - +
- + ]
986 [ - + + - : 64926 : ABI_MALLOC(my_sphere,(nfft_tmp))
+ - + - -
+ - + ]
987 : :
988 [ + + ]: 386606898 : approp(:)=app_remain
989 : : !First loop over atoms in unit cell : build appropriation function
990 : : !Second loop : compute forces
991 [ + + ]: 64926 : do iloop=1,2
992 : :
993 : : ! Take into account the remaining density
994 [ + + + + ]: 43284 : 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 [ + + ]: 157098 : do iatom=1,dtset%natom
1019 : :
1020 : : ! Get the covalent radius
1021 : 113814 : call atomdata_from_znucl(atom,znucl(dtset%typat(iatom)))
1022 : 113814 : rcov = atom%rcov
1023 : : ! PAW choose PAW radius instead...
1024 [ + + ]: 113814 : if (dtset%usepaw==1) rcov=max(rcov,pawtab(dtset%typat(iatom))%rpaw)
1025 : :
1026 : : ! Set search range
1027 : 113814 : rcov2=rcov**2
1028 : 113814 : range=2._dp*rcov
1029 : 113814 : range2=range**2
1030 : 113814 : 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 : 113814 : irange(1)=1+nint((range/scale(1))*dble(n1))
1035 : 113814 : irange(2)=1+nint((range/scale(2))*dble(n2))
1036 : 113814 : irange(3)=1+nint((range/scale(3))*dble(n3))
1037 : :
1038 : : ! Allocate ii and rrdiff
1039 [ + + ]: 455256 : mshift=2*maxval(irange(1:3))+1
1040 [ - + + - : 341442 : ABI_MALLOC(ii,(mshift,3))
+ - + - -
+ - + ]
1041 [ - + + - : 341442 : ABI_MALLOC(rrdiff,(mshift,3))
+ - + - -
+ - + ]
1042 : :
1043 : : ! Consider each component in turn
1044 [ + + ]: 455256 : do mu=1,3
1045 : :
1046 : : ! Convert reduced coord of given atom to [0,1)
1047 : 341442 : 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 : 341442 : 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 : 341442 : ishift=0
1056 [ + + ]: 7237788 : do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
1057 : 6896346 : ishift=ishift+1
1058 : 6896346 : ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
1059 : 7237788 : 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 [ + + + + ]: 455256 : 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 : 113814 : atmove=1
1085 [ + + + + ]: 113814 : 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 [ + + ]: 108129 : if(iloop==2) then
1098 [ + + + + ]: 1717183922 : atmrho(:,:)=0.0_dp
1099 [ + + ]: 1603803125 : my_sphere(:)=.false.
1100 : : end if
1101 : :
1102 : : ! Conduct triple loop over restricted range of grid points for iatom
1103 : :
1104 [ + + ]: 2397664 : do ishift3=1,1+2*irange(3)
1105 : : ! map back to [1,ngfft(3)] for usual fortran index in unit cell
1106 : 2283850 : i3=ii(ishift3,3)
1107 [ + + ]: 2283850 : i3m=i3-1 ; if(i3==1)i3m=n3
1108 [ + + ]: 2283850 : i3p=i3+1 ; if(i3==n3)i3p=1
1109 : :
1110 : : ! find vector from atom location to grid point (reduced)
1111 : 2283850 : rdiff3=rrdiff(ishift3,3)
1112 : :
1113 [ + + ]: 63108146 : do ishift2=1,1+2*irange(2)
1114 : 60710482 : i2=ii(ishift2,2)
1115 [ + + ]: 60710482 : i2m=i2-1 ; if(i2==1)i2m=n2
1116 [ + + ]: 60710482 : i2p=i2+1 ; if(i2==n2)i2p=1
1117 : 60710482 : index=n1*(i2-1+n2*(i3-1))
1118 : 60710482 : ind3m=n1*(i2-1+n2*(i3m-1))
1119 : 60710482 : ind3p=n1*(i2-1+n2*(i3p-1))
1120 : 60710482 : ind2m=n1*(i2m-1+n2*(i3-1))
1121 : 60710482 : ind2p=n1*(i2p-1+n2*(i3-1))
1122 : 60710482 : ind2p3p=n1*(i2p-1+n2*(i3p-1))
1123 : :
1124 : 60710482 : rdiff2=rrdiff(ishift2,2)
1125 : : ! Prepare the computation of difmag2
1126 : : difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
1127 : 60710482 : & +2.0_dp*rmet(3,2)*rdiff3*rdiff2
1128 : 60710482 : difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
1129 : :
1130 [ + + ]: 2137564762 : do ishift1=1,1+2*irange(1)
1131 : 2074570430 : rdiff1=rrdiff(ishift1,1)
1132 : :
1133 : : ! Compute (rgrid-tau-Rprim)**2
1134 : 2074570430 : 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 [ + + ]: 2135280912 : if (difmag2<range2) then
1139 : :
1140 : : ! Will compute contribution to appropriation function based on
1141 : : ! rcov2, range2 and difmag2
1142 : 697428866 : i1=ii(ishift1,1)
1143 : 697428866 : ifft=i1+index
1144 : :
1145 [ + + ]: 697428866 : if(iloop==1)then
1146 : :
1147 : : ! Build appropriation function
1148 [ + + ]: 348714433 : if (difmag2<rcov2)then
1149 : 43601822 : approp(ifft)=approp(ifft)+1.0_dp
1150 : : else
1151 : 305112611 : difmag=sqrt(difmag2)
1152 : 305112611 : xx=difmag*rcovm1
1153 : : ! The following function is 1. at xx=1, 0. at xx=2, with vanishing
1154 : : ! derivatives at these points.
1155 : 305112611 : 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 [ + + ]: 348714433 : if (difmag2<rcov2) then
1161 : : fact=one
1162 : : else
1163 : 305112611 : difmag=sqrt(difmag2)
1164 : 305112611 : xx=difmag*rcovm1
1165 : 305112611 : 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 [ - + ]: 348714433 : 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 [ + + ]: 753266732 : & +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 [ + + ]: 348714433 : i1p=i1+1 ; if(i1==n1)i1p=1
1180 [ + + ]: 348714433 : if(option==1)then
1181 [ + + ]: 317335601 : i1m=i1-1 ; if(i1==1)i1m=n1
1182 : 317335601 : my_sphere(ifft)=.true.
1183 : 317335601 : my_sphere(i1p+index)=.true. ; my_sphere(i1m+index)=.true.
1184 : 317335601 : my_sphere(i1+ind2p)=.true. ; my_sphere(i1+ind2m)=.true.
1185 : 317335601 : 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 [ + + + + : 113814 : if(iloop==2 .and. (option==1 .or. atmove==1) )then
+ + ]
1216 : :
1217 [ + + ]: 1384091 : do i3=1,n3
1218 [ + + ]: 1328636 : i3m=i3-1 ; if(i3==1)i3m=n3
1219 [ + + ]: 1328636 : i3p=i3+1 ; if(i3==n3)i3p=1
1220 : : ! note: diff_igrid is only set if(option==2 .and. iloop==2)
1221 [ + + ]: 1328636 : i3_new=i3+diff_igrid(3) ; if(i3_new > n3)i3_new=i3_new-n3
1222 [ + + ]: 41609017 : do i2=1,n2
1223 [ + + ]: 40224926 : i2m=i2-1 ; if(i2==1)i2m=n2
1224 [ + + ]: 40224926 : i2p=i2+1 ; if(i2==n2)i2p=1
1225 [ + + ]: 40224926 : i2_new=i2+diff_igrid(2) ; if(i2_new > n2)i2_new=i2_new-n2
1226 : 40224926 : index=n1*(i2-1+n2*(i3-1))
1227 : 40224926 : index_new=n1*(i2_new-1+n2*(i3_new-1))
1228 : 40224926 : ind3m=n1*(i2-1+n2*(i3m-1))
1229 : 40224926 : ind3p=n1*(i2-1+n2*(i3p-1))
1230 : 40224926 : ind2m=n1*(i2m-1+n2*(i3-1))
1231 : 40224926 : ind2p=n1*(i2p-1+n2*(i3-1))
1232 : 40224926 : ind2m3m=n1*(i2m-1+n2*(i3m-1))
1233 [ + + ]: 1635638296 : do i1=1,n1
1234 : 1594084734 : ifft=i1+index
1235 [ + + ]: 1634309660 : if(my_sphere(ifft))then
1236 : :
1237 [ + + ]: 316804670 : i1m=i1-1 ; if(i1==1)i1m=n1
1238 : :
1239 [ + + ]: 316804670 : if(option==1)then
1240 : : ! Treat option 1 : computation of residual forces
1241 [ + + ]: 294920864 : i1p=i1+1 ; if(i1==n1)i1p=1
1242 : : ! Distinguish spin-unpolarized and spin-polarized
1243 [ + + ]: 294920864 : 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 : 262754595 : drho1=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1247 : 262754595 : drho2=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1248 : 262754595 : drho3=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1249 [ - + ]: 262754595 : if(mpi_enreg%paral_kgb==1) then
1250 : 0 : vresid1=work_tot(ifft,1)
1251 : : else
1252 : 262754595 : vresid1=work(ifft,1)
1253 : : end if
1254 : 262754595 : gresid(1,iatom)=gresid(1,iatom)+drho1*vresid1
1255 : 262754595 : gresid(2,iatom)=gresid(2,iatom)+drho2*vresid1
1256 : 262754595 : gresid(3,iatom)=gresid(3,iatom)+drho3*vresid1
1257 [ + + ]: 32166269 : else if(dtset%nspden==2) then ! Collinear magnetism
1258 : 29090378 : drho1tot=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1259 : 29090378 : drho2tot=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1260 : 29090378 : drho3tot=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1261 : 29090378 : drho1up=atmrho(i1p+index,2)-atmrho(i1m+index,2)
1262 : 29090378 : drho2up=atmrho(i1+ind2p,2) -atmrho(i1+ind2m,2)
1263 : 29090378 : drho3up=atmrho(i1+ind3p,2) -atmrho(i1+ind3m,2)
1264 : 29090378 : drho1dn=drho1tot-drho1up
1265 : 29090378 : drho2dn=drho2tot-drho2up
1266 : 29090378 : drho3dn=drho3tot-drho3up
1267 [ - + ]: 29090378 : if(mpi_enreg%paral_kgb==1) then
1268 : 0 : vresid1=work_tot(ifft,1)
1269 : 0 : vresid2=work_tot(ifft,2)
1270 : : else
1271 : 29090378 : vresid1=work(ifft,1)
1272 : 29090378 : vresid2=work(ifft,2)
1273 : : end if
1274 : 29090378 : gresid(1,iatom)=gresid(1,iatom)+drho1up*vresid1+drho1dn*vresid2
1275 : 29090378 : gresid(2,iatom)=gresid(2,iatom)+drho2up*vresid1+drho2dn*vresid2
1276 : 29090378 : gresid(3,iatom)=gresid(3,iatom)+drho3up*vresid1+drho3dn*vresid2
1277 : : else ! Non-collinear magnetism
1278 : 3075891 : drho1tot=atmrho(i1p+index,1)-atmrho(i1m+index,1)
1279 : 3075891 : drho1mx =atmrho(i1p+index,2)-atmrho(i1m+index,2)
1280 : 3075891 : drho1my =atmrho(i1p+index,3)-atmrho(i1m+index,3)
1281 : 3075891 : drho1mz =atmrho(i1p+index,4)-atmrho(i1m+index,4)
1282 : 3075891 : drho2tot=atmrho(i1+ind2p,1) -atmrho(i1+ind2m,1)
1283 : 3075891 : drho2mx =atmrho(i1+ind2p,2) -atmrho(i1+ind2m,2)
1284 : 3075891 : drho2my =atmrho(i1+ind2p,3) -atmrho(i1+ind2m,3)
1285 : 3075891 : drho2mz =atmrho(i1+ind2p,4) -atmrho(i1+ind2m,4)
1286 : 3075891 : drho3tot=atmrho(i1+ind3p,1) -atmrho(i1+ind3m,1)
1287 : 3075891 : drho3mx =atmrho(i1+ind3p,2) -atmrho(i1+ind3m,2)
1288 : 3075891 : drho3my =atmrho(i1+ind3p,3) -atmrho(i1+ind3m,3)
1289 : 3075891 : drho3mz =atmrho(i1+ind3p,4) -atmrho(i1+ind3m,4)
1290 : 3075891 : drho11=half*(drho1tot+drho1mz)
1291 : 3075891 : drho12=half*(drho1tot-drho1mz)
1292 : 3075891 : drho13= half*drho1mx
1293 : 3075891 : drho14=-half*drho1my
1294 : 3075891 : drho21=half*(drho2tot+drho2mz)
1295 : 3075891 : drho22=half*(drho2tot-drho2mz)
1296 : 3075891 : drho23= half*drho2mx
1297 : 3075891 : drho24=-half*drho2my
1298 : 3075891 : drho31=half*(drho3tot+drho3mz)
1299 : 3075891 : drho32=half*(drho3tot-drho3mz)
1300 : 3075891 : drho33= half*drho3mx
1301 : 3075891 : drho34=-half*drho3my
1302 [ - + ]: 3075891 : 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 : 3075891 : vresid1=work(ifft,1)
1309 : 3075891 : vresid2=work(ifft,2)
1310 : 3075891 : vresid3=work(ifft,3)
1311 : 3075891 : vresid4=work(ifft,4)
1312 : : end if
1313 : 3075891 : gresid(1,iatom)=gresid(1,iatom)+drho11*vresid1+drho12*vresid2+two*(drho13*vresid3+drho14*vresid4)
1314 : 3075891 : gresid(2,iatom)=gresid(2,iatom)+drho21*vresid1+drho22*vresid2+two*(drho23*vresid3+drho24*vresid4)
1315 : 3075891 : 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 : 55455 : fact=-ucvol/dble(nfftot)
1350 : 55455 : gresid(1,iatom)=gresid(1,iatom)*dble(n1)*.5_dp*fact
1351 : 55455 : gresid(2,iatom)=gresid(2,iatom)*dble(n2)*.5_dp*fact
1352 : 55455 : gresid(3,iatom)=gresid(3,iatom)*dble(n3)*.5_dp*fact
1353 : : end if
1354 : :
1355 : : ! Update work if the atom is fixed.
1356 [ + + + + ]: 56907 : 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 [ - + ]: 113814 : ABI_FREE(ii)
1380 [ - + ]: 270912 : 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 [ + + + + ]: 386650182 : 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 [ + + ]: 21642 : 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 [ - + ]: 21642 : ABI_FREE(approp)
1453 [ - + ]: 21642 : ABI_FREE(atmrho)
1454 [ - + ]: 21642 : ABI_FREE(my_sphere)
1455 [ - + - + : 21642 : 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 : 64926 : 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 : 64926 : cross_fr=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1474 : 64926 : 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 : : !!***
|
