Branch data Line data Source code
1 : : !!****m* ABINIT/m_mklocl_realspace
2 : : !! NAME
3 : : !! m_mklocl_realspace
4 : : !!
5 : : !! FUNCTION
6 : : !! Routines related to the local part of the pseudopotentials.
7 : : !! Computation is done in real space (useful for isolated systems).
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2013-2026 ABINIT group (TRangel, MT, DC)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! TODO
16 : : !! This module could be merged with m_mklocl
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : module m_mklocl_realspace
27 : :
28 : : use defs_basis
29 : : use defs_wvltypes
30 : : use m_xmpi
31 : : use m_abicore
32 : : use m_errors
33 : :
34 : : use defs_datatypes, only : pseudopotential_type
35 : : use defs_abitypes, only : MPI_type
36 : : use m_time, only : timab
37 : : use m_geometry, only : xred2xcart
38 : : use m_fft_mesh, only : mkgrid_fft
39 : : use m_mpinfo, only : ptabs_fourdp
40 : : use m_pawtab, only : pawtab_type
41 : : use m_paw_numeric, only : paw_splint, paw_splint_der
42 : : use m_psolver, only : psolver_hartree, psolver_kernel
43 : : use m_abi2big, only : wvl_rhov_abi2big
44 : : use m_wvl_wfs, only : derf_ab
45 : : use m_fft, only : fourdp
46 : :
47 : : implicit none
48 : :
49 : : private
50 : : !!***
51 : :
52 : : public :: mklocl_realspace
53 : : public :: mklocl_wavelets
54 : : !!***
55 : :
56 : : contains
57 : : !!***
58 : :
59 : : !!****f* ABINIT/mklocl_realspace
60 : : !! NAME
61 : : !! mklocl_realspace
62 : : !!
63 : : !! FUNCTION
64 : : !! This method is equivalent to mklocl_recipspace except that
65 : : !! it uses real space pseudo-potentials. It is useful for isolated
66 : : !! systems. Then the option 3 and 4 are not available for this implementation.
67 : : !!
68 : : !! Optionally compute:
69 : : !! option=1 : local ionic potential throughout unit cell
70 : : !! option=2 : contribution of local ionic potential to E gradient wrt xred
71 : : !!
72 : : !! INPUTS
73 : : !! dtset <type(dataset_type)>=all input variables in this dataset
74 : : !! mpi_enreg=information about MPI parallelization
75 : : !! natom=number of atoms in unit cell.
76 : : !! nattyp(ntypat)=number of atoms of each type in cell.
77 : : !! nfft=(effective) number of FFT grid points (for this processor)
78 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
79 : : !! nspden=number of spin-density components
80 : : !! ntypat=number of types of atoms.
81 : : !! option= (see above)
82 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
83 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
84 : : !! rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$)
85 : : !! rhor(nfft,nspden)=electron density in electrons/bohr**3.
86 : : !! (needed if option==2 or if option==4)
87 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
88 : : !! ucvol=unit cell volume ($\textrm{Bohr}^3$).
89 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
90 : : !!
91 : : !! OUTPUT
92 : : !! (if option==1) vpsp(nfft)=local crystal pseudopotential in real space.
93 : : !! (if option==2) grtn(3,natom)=grads of Etot wrt tn. These gradients are in
94 : : !! reduced coordinates. Multiply them by rprimd to get
95 : : !! gradients in cartesian coordinates.
96 : : !!
97 : : !! SOURCE
98 : :
99 : 0 : subroutine mklocl_realspace(grtn,icoulomb,mpi_enreg,natom,nattyp,nfft,ngfft,nscforder, &
100 : 0 : & nspden,ntypat,option,pawtab,psps,rhog,rhor,rprimd,typat,&
101 : 0 : & ucvol,usewvl,vpsp,xred)
102 : :
103 : : #if defined HAVE_BIGDFT
104 : : use BigDFT_API, only : coulomb_operator,deallocate_coulomb_operator
105 : : use defs_PSolver
106 : : #else
107 : : use defs_wvltypes, only : coulomb_operator
108 : : #endif
109 : :
110 : : !Arguments ------------------------------------
111 : : !scalars
112 : : integer,intent(in) :: natom,nfft,nspden,ntypat,option
113 : : real(dp),intent(in) :: ucvol
114 : : type(MPI_type),intent(in) :: mpi_enreg
115 : : type(pseudopotential_type),intent(in) :: psps
116 : : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
117 : : !arrays
118 : : integer,intent(in) :: icoulomb,nscforder,usewvl
119 : : integer,intent(in) :: nattyp(ntypat),ngfft(18),typat(natom)
120 : : real(dp),intent(in) :: rhog(2,nfft)
121 : : real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3)
122 : : real(dp),intent(in) :: xred(3,natom)
123 : : real(dp),intent(out) :: grtn(3,natom),vpsp(nfft)
124 : :
125 : : !Local variables-------------------------------
126 : : character(len=1) :: geocode
127 : : !testing variables
128 : : !scalars
129 : : integer,parameter :: nStep=2
130 : : integer :: comm_fft,countParSeconde,i1,i2,i3
131 : : integer :: ia,ia1,ia2,igeo,ii,ind,itypat,ix,iy,iz,jj
132 : : integer :: kk,me_fft,n1,n2,n3,n3d,n_interpol
133 : : integer :: nproc_fft,tpsStart,tpsStop
134 : : real(dp),parameter :: min_rho_value=1.0d-12
135 : : real(dp) :: aa,bb,cc,dd,delta,deltaV,dr,dr2div6,invdr,r,vol_interpol,x,y,z,hgx,hgy,hgz,entmp
136 : : logical,parameter :: customRho=.false.,finiteDiff=.false.,testing=.false.
137 : : logical :: doIt
138 : : character(len=500) :: message
139 : : !arrays
140 : : integer :: ngfft_interpol(18)
141 : 0 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
142 : 0 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
143 : : real(dp) :: coord(3),coordXYZ(3),refValue(3),tsec(2)
144 : 0 : real(dp),allocatable :: coordCart_interpol(:,:),coordRed_interpol(:,:)
145 : 0 : real(dp),allocatable :: gridcart(:,:)
146 : 0 : real(dp),allocatable :: grtn_cart_interpol(:,:),grtn_diff(:,:)
147 : 0 : real(dp),allocatable :: rhog_interpol(:,:),rhog_testing(:,:),rhor_interpol(:)
148 : 0 : real(dp),allocatable :: rhor_testing(:),rhor_work(:),xcart(:,:),vhartr(:),gxyz(:,:)
149 : : type(coulomb_operator):: kernel
150 : :
151 : : ! *************************************************************************
152 : :
153 : : !Keep track of total time spent here
154 : 0 : if (option==2) then
155 : 0 : call timab(72,1,tsec)
156 : : end if
157 : :
158 : : !Several constants (FFT sizes and parallelism)
159 : 0 : n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
160 : 0 : nproc_fft = ngfft(10) ; me_fft = ngfft(11)
161 : 0 : n3d = ngfft(13) !for parallel runs
162 [ # # ]: 0 : if (nproc_fft==1) n3d=n3 !for serial runs
163 : 0 : comm_fft=mpi_enreg%comm_fft
164 [ # # # # ]: 0 : if(me_fft /= mpi_enreg%me_fft .or. nproc_fft /= mpi_enreg%nproc_fft) then
165 : 0 : ABI_BUG("mpi_enreg%x_fft not equal to the corresponding values in ngfft")
166 : : end if
167 : :
168 : : !Conditions for periodicity in the three directions
169 : 0 : geocode='P'
170 [ # # ]: 0 : if (icoulomb==1) geocode='F'
171 [ # # ]: 0 : if (icoulomb==2) geocode='S'
172 : :
173 : : !Get the distrib associated with this fft_grid
174 : 0 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
175 : :
176 : : !Store xcart for each atom
177 [ # # # # : 0 : ABI_MALLOC(xcart,(3, natom))
# # # # #
# # # ]
178 : 0 : call xred2xcart(natom, rprimd, xcart, xred)
179 : : !Store cartesian coordinates for each grid points
180 [ # # # # : 0 : ABI_MALLOC(gridcart,(3, nfft))
# # # # #
# # # ]
181 : 0 : call mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd)
182 : :
183 : : !Check whether all the PSP considered are of type GTH-HGH or PAW
184 : 0 : doIt=.true.
185 [ # # ]: 0 : do ii=1,psps%npsp
186 : : doIt=doIt .and.&
187 [ # # # # ]: 0 : (psps%pspcod(ii)==2.or.psps%pspcod(ii)==3.or.psps%pspcod(ii)==10.or.psps%pspcod(ii)==7)
188 : : end do
189 : :
190 : : !HGH-GTH/PAW treatment presumably starts here
191 [ # # ]: 0 : if (doIt) then
192 : :
193 : : ! Definition of the grid spacings as in the kernel routine
194 : 0 : hgx = rprimd(1,1)/(n1)
195 : 0 : hgy = rprimd(2,2)/(n2)
196 : 0 : hgz = rprimd(3,3)/(n3)
197 : :
198 : :
199 : : !----------------------------------------------------------------------
200 : : ! ----- Option 1: compute local ionic potential -----
201 : : !----------------------------------------------------------------------
202 [ # # ]: 0 : if (option==1) then
203 : :
204 : : call psolver_kernel( (/ hgx, hgy, hgz /), 2, icoulomb, me_fft, kernel, comm_fft, &
205 [ # # # # ]: 0 : & (/n1,n2,n3/), nproc_fft, nscforder)
206 : :
207 : : call createIonicPotential_new(fftn3_distrib,ffti3_local,&
208 : : & geocode,me_fft, nproc_fft, natom, &
209 : : & ntypat, typat, psps%gth_params%psppar, &
210 : : & int(psps%ziontypat), xcart,gridcart, hgx,hgy,hgz, &
211 [ # # # # : 0 : & n1,n2,n3d,n3, kernel, vpsp, comm_fft,pawtab,psps%usepaw)
# # ]
212 : :
213 : : !----------------------------------------------------------------------
214 : : ! ----- Option 2: compute forces induced by local ionic potential -----
215 : : !----------------------------------------------------------------------
216 [ # # ]: 0 : else if (option == 2) then
217 : :
218 : : ! Compute Hartree potential from rhor
219 [ # # # # : 0 : ABI_MALLOC(vhartr,(nfft))
# # # # #
# # # ]
220 : : call psolver_hartree(entmp, (/ hgx, hgy, hgz /), icoulomb, me_fft, comm_fft, nfft, &
221 [ # # # # ]: 0 : & (/n1,n2,n3/), nproc_fft, nscforder, nspden, rhor, vhartr, usewvl)
222 : :
223 : : ! Allocate temporary array for forces
224 [ # # # # : 0 : ABI_MALLOC(gxyz,(3, natom))
# # # # ]
225 : :
226 : : ! Calculate local part of the forces grtn (inspired from BigDFT routine)
227 : : call local_forces_new(fftn3_distrib,ffti3_local,geocode,me_fft, ntypat, natom, &
228 : : & typat, xcart, gridcart, psps%gth_params%psppar, int(psps%ziontypat), &
229 [ # # # # : 0 : & hgx,hgy,hgz, n1,n2,n3,n3d, rhor,vhartr, gxyz, pawtab,psps%usepaw)
# # ]
230 : :
231 : : ! Forces should be in reduced coordinates.
232 [ # # ]: 0 : do ia = 1, natom, 1
233 [ # # ]: 0 : do igeo = 1, 3, 1
234 : : grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) &
235 : : & - rprimd(2, igeo) * gxyz(2, ia) &
236 : 0 : & - rprimd(3, igeo) * gxyz(3, ia)
237 : : end do
238 : : end do
239 : :
240 : : ! Deallocate local variables
241 [ # # ]: 0 : ABI_FREE(vhartr)
242 [ # # ]: 0 : ABI_FREE(gxyz)
243 : : end if
244 : :
245 : : !----------------------------------------------------------------------
246 : : ! ----- Section for the non-HGH/GTH/PAW pseudopotentials (testing) ----
247 : : !----------------------------------------------------------------------
248 : : else
249 : :
250 : : if (testing) then
251 : : call system_clock(count_rate = countParSeconde)
252 : : call system_clock(tpsStart, count_rate = countParSeconde)
253 : : end if
254 : :
255 : : ! dr is the r step in the sampling psps%vlspl
256 : 0 : dr = psps%qgrid_vl(2)
257 : 0 : invdr = 1._dp / dr
258 : 0 : dr2div6 = dr * dr / 6._dp
259 : :
260 [ # # ]: 0 : if (option == 1) then
261 : : ! Set 0 in vpsp before summing
262 [ # # ]: 0 : vpsp(:) = 0._dp
263 [ # # ]: 0 : else if (option == 2) then
264 : : ! Allocate array to store cartesian gradient computed with
265 : : ! an interpolation of rhor
266 [ # # # # : 0 : ABI_MALLOC(grtn_cart_interpol,(3, natom))
# # # # ]
267 [ # # # # ]: 0 : grtn_cart_interpol(:, :) = 0._dp
268 : :
269 : 0 : n_interpol = nStep ** 3
270 [ # # ]: 0 : ABI_MALLOC(coordRed_interpol,(3, nStep ** 3))
271 [ # # ]: 0 : ABI_MALLOC(coordCart_interpol,(3, nStep ** 3))
272 : :
273 : : if (testing .and. customRho) then
274 : : ! Use a custom rho instead of the self-consistent one.
275 : : ABI_MALLOC(rhor_testing,(nfft))
276 : : ABI_MALLOC(rhog_testing,(2, nfft))
277 : : end if
278 : :
279 [ # # # # : 0 : ABI_MALLOC(rhor_interpol,(nfft * n_interpol))
# # # # #
# # # ]
280 [ # # # # : 0 : ABI_MALLOC(rhor_work,(nfft * n_interpol))
# # # # ]
281 [ # # # # : 0 : ABI_MALLOC(rhog_interpol,(2, nfft * n_interpol))
# # # # #
# # # ]
282 : :
283 : : if (testing .and. customRho) then
284 : : ! Testing only, changing rho with a centered gaussian
285 : : do ii = 1, nfft, 1
286 : : ! using the position of the first atom as center.
287 : : r = (gridcart(1, ii) - xcart(1, 1)) ** 2 + &
288 : : & (gridcart(2, ii) - xcart(2, 1)) ** 2 + &
289 : : & (gridcart(3, ii) - xcart(3, 1)) ** 2
290 : : rhor_testing(ii) = exp(-r/4._dp)
291 : : end do
292 : : ! Testing only, compute rhog_testing from rhor_testing
293 : : call fourdp(1,rhog_testing,rhor_testing,-1,mpi_enreg,nfft,1,ngfft,0)
294 : : end if
295 : :
296 : : ! Compute the interpolation of rho, using a fourier transform
297 [ # # # # ]: 0 : rhog_interpol(:, :) = 0._dp
298 : : ii = 0
299 [ # # ]: 0 : do i3 = 1, n3, 1
300 [ # # ]: 0 : if (i3 <= n3 / 2) then
301 : : iz = i3
302 : : else
303 : 0 : iz = n3 * nStep - n3 + i3
304 : : end if
305 [ # # ]: 0 : do i2 = 1, n2, 1
306 [ # # ]: 0 : if (i2 <= n2 / 2) then
307 : : iy = i2
308 : : else
309 : 0 : iy = n2 * nStep - n2 + i2
310 : : end if
311 [ # # ]: 0 : do i1 = 1, n1, 1
312 : 0 : ii = ii + 1
313 [ # # ]: 0 : if (i1 <= n1 / 2) then
314 : : ix = i1
315 : : else
316 : 0 : ix = n1 * nStep - n1 + i1
317 : : end if
318 : 0 : jj = (iz - 1) * n2 * n1 * nStep ** 2 + (iy - 1) * n3 * nStep + ix
319 : 0 : if (testing .and. customRho) then
320 : : rhog_interpol(:, jj) = rhog_testing(:, ii)
321 : : else
322 [ # # ]: 0 : rhog_interpol(:, jj) = rhog(:, ii)
323 : : end if
324 : : end do
325 : : end do
326 : : end do
327 : :
328 : : ! Compute the interpolation of rho from the Fourier transformation
329 : 0 : ngfft_interpol(:) = ngfft(:)
330 [ # # ]: 0 : ngfft_interpol(1:3) = (/ n1 * nStep, n2 * nStep, n3 * nStep /)
331 [ # # ]: 0 : ngfft_interpol(4:6) = (/ n1 * nStep + 1, n2 * nStep + 1, n3 * nStep /)
332 : 0 : call fourdp(1,rhog_interpol,rhor_work,1,mpi_enreg,nfft*n_interpol,1,ngfft_interpol,0)
333 : :
334 : : ! Reorder rhor_interpol to be able to read it linearly
335 : 0 : jj = 0
336 [ # # ]: 0 : do i3 = 1, n3, 1
337 [ # # ]: 0 : do i2 = 1, n2, 1
338 [ # # ]: 0 : do i1 = 1, n1, 1
339 [ # # ]: 0 : do iz = 1, nStep, 1
340 [ # # ]: 0 : do iy = 1, nStep, 1
341 [ # # ]: 0 : do ix = 1, nStep, 1
342 : 0 : jj = jj + 1
343 : 0 : kk = ((i3 - 1) * nStep + iz - 1) ! z coordinate in the interpolated grid
344 : 0 : kk = kk * n1 * n2 * nStep ** 2
345 : 0 : kk = kk + ((i2 - 1) * nStep + iy - 1) * n1 * nStep ! adding y coordinate
346 : 0 : kk = kk + (i1 - 1) * nStep + ix ! adding x coordinate
347 : 0 : rhor_interpol(jj) = rhor_work(kk)
348 : : end do
349 : : end do
350 : : end do
351 : : end do
352 : : end do
353 : : end do
354 [ # # ]: 0 : ABI_FREE(rhor_work)
355 : :
356 : : ! Compute grid access in the interpolated volume
357 : 0 : ii = 0
358 [ # # ]: 0 : do iz = 1, nStep, 1
359 : 0 : z = real(iz - 1, dp) / real(nStep, dp)
360 [ # # ]: 0 : do iy = 1, nStep, 1
361 : 0 : y = real(iy - 1, dp) / real(nStep, dp)
362 [ # # ]: 0 : do ix = 1, nStep, 1
363 : 0 : x = real(ix - 1, dp) / real(nStep, dp)
364 : 0 : ii = ii + 1
365 [ # # ]: 0 : coordRed_interpol(:, ii) = (/ x, y, z /)
366 : : ! Assuming orthogonal box (should be change later)
367 : : coordCart_interpol(:, ii) = (/ x * rprimd(1, 1) / real(n1, dp), &
368 : : & y * rprimd(2, 2) / real(n2, dp), &
369 [ # # ]: 0 : & z * rprimd(3, 3) / real(n3, dp) /)
370 : : end do
371 : : end do
372 : : end do
373 : :
374 : : vol_interpol = 1._dp / real(nStep, dp) ** 3
375 : : ! Compute the coordinates (integer) of each atom and deduce
376 : : ! the max extens of the integral summation.
377 : : ! !! do ia = 1, natom, 1
378 : : ! !! coordAtom(1, ia) = int(xred(1, ia) * n1) + 1
379 : : ! !! coordAtom(2, ia) = int(xred(2, ia) * n2) + 1
380 : : ! !! coordAtom(3, ia) = int(xred(3, ia) * n3) + 1
381 : : ! !! end do
382 : : end if
383 : :
384 : : if (testing .and. option == 2) then
385 : : call system_clock(tpsStop, count_rate = countParSeconde)
386 : : write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde)
387 : : end if
388 : :
389 : 0 : ia1=1
390 [ # # ]: 0 : do itypat = 1, ntypat, 1
391 : : ! ia1,ia2 sets range of loop over atoms:
392 : 0 : ia2 = ia1 + nattyp(itypat) - 1
393 : :
394 [ # # ]: 0 : do ii = 1, nfft, 1
395 [ # # ]: 0 : do ia = ia1, ia2, 1
396 [ # # ]: 0 : if (option == 1) then
397 : : ! Compute the potential
398 : : ! r is the distance between grid point and atom
399 : : r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + &
400 : : & (gridcart(2, ii) - xcart(2, ia)) ** 2 + &
401 : 0 : & (gridcart(3, ii) - xcart(3, ia)) ** 2)
402 : :
403 : : ! Coefficients needed to compute the spline.
404 : 0 : jj = int(r * invdr) + 1
405 [ # # ]: 0 : if (jj > psps%mqgrid_vl - 2) then
406 : : write(message, '(3a,i0,a,i0,a,a)' )&
407 : 0 : & ' pseudo-potential local part sampling is not wide enough', ch10, &
408 : 0 : & ' want to access position ', jj, ' whereas mqgrid_vl = ', psps%mqgrid_vl, ch10, &
409 : 0 : & ' Action : no idea, contact developpers...'
410 : 0 : ABI_ERROR(message)
411 : : end if
412 : 0 : delta = r - psps%qgrid_vl(jj)
413 : 0 : bb = delta * invdr
414 : 0 : aa = 1._dp - bb
415 : 0 : cc = aa * (aa ** 2 - 1._dp) * dr2div6
416 : 0 : dd = bb * (bb ** 2 - 1._dp) * dr2div6
417 : :
418 : : ! compute V(r) from the spline, jj and jj + 1 is braketting r in
419 : : ! the sampling
420 : : deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + &
421 : 0 : & cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat)
422 : : ! Add on grid point ii the contribution of atom ia
423 : 0 : vpsp(ii) = vpsp(ii) + deltaV
424 [ # # ]: 0 : else if (option == 2) then
425 : : ! Compute the forces, as gradient of energy (V(r).rho(r))
426 : :
427 : : ! Testing only - reference points
428 : : if (.false.) then
429 : : ! r is the distance between grid point and atom
430 : : r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + &
431 : : & (gridcart(2, ii) - xcart(2, ia)) ** 2 + &
432 : : & (gridcart(3, ii) - xcart(3, ia)) ** 2)
433 : :
434 : : ! Coefficients needed to compute the spline.
435 : : jj = int(r * invdr) + 1
436 : : delta = r - psps%qgrid_vl(jj)
437 : : bb = delta * invdr
438 : : aa = 1._dp - bb
439 : : cc = aa * (aa ** 2 - 1._dp) * dr2div6
440 : : dd = bb * (bb ** 2 - 1._dp) * dr2div6
441 : :
442 : : ! When mesh position is on a node, forces are null.
443 : : if (r /= 0._dp) then
444 : : ! This value deltaV is the first derivative of V(r) taken at r.
445 : : deltaV = aa * psps%dvlspl(jj, 1, itypat) + bb * psps%dvlspl(jj + 1, 1, itypat) + &
446 : : & cc * psps%dvlspl(jj, 2, itypat) + dd * psps%dvlspl(jj + 1, 2, itypat)
447 : : ! We multiply by rho(r) to have an energy.
448 : : deltaV = deltaV * rhor(ii, 1) / r
449 : : refValue(:) = - deltaV * (gridcart(:, ii) - xcart(:, ia))
450 : : grtn_cart_interpol(:, ia) = grtn_cart_interpol(:, ia) + refValue(:)
451 : : end if
452 : : end if
453 : :
454 : : ! Compute the interpolation for the point ii
455 : 0 : ind = (ii - 1) * n_interpol
456 [ # # ]: 0 : do kk = 1, n_interpol, 1
457 : 0 : ind = ind + 1
458 : :
459 [ # # ]: 0 : if (rhor_interpol(ind) > min_rho_value) then
460 : : ! Assume orthogonal box...
461 : 0 : coordXYZ(1) = gridcart(1, ii) - xcart(1, ia) + coordCart_interpol(1, kk)
462 : 0 : coordXYZ(2) = gridcart(2, ii) - xcart(2, ia) + coordCart_interpol(2, kk)
463 : 0 : coordXYZ(3) = gridcart(3, ii) - xcart(3, ia) + coordCart_interpol(3, kk)
464 : 0 : r = coordXYZ(1) ** 2 + coordXYZ(2) ** 2 + coordXYZ(3) ** 2
465 : :
466 [ # # ]: 0 : if (r /= 0._dp) then
467 : 0 : r = sqrt(r)
468 : : ! Coefficients needed to compute the spline.
469 : 0 : jj = int(r * invdr) + 1
470 : 0 : delta = r - psps%qgrid_vl(jj)
471 : 0 : bb = delta * invdr
472 : 0 : aa = 1._dp - bb
473 : 0 : cc = aa * (aa ** 2 - 1._dp) * dr2div6
474 : 0 : dd = bb * (bb ** 2 - 1._dp) * dr2div6
475 : : deltaV = aa * psps%dvlspl(jj, 1, itypat) + &
476 : : & bb * psps%dvlspl(jj + 1, 1, itypat) + &
477 : : & cc * psps%dvlspl(jj, 2, itypat) + &
478 : 0 : & dd * psps%dvlspl(jj + 1, 2, itypat)
479 : 0 : deltaV = deltaV * rhor_interpol(ind) / r
480 : 0 : grtn_cart_interpol(1, ia) = grtn_cart_interpol(1, ia) - deltaV * coordXYZ(1)
481 : 0 : grtn_cart_interpol(2, ia) = grtn_cart_interpol(2, ia) - deltaV * coordXYZ(2)
482 : 0 : grtn_cart_interpol(3, ia) = grtn_cart_interpol(3, ia) - deltaV * coordXYZ(3)
483 : : ! do igeo = 1, 3, 1
484 : : ! grtn_cart_interpol(igeo, ia) = grtn_cart_interpol(igeo, ia) - deltaV * coordXYZ(igeo)
485 : : ! end do
486 : : end if
487 : : end if
488 : : end do
489 : :
490 : : ! =============
491 : : ! Testing only
492 : : ! =============
493 : : ! use of finite differences
494 : : if (finiteDiff) then
495 : : do igeo = 1, 3, 1
496 : : coord(:) = 0._dp
497 : : coord(igeo) = dr / 2.0_dp
498 : : r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + &
499 : : & (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + &
500 : : & (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2)
501 : :
502 : : ! Coefficients needed to compute the spline.
503 : : jj = int(r * invdr) + 1
504 : : delta = r - psps%qgrid_vl(jj)
505 : : bb = delta * invdr
506 : : aa = 1._dp - bb
507 : : cc = aa * (aa ** 2 - 1._dp) * dr2div6
508 : : dd = bb * (bb ** 2 - 1._dp) * dr2div6
509 : :
510 : : deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + &
511 : : & cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat)
512 : :
513 : :
514 : : coord(:) = 0._dp
515 : : coord(igeo) = -dr / 2.0_dp
516 : : r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + &
517 : : & (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + &
518 : : & (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2)
519 : :
520 : : ! Coefficients needed to compute the spline.
521 : : jj = int(r * invdr) + 1
522 : : delta = r - psps%qgrid_vl(jj)
523 : : bb = delta * invdr
524 : : aa = 1._dp - bb
525 : : cc = aa * (aa ** 2 - 1._dp) * dr2div6
526 : : dd = bb * (bb ** 2 - 1._dp) * dr2div6
527 : :
528 : : deltaV = deltaV - (aa * psps%vlspl(jj, 1, itypat) + &
529 : : & bb * psps%vlspl(jj + 1, 1, itypat) + &
530 : : & cc * psps%vlspl(jj, 2, itypat) + &
531 : : & dd * psps%vlspl(jj + 1, 2, itypat))
532 : : grtn_diff(igeo, ia) = grtn_diff(igeo, ia) - deltaV * rhor(ii, 1) / dr
533 : : end do
534 : : end if
535 : : ! =============
536 : : ! Testing only
537 : : ! =============
538 : :
539 : : end if
540 : : end do
541 : : ! End loop over atoms of type itypat
542 : : end do
543 : : ! End loop over real space grid points
544 : :
545 : 0 : ia1 = ia2 + 1
546 : : end do
547 : : ! End loop over type of atoms
548 : :
549 [ # # ]: 0 : if(option==2)then
550 : : ! multiply the forces by the volume of a single box mesh.
551 : : grtn_cart_interpol(:, :) = grtn_cart_interpol(:, :) * &
552 [ # # # # ]: 0 : & ucvol / real(n1 * n2 * n3, dp) * vol_interpol
553 : : ! Transform cartesian forces to reduce coordinates
554 [ # # ]: 0 : do ia = 1, natom, 1
555 [ # # ]: 0 : do igeo = 1, 3, 1
556 : : grtn(igeo, ia) = rprimd(1, igeo) * grtn_cart_interpol(1, ia) + &
557 : : & rprimd(2, igeo) * grtn_cart_interpol(2, ia) + &
558 : 0 : & rprimd(3, igeo) * grtn_cart_interpol(3, ia)
559 : : end do
560 : : end do
561 [ # # ]: 0 : ABI_FREE(rhor_interpol)
562 [ # # ]: 0 : ABI_FREE(rhog_interpol)
563 [ # # ]: 0 : ABI_FREE(coordRed_interpol)
564 [ # # ]: 0 : ABI_FREE(coordCart_interpol)
565 : : if (testing .and. customRho) then
566 : : ABI_FREE(rhor_testing)
567 : : ABI_FREE(rhog_testing)
568 : : end if
569 : :
570 : : if (testing) then
571 : : call system_clock(tpsStop, count_rate = countParSeconde)
572 : : write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde)
573 : : write(std_out,*) grtn_cart_interpol
574 : : ABI_ERROR("Testing section!")
575 : : end if
576 : :
577 : : end if
578 : :
579 : : !-------------------------
580 : : end if ! GTH/HGH/PAW psps
581 : :
582 : : !Release temporary memory
583 [ # # ]: 0 : ABI_FREE(xcart)
584 [ # # ]: 0 : ABI_FREE(gridcart)
585 : :
586 : : !Close timing counters
587 [ # # ]: 0 : if (option==2)then
588 : 0 : call timab(72,2,tsec)
589 : : end if
590 : :
591 [ # # # # : 0 : end subroutine mklocl_realspace
# # # # #
# # # ]
592 : : !!***
593 : :
594 : : !----------------------------------------------------------------------
595 : :
596 : : !!****f* mklocl_realspace/createIonicPotential_new
597 : : !! NAME
598 : : !! createIonicPotential_new
599 : : !!
600 : : !! FUNCTION
601 : : !!
602 : : !! INPUTS
603 : : !!
604 : : !! OUTPUT
605 : : !!
606 : : !! SOURCE
607 : :
608 : 0 : subroutine createIonicPotential_new(fftn3_distrib,ffti3_local,geocode,iproc,&
609 : 0 : & nproc,nat,ntypes,iatype,psppar,nelpsp,rxyz,gridcart,&
610 : 0 : & hxh,hyh,hzh,n1i,n2i,n3d,n3i,kernel,pot_ion,spaceworld,pawtab,usepaw)
611 : :
612 : : use defs_wvltypes, only : coulomb_operator
613 : : #if defined HAVE_BIGDFT
614 : : use BigDFT_API, only : xc_info,xc_init,xc_end,XC_ABINIT
615 : : #endif
616 : :
617 : : !Arguments -------------------------------
618 : : !scalars
619 : : integer, intent(in) :: iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw
620 : : real(dp), intent(in) :: hxh,hyh,hzh
621 : : character(len=1), intent(in) :: geocode
622 : : type(coulomb_operator), intent(in) :: kernel
623 : : !arrays
624 : : integer, dimension(nat), intent(in) :: iatype
625 : : integer, dimension(ntypes), intent(in) :: nelpsp
626 : : integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local
627 : : real(dp), dimension(3,n1i*n2i*n3d), intent(in) :: gridcart
628 : : real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar
629 : : real(dp), dimension(3,nat), intent(in) :: rxyz
630 : : real(dp), dimension(*), intent(inout) :: pot_ion
631 : : type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw)
632 : :
633 : : !Local variables -------------------------
634 : : #if defined HAVE_BIGDFT
635 : : !scalars
636 : : integer :: iat,i1,i2,i3,j1,j2,j3,isx,isy,isz,iex,iey,iez,ierr,ityp
637 : : integer :: ind,nloc,iloc,i3loc,msz
638 : : logical :: gox,goy,goz,perx,pery,perz
639 : : real(dp) :: arg,charge,cutoff,ehart,eexcu,rholeaked,rholeaked_tot,rloc
640 : : real(dp) :: rx,ry,rz,rzero,r2,tt,tt_tot,vexcu,vhgh,x,xp,y,z
641 : : type(xc_info) :: xc
642 : : !arrays
643 : : real(dp) :: rr(1),vpaw(1)
644 : : real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
645 : : #endif
646 : :
647 : : ! *********************************************************************
648 : :
649 : : #if defined HAVE_BIGDFT
650 : :
651 : : if(nproc<0)then
652 : : ABI_ERROR('nproc should not be negative')
653 : : end if
654 : :
655 : : !Ionic charge (must be calculated for the PS active processes)
656 : : rholeaked=0._dp
657 : : !Ionic energy (can be calculated for all the processors)
658 : :
659 : : !here we should insert the calculation of the ewald energy for the periodic BC case
660 : : !!! eion=0._dp
661 : : !!! do iat=1,nat
662 : : !!! ityp=iatype(iat)
663 : : !!! rx=rxyz(1,iat)
664 : : !!! ry=rxyz(2,iat)
665 : : !!! rz=rxyz(3,iat)
666 : : !!! ! ion-ion interaction
667 : : !!! do jat=1,iat-1
668 : : !!! dist=sqrt( (rx-rxyz(1,jat))**2+(ry-rxyz(2,jat))**2+(rz-rxyz(3,jat))**2 )
669 : : !!! jtyp=iatype(jat)
670 : : !!! eion=eion+real(nelpsp(jtyp)*nelpsp(ityp),kind=dp)/dist
671 : : !!! enddo
672 : : !!! end do
673 : : !!! if (iproc.eq.0) write(std_out,'(1x,a,1pe22.14)') 'ion-ion interaction energy',eion
674 : :
675 : : !Creates charge density arising from the ionic PSP cores
676 : : !the n3pi dimension indicates the number of planes trated by each processor in the FFT parallelisation
677 : : !for a plane wave treatment this value depends on whether the direct space is divided in planes or not
678 : : !I don't know this variable, which in the future must be inserted at the place of n3pi (LG)
679 : : !if n3pi=0 this means that the processors doesn't calculate anything
680 : : !if (n3pi >0 ) then
681 : :
682 : : !conditions for periodicity in the three directions
683 : : perx=(geocode /= 'F')
684 : : pery=(geocode == 'P')
685 : : perz=(geocode /= 'F')
686 : :
687 : : !this initialise the array to zero, it will work only if bigdft library is enabled
688 : : pot_ion(1:n1i*n2i*n3d)=zero
689 : :
690 : : do iat=1,nat
691 : : ityp=iatype(iat)
692 : : rx=rxyz(1,iat)
693 : : ry=rxyz(2,iat)
694 : : rz=rxyz(3,iat)
695 : :
696 : : rloc=psppar(0,0,ityp)
697 : : charge=real(nelpsp(ityp),kind=dp)/(2._dp*pi*sqrt(2._dp*pi)*rloc**3)
698 : : cutoff=10._dp*rloc
699 : :
700 : : isx=floor((rx-cutoff)/hxh)
701 : : isy=floor((ry-cutoff)/hyh)
702 : : isz=floor((rz-cutoff)/hzh)
703 : :
704 : : iex=ceiling((rx+cutoff)/hxh)
705 : : iey=ceiling((ry+cutoff)/hyh)
706 : : iez=ceiling((rz+cutoff)/hzh)
707 : :
708 : : ! Calculate Ionic Density
709 : : ! using HGH parameters.
710 : : ! Eq. 1.104, T. Deutsch and L. Genovese, JDN. 12, 2011
711 : : do i3=isz,iez
712 : : z=real(i3,kind=dp)*hzh-rz
713 : : call ind_positions_mklocl(perz,i3,n3i,j3,goz)
714 : : if(fftn3_distrib(j3)==iproc) then
715 : : i3loc=ffti3_local(j3)
716 : : do i2=isy,iey
717 : : y=real(i2,kind=dp)*hyh-ry
718 : : call ind_positions_mklocl(pery,i2,n2i,j2,goy)
719 : : do i1=isx,iex
720 : : x=real(i1,kind=dp)*hxh-rx
721 : : call ind_positions_mklocl(perx,i1,n1i,j1,gox)
722 : : r2=x**2+y**2+z**2
723 : : if (goz .and. goy .and. gox ) then
724 : : ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i
725 : : r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2
726 : : end if
727 : : arg=r2/rloc**2
728 : : xp=exp(-.5d0*arg)
729 : : if (goz .and. goy .and. gox ) then
730 : : pot_ion(ind)=pot_ion(ind)-xp*charge
731 : : else
732 : : rholeaked=rholeaked+xp*charge
733 : : end if
734 : : end do
735 : : end do
736 : : end if
737 : : end do
738 : :
739 : : end do
740 : :
741 : : !Check
742 : : tt=0._dp
743 : : do j3= 1,n3d
744 : : do i2= 1,n2i
745 : : do i1= 1,n1i
746 : : ind=i1+(i2-1)*n1i+(j3-1)*n1i*n2i
747 : : tt=tt+pot_ion(ind)
748 : : end do
749 : : end do
750 : : end do
751 : :
752 : : tt=tt*hxh*hyh*hzh
753 : : rholeaked=rholeaked*hxh*hyh*hzh
754 : :
755 : : call xmpi_sum(tt,tt_tot,spaceworld,ierr)
756 : : call xmpi_sum(rholeaked,rholeaked_tot,spaceworld,ierr)
757 : :
758 : : if (iproc.eq.0) then
759 : : write(std_out,'(1x,a,f26.12,2x,1pe10.3)') &
760 : : & 'total ionic charge, leaked charge ',tt_tot,rholeaked_tot
761 : : end if
762 : :
763 : : !Here the value of the datacode must be kept fixed
764 : : !there can be some problems when running this stuff in parallel,
765 : : ! if the ionic potential distribution does not agree with the
766 : : ! plane distribution which is supposed to hold for the Poisson Solver
767 : : !No exchange-correlation here: PSolver expects an xc_info object, set to ixc=0
768 : : call xc_init(xc,0,XC_ABINIT,1)
769 : : call psolver(geocode,'D',iproc,nproc,n1i,n2i,n3i,xc,hxh,hyh,hzh,&
770 : : & pot_ion,kernel%kernel,pot_ion,ehart,eexcu,vexcu,0._dp,.false.,1)
771 : : call xc_end(xc)
772 : :
773 : : !Add the remaining short-range local terms
774 : : do iat=1,nat
775 : : ityp=iatype(iat)
776 : :
777 : : rx=rxyz(1,iat)
778 : : ry=rxyz(2,iat)
779 : : rz=rxyz(3,iat)
780 : :
781 : : ! determine number of local terms
782 : : rloc=psppar(0,0,ityp)
783 : : cutoff=10._dp*rloc
784 : : charge=real(nelpsp(ityp),kind=dp)
785 : :
786 : : ! determine number of local terms (HGH pot)
787 : : nloc=0
788 : : do iloc=1,4
789 : : if (psppar(0,iloc,ityp).ne.0._dp) nloc=iloc
790 : : end do
791 : :
792 : : ! PAW specifics
793 : : if (usepaw==1) then
794 : : msz=pawtab(ityp)%wvl%rholoc%msz
795 : : rad => pawtab(ityp)%wvl%rholoc%rad(1:msz)
796 : : vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
797 : : d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
798 : : rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2)
799 : : end if
800 : :
801 : : isx=floor((rx-cutoff)/hxh)
802 : : isy=floor((ry-cutoff)/hyh)
803 : : isz=floor((rz-cutoff)/hzh)
804 : :
805 : : iex=ceiling((rx+cutoff)/hxh)
806 : : iey=ceiling((ry+cutoff)/hyh)
807 : : iez=ceiling((rz+cutoff)/hzh)
808 : :
809 : : do i3=isz,iez
810 : : z=real(i3,kind=dp)*hzh-rz
811 : : call ind_positions_mklocl(perz,i3,n3i,j3,goz)
812 : : if(fftn3_distrib(j3) == iproc .and. goz) then !MPI
813 : : i3loc=ffti3_local(j3)
814 : : if (goz) then
815 : : do i2=isy,iey
816 : : y=real(i2,kind=dp)*hyh-ry
817 : : call ind_positions_mklocl(pery,i2,n2i,j2,goy)
818 : : if (goy) then
819 : : do i1=isx,iex
820 : : x=real(i1,kind=dp)*hxh-rx
821 : : call ind_positions_mklocl(perx,i1,n1i,j1,gox)
822 : : if (gox) then
823 : : ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i
824 : : r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2
825 : : ! r2=x**2+y**2+z**2
826 : :
827 : : ! HGH: V_S=gaussian potential of Eq. (9) in JCP 129, 014109(2008)
828 : : if (usepaw==0) then
829 : : if (nloc /= 0) then
830 : : arg=r2/rloc**2
831 : : xp=exp(-.5d0*arg)
832 : : tt=psppar(0,nloc,ityp)
833 : : do iloc=nloc-1,1,-1
834 : : tt=arg*tt+psppar(0,iloc,ityp)
835 : : end do
836 : : pot_ion(ind)=pot_ion(ind)+xp*tt
837 : : end if
838 : :
839 : : ! PAW: V_PAW-V_L^HGH
840 : : else
841 : : rr(1)=sqrt(r2)
842 : : if (rr(1)>=rzero) then
843 : : call paw_splint(msz,rad,vloc,d2vloc,1,rr,vpaw)
844 : : call calcVloc_mklocl(vhgh,rr(1),rloc,charge)
845 : : pot_ion(ind)=pot_ion(ind)+vpaw(1)-vhgh
846 : : else
847 : : pot_ion(ind)=pot_ion(ind)+vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
848 : : end if
849 : : end if
850 : :
851 : : end if
852 : : end do
853 : : end if
854 : : end do
855 : : end if
856 : : end if
857 : : end do
858 : :
859 : : end do !iat
860 : : #else
861 : 0 : BIGDFT_NOTENABLED_ERROR()
862 : : if (.false.) write(std_out,*) geocode,iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw,&
863 : : & hxh,hyh,hzh,iatype(1),nelpsp(1),fftn3_distrib(1),ffti3_local(1),gridcart(1,1),psppar(1,1,1),&
864 : : & rxyz(1,1),pot_ion(1),pawtab(1)%mesh_size,kernel%co
865 : : #endif
866 : :
867 : : CONTAINS
868 : : !!***
869 : :
870 : : !----------------------------------------------------------------------
871 : :
872 : : !!****f* mklocl_realspace/calcVloc_mklocl
873 : : !! NAME
874 : : !! calcVloc_mklocl
875 : : !!
876 : : !! FUNCTION
877 : : !!
878 : : !! INPUTS
879 : : !!
880 : : !! OUTPUT
881 : : !!
882 : : !! SOURCE
883 : :
884 : : subroutine calcVloc_mklocl(yy,xx,rloc,Z)
885 : :
886 : : !Arguments ------------------------------------
887 : : !scalars
888 : : real(dp),intent(in) :: xx,rloc,Z
889 : : real(dp),intent(out) :: yy
890 : :
891 : : !Local variables-------------------------------
892 : : !scalars
893 : : real(dp):: arg,tt
894 : :
895 : : ! *************************************************************************
896 : :
897 : : arg=xx/(sqrt(2.0)*rloc)
898 : : call derf_ab(tt,arg)
899 : : yy=-Z/xx*tt
900 : :
901 : : end subroutine calcVloc_mklocl
902 : : !!***
903 : :
904 : : !----------------------------------------------------------------------
905 : :
906 : : !!****f* mklocl_realspace/vloc_zero_mklocl
907 : : !! NAME
908 : : !! vloc_zero_mklocl
909 : : !!
910 : : !! FUNCTION
911 : : !! Use a quadratic interpolation to get limit of Vloc(x) at x->0
912 : : !!
913 : : !! INPUTS
914 : : !!
915 : : !! OUTPUT
916 : : !!
917 : : !! SOURCE
918 : :
919 : : function vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
920 : :
921 : :
922 : : !Arguments ------------------------------------
923 : : !scalars
924 : : integer,intent(in) :: msz
925 : : real(dp) :: vloc_zero_mklocl
926 : : real(dp),intent(in) :: charge,rloc
927 : : !arrays
928 : : real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
929 : :
930 : : !Local variables-------------------------------
931 : : !scalars
932 : : real(dp) :: y1,y2,y3,zz=0._dp
933 : : !arrays
934 : : real(dp) :: ll(3),xx(3),yy(3)
935 : :
936 : : ! *************************************************************************
937 : :
938 : : !Select 3 points x1,x2,x3 near 0
939 : : if (rad(1)>1.d-10) then
940 : : xx(1:3)=rad(1:3)
941 : : else
942 : : xx(1:3)=rad(2:4)
943 : : end if
944 : :
945 : : !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
946 : : call paw_splint(msz,rad,vloc,d2vloc,3,xx,yy)
947 : : call calcVloc_mklocl(y1,xx(1),rloc,charge)
948 : : call calcVloc_mklocl(y2,xx(2),rloc,charge)
949 : : call calcVloc_mklocl(y3,xx(3),rloc,charge)
950 : : yy(1)= yy(1)-y1
951 : : yy(2)= yy(2)-y2
952 : : yy(3)= yy(3)-y3
953 : :
954 : : !Find a polynomial of the form (z=0):
955 : : !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
956 : :
957 : : !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
958 : : ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
959 : : !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
960 : : ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
961 : : !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
962 : : ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
963 : :
964 : : vloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
965 : :
966 : : end function vloc_zero_mklocl
967 : : !!***
968 : :
969 : : end subroutine createIonicPotential_new
970 : : !!***
971 : :
972 : : !----------------------------------------------------------------------
973 : :
974 : : !!****f* mklocl_realspace/local_forces_new
975 : : !! NAME
976 : : !! local_forces_new
977 : : !!
978 : : !! FUNCTION
979 : : !!
980 : : !! INPUTS
981 : : !!
982 : : !! OUTPUT
983 : : !!
984 : : !! SOURCE
985 : :
986 : 0 : subroutine local_forces_new(fftn3_distrib,ffti3_local,&
987 : 0 : geocode,iproc,ntypes,nat,iatype,rxyz,gridcart,psppar,nelpsp,hxh,hyh,hzh,&
988 : 0 : n1,n2,n3,n3d,rho,pot,floc,pawtab,usepaw)
989 : :
990 : :
991 : : !Arguments -------------------------------
992 : : !scalars
993 : : integer, intent(in) :: iproc,ntypes,nat,n1,n2,n3,n3d,usepaw
994 : : character(len=1), intent(in) :: geocode
995 : : real(dp), intent(in) :: hxh,hyh,hzh
996 : : !arrays
997 : : integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local
998 : : integer, dimension(nat), intent(in) :: iatype
999 : : integer, dimension(ntypes), intent(in) :: nelpsp
1000 : : real(dp), dimension(3,n1*n2*n3d), intent(in) :: gridcart
1001 : : real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar
1002 : : real(dp), dimension(3,nat), intent(in) :: rxyz
1003 : : real(dp), dimension(*), intent(in) :: rho,pot
1004 : : real(dp), dimension(3,nat), intent(out) :: floc
1005 : : type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw)
1006 : :
1007 : : !Local variables -------------------------
1008 : : !scalars
1009 : : integer :: isx,isy,isz,iex,iey,iez,i1,i2,i3,j1,j2,j3,ind,iat,ityp,iloc,i3loc,msz,nloc
1010 : : logical :: perx,pery,perz,gox,goy,goz
1011 : : real(dp) :: arg,charge,cutoff,dvhgh,fxerf,fyerf,fzerf,fxgau,fygau,fzgau,forceleaked
1012 : : real(dp) :: forceloc,prefactor,rloc,rloc2,rhoel,rx,ry,rz,rzero,r2,tt,x,xp,y,z,Vel
1013 : : real(dp), dimension(4) :: cprime
1014 : : !arrays
1015 : : real(dp) :: dvpawdr(1),rr(1)
1016 : 0 : real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
1017 : :
1018 : : ! *********************************************************************
1019 : :
1020 : 0 : if (iproc == 0) write(std_out,'(1x,a)',advance='no')'Calculate local forces...'
1021 : :
1022 : : !Conditions for periodicity in the three directions
1023 : 0 : perx=(geocode /= 'F')
1024 : 0 : pery=(geocode == 'P')
1025 : 0 : perz=(geocode /= 'F')
1026 : :
1027 : 0 : forceleaked=zero
1028 : :
1029 [ # # ]: 0 : do iat=1,nat
1030 : 0 : ityp=iatype(iat)
1031 : : ! Coordinates of the center
1032 : 0 : rx=rxyz(1,iat)
1033 : 0 : ry=rxyz(2,iat)
1034 : 0 : rz=rxyz(3,iat)
1035 : :
1036 : : ! Initialization of the forces
1037 : : ! ion-electron term, error function part
1038 : 0 : fxerf=zero
1039 : 0 : fyerf=zero
1040 : 0 : fzerf=zero
1041 : : ! ion-electron term, gaussian part
1042 : 0 : fxgau=zero
1043 : 0 : fygau=zero
1044 : 0 : fzgau=zero
1045 : :
1046 : : ! Building array of coefficients of the derivative of the gaussian part
1047 : 0 : cprime(1)=2._dp*psppar(0,2,ityp)-psppar(0,1,ityp)
1048 : 0 : cprime(2)=4._dp*psppar(0,3,ityp)-psppar(0,2,ityp)
1049 : 0 : cprime(3)=6._dp*psppar(0,4,ityp)-psppar(0,3,ityp)
1050 : 0 : cprime(4)=-psppar(0,4,ityp)
1051 : :
1052 : : ! Determine number of local terms (HGH pot)
1053 : 0 : nloc=0
1054 [ # # ]: 0 : do iloc=1,4
1055 [ # # ]: 0 : if (psppar(0,iloc,ityp).ne.zero) nloc=iloc
1056 : : end do
1057 : :
1058 : : ! Some constants depending on the atom type
1059 : 0 : rloc=psppar(0,0,ityp) ; rloc2=rloc**2
1060 : 0 : charge=real(nelpsp(ityp),kind=dp)
1061 : 0 : prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1062 : :
1063 : : ! PAW specifics
1064 [ # # ]: 0 : if (usepaw==1) then
1065 : 0 : msz=pawtab(ityp)%wvl%rholoc%msz
1066 : 0 : rad => pawtab(ityp)%wvl%rholoc%rad(1:msz)
1067 : 0 : vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
1068 : 0 : d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
1069 [ # # ]: 0 : rzero=rad(1);if (rad(1)<=1.d-10) rzero=rad(2)
1070 : : end if
1071 : :
1072 : : ! Maximum extension of the gaussian
1073 : 0 : cutoff=10._dp*rloc
1074 [ # # ]: 0 : isx=floor((rx-cutoff)/hxh)
1075 [ # # ]: 0 : isy=floor((ry-cutoff)/hyh)
1076 [ # # ]: 0 : isz=floor((rz-cutoff)/hzh)
1077 [ # # ]: 0 : iex=ceiling((rx+cutoff)/hxh)
1078 [ # # ]: 0 : iey=ceiling((ry+cutoff)/hyh)
1079 [ # # ]: 0 : iez=ceiling((rz+cutoff)/hzh)
1080 : :
1081 : : ! Calculate the forces near the atom due to the gaussian
1082 : : ! and error function parts of the potential
1083 [ # # ]: 0 : do i3=isz,iez
1084 : 0 : z=real(i3,kind=dp)*hzh-rz
1085 : 0 : call ind_positions_mklocl(perz,i3,n3,j3,goz)
1086 [ # # ]: 0 : if(fftn3_distrib(j3)==iproc) then
1087 : 0 : i3loc=ffti3_local(j3)
1088 [ # # ]: 0 : do i2=isy,iey
1089 : 0 : y=real(i2,kind=dp)*hyh-ry
1090 : 0 : call ind_positions_mklocl(pery,i2,n2,j2,goy)
1091 [ # # ]: 0 : do i1=isx,iex
1092 : 0 : x=real(i1,kind=dp)*hxh-rx
1093 : 0 : call ind_positions_mklocl(perx,i1,n1,j1,gox)
1094 : :
1095 [ # # # # : 0 : if (goz.and.goy.and.gox) then
# # ]
1096 : 0 : ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2
1097 : 0 : x=(gridcart(1,ind)-rx)
1098 : 0 : y=(gridcart(2,ind)-ry)
1099 : 0 : z=(gridcart(3,ind)-rz)
1100 : 0 : r2=x**2+y**2+z**2
1101 : 0 : xp=exp(-0.5_dp*r2/rloc**2)
1102 : :
1103 : : ! Short range part
1104 : 0 : rhoel=rho(ind)
1105 : : ! HGH: V_S^prime=gaussian
1106 [ # # ]: 0 : if (usepaw==0) then
1107 [ # # ]: 0 : if (nloc/=0) then
1108 : 0 : arg=r2/rloc**2
1109 : 0 : tt=cprime(nloc)
1110 [ # # ]: 0 : do iloc=nloc-1,1,-1
1111 : 0 : tt=arg*tt+cprime(iloc)
1112 : : end do
1113 : 0 : forceloc=xp*tt*rhoel
1114 : : else
1115 : : forceloc=zero
1116 : : end if
1117 : : ! PAW: V_PAW^prime-V_L^prime
1118 : : else
1119 : 0 : rr(1)=sqrt(r2)
1120 [ # # ]: 0 : if (rr(1)>=rzero) then
1121 [ # # # # : 0 : call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
# # ]
1122 : 0 : call calcdVloc_mklocl(dvhgh,rr(1),rloc,charge)
1123 : 0 : forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1124 : : else
1125 [ # # # # : 0 : forceloc=rhoel*rloc2*dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
# # ]
1126 : : end if
1127 : : end if
1128 : :
1129 : 0 : fxgau=fxgau+forceloc*x
1130 : 0 : fygau=fygau+forceloc*y
1131 : 0 : fzgau=fzgau+forceloc*z
1132 : :
1133 : : ! Long range part: error function
1134 : 0 : Vel=pot(ind)
1135 : 0 : fxerf=fxerf+xp*Vel*x
1136 : 0 : fyerf=fyerf+xp*Vel*y
1137 : 0 : fzerf=fzerf+xp*Vel*z
1138 : :
1139 [ # # ]: 0 : else if (nloc>0) then
1140 : 0 : r2=x**2+y**2+z**2
1141 : 0 : arg=r2/rloc**2
1142 : 0 : xp=exp(-0.5_dp*arg)
1143 : 0 : tt=cprime(nloc)
1144 [ # # ]: 0 : do iloc=nloc-1,1,-1
1145 : 0 : tt=arg*tt+cprime(iloc)
1146 : : end do
1147 : 0 : forceleaked=forceleaked+xp*(1._dp+tt)
1148 : : end if
1149 : : end do
1150 : : end do
1151 : : end if
1152 : : end do
1153 : :
1154 : : ! Final result of the forces
1155 : 0 : floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc**2)*fxgau
1156 : 0 : floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc**2)*fygau
1157 : 0 : floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc**2)*fzgau
1158 : :
1159 : : end do
1160 : :
1161 : 0 : forceleaked=forceleaked*prefactor*hxh*hyh*hzh
1162 [ # # # # ]: 0 : if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1163 : :
1164 : : CONTAINS
1165 : : !!***
1166 : :
1167 : : !----------------------------------------------------------------------
1168 : :
1169 : : !!****f* mklocl_realspace/calcdVloc_mklocl
1170 : : !! NAME
1171 : : !! calcdVloc_mklocl
1172 : : !!
1173 : : !! FUNCTION
1174 : : !! Compute 1st-derivative of long-range HGH local ionic potential (derf)
1175 : : !!
1176 : : !! INPUTS
1177 : : !!
1178 : : !! OUTPUT
1179 : : !!
1180 : : !! SOURCE
1181 : :
1182 : 0 : subroutine calcdVloc_mklocl(yy,xx,rloc,Z)
1183 : :
1184 : : !Arguments ------------------------------------
1185 : : !scalars
1186 : : real(dp),intent(in) :: xx,rloc,Z
1187 : : real(dp),intent(out) :: yy
1188 : :
1189 : : !Local variables-------------------------------
1190 : : !scalars
1191 : : real(dp):: arg,tt
1192 : :
1193 : : ! *************************************************************************
1194 : :
1195 : 0 : arg=xx/(sqrt(2._dp)*rloc)
1196 : 0 : call derf_ab(tt,arg)
1197 : 0 : yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1198 : :
1199 : 0 : end subroutine calcdVloc_mklocl
1200 : : !!***
1201 : :
1202 : : !----------------------------------------------------------------------
1203 : :
1204 : : !!****f* mklocl_wavelets/dvloc_zero_mklocl
1205 : : !! NAME
1206 : : !! dvloc_zero_mklocl
1207 : : !!
1208 : : !! FUNCTION
1209 : : !! Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0
1210 : : !!
1211 : : !! INPUTS
1212 : : !!
1213 : : !! OUTPUT
1214 : : !!
1215 : : !! SOURCE
1216 : :
1217 : 0 : function dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
1218 : :
1219 : :
1220 : : !Arguments ------------------------------------
1221 : : !scalars
1222 : : integer,intent(in) :: msz
1223 : : real(dp) :: dvloc_zero_mklocl
1224 : : real(dp),intent(in) :: charge,rloc
1225 : : !arrays
1226 : : real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
1227 : :
1228 : : !Local variables-------------------------------
1229 : : !scalars
1230 : : real(dp) :: y1,y2,y3,zz=0._dp
1231 : : !arrays
1232 : : real(dp) :: ll(3),xx(3),yy(3)
1233 : :
1234 : : ! *************************************************************************
1235 : :
1236 : : !Select 3 points x1,x2,x3 near 0
1237 [ # # ]: 0 : if (rad(1)>1.d-10) then
1238 [ # # ]: 0 : xx(1:3)=rad(1:3)
1239 : : else
1240 [ # # ]: 0 : xx(1:3)=rad(2:4)
1241 : : end if
1242 : :
1243 : : !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
1244 : 0 : call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1245 : 0 : call calcdVloc_mklocl(y1,xx(1),rloc,charge)
1246 : 0 : call calcdVloc_mklocl(y2,xx(2),rloc,charge)
1247 : 0 : call calcdVloc_mklocl(y3,xx(3),rloc,charge)
1248 : 0 : yy(1)=(yy(1)-y1)/xx(1)
1249 : 0 : yy(2)=(yy(2)-y2)/xx(2)
1250 : 0 : yy(3)=(yy(3)-y3)/xx(3)
1251 : :
1252 : : !Find a polynomial of the form (z=0):
1253 : : !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1254 : :
1255 : : !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1256 : 0 : ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1257 : : !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1258 : 0 : ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1259 : : !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1260 : 0 : ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1261 : :
1262 : 0 : dvloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1263 : :
1264 : 0 : end function dvloc_zero_mklocl
1265 : : !!***
1266 : :
1267 : : end subroutine local_forces_new
1268 : : !!***
1269 : :
1270 : : !----------------------------------------------------------------------
1271 : :
1272 : : !!****f* mklocl_realspace/ind_positions_mklocl
1273 : : !! NAME
1274 : : !! ind_positions_mklocl
1275 : : !!
1276 : : !! FUNCTION
1277 : : !! determine the index in which the potential must be inserted, following the BC
1278 : : !! determine also whether the index is inside or outside the box for free BC
1279 : : !!
1280 : : !! INPUTS
1281 : : !!
1282 : : !! OUTPUT
1283 : : !!
1284 : : !! SOURCE
1285 : :
1286 : 0 : subroutine ind_positions_mklocl(periodic,i,n,j,go)
1287 : :
1288 : :
1289 : : !Arguments -------------------------------
1290 : : logical, intent(in) :: periodic
1291 : : integer, intent(in) :: i,n
1292 : : logical, intent(out) :: go
1293 : : integer, intent(out) :: j
1294 : :
1295 : : ! *********************************************************************
1296 : :
1297 [ # # ]: 0 : if (periodic) then
1298 : 0 : go=.true.
1299 : 0 : j=modulo(i-1,n)+1
1300 : : else
1301 : 0 : j=i
1302 [ # # # # ]: 0 : go=(i >= 1 .and. i <= n)
1303 : : end if
1304 : :
1305 : 0 : end subroutine ind_positions_mklocl
1306 : : !!***
1307 : :
1308 : : !!****f* ABINIT/mklocl_wavelets
1309 : : !!
1310 : : !! NAME
1311 : : !! mklocl_wavelets
1312 : : !!
1313 : : !! FUNCTION
1314 : : !! Compute the ionic local potential when the pseudo-potentials are GTH, using
1315 : : !! the special decomposition of these pseudo. The resulting potential is computed with
1316 : : !! free boundary conditions. It gives the same result than mklocl_realspace for the
1317 : : !! GTH pseudo only with a different way to compute the potential.
1318 : : !!
1319 : : !! Optionally compute :
1320 : : !! option=1 : local ionic potential throughout unit cell
1321 : : !! option=2 : contribution of local ionic potential to E gradient wrt xred
1322 : : !!
1323 : : !! INPUTS
1324 : : !! efield (3)=external electric field
1325 : : !! mpi_enreg=information about MPI parallelization
1326 : : !! natom=number of atoms
1327 : : !! nfft=size of vpsp (local potential)
1328 : : !! nspden=number of spin-density components
1329 : : !! option=type of calculation (potential, forces, ...)
1330 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1331 : : !! xcart(3,natom)=cartesian atomic coordinates.
1332 : : !! wvl_den=density-potential BigDFT object
1333 : : !! wvl_descr=wavelet BigDFT object
1334 : : !!
1335 : : !! OUTPUT
1336 : : !! (if option==1) vpsp(nfft)=the potential resulting from the ionic
1337 : : !! density of charge.
1338 : : !! (if option==2) grtn(3,natom)=grads of Etot wrt tn. These gradients are in
1339 : : !! reduced coordinates. Multiply them by rprimd to get
1340 : : !! gradients in cartesian coordinates.
1341 : : !!
1342 : : !! SOURCE
1343 : :
1344 : 0 : subroutine mklocl_wavelets(efield, grtn, mpi_enreg, natom, nfft, &
1345 : 0 : & nspden, option, rprimd, vpsp, wvl_den, wvl_descr, xcart)
1346 : :
1347 : : use defs_wvltypes
1348 : : #if defined HAVE_BIGDFT
1349 : : use BigDFT_API, only : ELECTRONIC_DENSITY,createIonicPotential,local_forces
1350 : : use poisson_solver, only : H_potential
1351 : : #endif
1352 : :
1353 : : !Arguments ------------------------------------
1354 : : !scalars
1355 : : integer,intent(in) :: option, natom, nfft, nspden
1356 : : type(MPI_type),intent(in) :: mpi_enreg
1357 : : type(wvl_denspot_type), intent(inout) :: wvl_den
1358 : : type(wvl_internal_type), intent(in) :: wvl_descr
1359 : : !arrays
1360 : : real(dp),intent(in) :: rprimd(3,3),efield(3)
1361 : : real(dp),intent(inout) :: grtn(3,natom)
1362 : : real(dp), intent(inout) :: vpsp(nfft)
1363 : : real(dp),intent(inout) :: xcart(3,natom)
1364 : :
1365 : : !Local variables-------------------------------
1366 : : #if defined HAVE_BIGDFT
1367 : : !scalars
1368 : : integer :: i,i1,i2,i3,ia,ierr,igeo,me,nproc,shift,comm
1369 : : real(dp) :: energ
1370 : : character(len=500) :: message
1371 : : !arrays
1372 : : real(dp) :: epot(3)
1373 : : real(dp) :: elecfield(3)=(/zero,zero,zero/) ! Not used here
1374 : : real(dp),allocatable :: gxyz(:,:),vhartr(:),rhov(:,:)
1375 : : #endif
1376 : :
1377 : : ! *********************************************************************
1378 : :
1379 : : #if defined HAVE_BIGDFT
1380 : :
1381 : : elecfield=zero !not used here
1382 : :
1383 : : !Manage parallelism
1384 : : comm=mpi_enreg%comm_wvl
1385 : : nproc=xmpi_comm_size(comm)
1386 : : me=xmpi_comm_rank(comm)
1387 : :
1388 : : !----------------------------------------------------------------------
1389 : : ! ----- Option 1: compute local ionic potential -----
1390 : : !----------------------------------------------------------------------
1391 : : if (option == 1) then
1392 : :
1393 : : ! We get the kernel for the Poisson solver (used to go from the ionic
1394 : : ! charge to the potential).If the kernel is uncomputed, it does it now.
1395 : : ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 2, icoulomb, me, kernel, &
1396 : : !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder)
1397 : : ! if (.not.associated(kernel%co%kernel)) then
1398 : : ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 1, icoulomb, me, kernel, &
1399 : : !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder)
1400 : : ! end if
1401 : :
1402 : : message=ch10//' mklocl_wavelets: Create local potential from ions.'
1403 : : call wrtout(std_out,message,'COLL')
1404 : :
1405 : : shift = 1 + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1406 : : & * wvl_den%denspot%dpbox%i3xcsh
1407 : :
1408 : : ! Call the BigDFT routine
1409 : : call createIonicPotential(wvl_descr%atoms%astruct%geocode, me, nproc, (me == 0), wvl_descr%atoms, &
1410 : : & xcart, wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), &
1411 : : & wvl_den%denspot%dpbox%hgrids(3), &
1412 : : & elecfield, wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, &
1413 : : & wvl_den%denspot%dpbox%n3pi, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, &
1414 : : & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), &
1415 : : & wvl_den%denspot%dpbox%ndims(3), wvl_den%denspot%pkernel, vpsp(shift), 0.d0,wvl_descr%rholoc)
1416 : :
1417 : : ! Copy vpsp into bigdft object:
1418 : : call wvl_rhov_abi2big(1,vpsp,wvl_den%denspot%v_ext,shift=(shift-1))
1419 : :
1420 : : ! Eventually add the electric field
1421 : : if (maxval(efield) > tol12) then
1422 : : message=ch10//'mklocl_wavelets: Add the electric field.'
1423 : : call wrtout(std_out,message,'COLL')
1424 : : ! We add here the electric field since in BigDFT, the field must be on x...
1425 : : epot(:) = real(0.5, dp) * efield(:) * wvl_den%denspot%dpbox%hgrids(:)
1426 : : do i3 = 1, wvl_den%denspot%dpbox%n3pi, 1
1427 : : ia = (i3 - 1) * wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2)
1428 : : do i2 = -14, 2 * wvl_descr%Glr%d%n2 + 16, 1
1429 : : i = ia + (i2 + 14) * wvl_den%denspot%dpbox%ndims(1)
1430 : : do i1 = -14, 2 * wvl_descr%Glr%d%n1 + 16, 1
1431 : : i = i + 1
1432 : : vpsp(shift + i) = vpsp(shift + i) + &
1433 : : & epot(1) * real(i1 - wvl_descr%Glr%d%n1, dp) + &
1434 : : & epot(2) * real(i2 - wvl_descr%Glr%d%n2, dp) + &
1435 : : & epot(3) * real(i3 - wvl_descr%Glr%d%n3, dp)
1436 : : end do
1437 : : end do
1438 : : end do
1439 : : end if
1440 : :
1441 : : !----------------------------------------------------------------------
1442 : : ! ----- Option 2: compute forces induced by local ionic potential -----
1443 : : !----------------------------------------------------------------------
1444 : : else if (option == 2) then
1445 : :
1446 : : message=ch10//' mklocl_wavelets: compute local forces.'
1447 : : call wrtout(std_out,message,'COLL')
1448 : :
1449 : : if (wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then
1450 : : message='denspot bigdft datstructure should contain rhor!'
1451 : : ABI_BUG(message)
1452 : : end if
1453 : :
1454 : : ! Extract density rhor from bigDFT datastructure
1455 : : ABI_MALLOC(rhov,(nfft, nspden))
1456 : : ABI_MALLOC(vhartr,(nfft))
1457 : : rhov=zero ; vhartr=zero
1458 : : shift = wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1459 : : & * wvl_den%denspot%dpbox%i3xcsh
1460 : : do i = 1, min(nfft,size(wvl_den%denspot%rhov)-shift)
1461 : : rhov(i, 1) = wvl_den%denspot%rhov(i + shift)
1462 : : vhartr(i) = wvl_den%denspot%rhov(i + shift)
1463 : : end do
1464 : : if (nspden == 2) then
1465 : : shift = shift + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1466 : : & * wvl_den%denspot%dpbox%n3d
1467 : : do i = 1, min(nfft,size(wvl_den%denspot%rhov)-shift)
1468 : : rhov(i, 2) = wvl_den%denspot%rhov(i + shift)
1469 : : vhartr(i) = vhartr(i) + wvl_den%denspot%rhov(i + shift)
1470 : : end do
1471 : : end if
1472 : :
1473 : : ! Compute Hartree potential from rhor
1474 : : call H_potential('D',wvl_den%denspot%pkernel,vhartr,vhartr,energ,zero,.false.)
1475 : :
1476 : : ! Allocate temporary array for forces
1477 : : ABI_MALLOC(gxyz,(3, natom))
1478 : :
1479 : : ! Calculate local part of the forces grtn (modified BigDFT routine)
1480 : : call local_forces_wvl(me,natom,xcart,&
1481 : : & wvl_den%denspot%dpbox%hgrids(1),&
1482 : : & wvl_den%denspot%dpbox%hgrids(2),&
1483 : : & wvl_den%denspot%dpbox%hgrids(3),&
1484 : : & wvl_descr%Glr%d%n1,wvl_descr%Glr%d%n2,wvl_descr%Glr%d%n3,&
1485 : : & wvl_den%denspot%dpbox%n3p,&
1486 : : & wvl_den%denspot%dpbox%i3s+wvl_den%denspot%dpbox%i3xcsh,&
1487 : : & wvl_den%denspot%dpbox%ndims(1),wvl_den%denspot%dpbox%ndims(2),&
1488 : : & rhov,vhartr,gxyz,wvl_descr)
1489 : : ! call local_forces(me, wvl_descr%atoms, xcart, &
1490 : : ! & wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), &
1491 : : ! & wvl_den%denspot%dpbox%hgrids(3), &
1492 : : ! & wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, &
1493 : : ! & wvl_den%denspot%dpbox%n3p, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, &
1494 : : ! & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), &
1495 : : ! & rhov, vhartr,gxyz,locstrten,charge)
1496 : :
1497 : : ! Pending: floc,locstrten and charge are not used here.
1498 : : ! Pending: put mpi_enreg%nscatterarr... in object denspot, initialize object, etc.
1499 : :
1500 : : if (nproc > 1) then
1501 : : call xmpi_sum(gxyz, comm, ierr)
1502 : : end if
1503 : :
1504 : : ! Forces should be in reduced coordinates.
1505 : : do ia = 1, natom, 1
1506 : : do igeo = 1, 3, 1
1507 : : grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) &
1508 : : & - rprimd(2, igeo) * gxyz(2, ia) &
1509 : : & - rprimd(3, igeo) * gxyz(3, ia)
1510 : : end do
1511 : : end do
1512 : :
1513 : : ! Deallocate local variables
1514 : : ABI_FREE(vhartr)
1515 : : ABI_FREE(rhov)
1516 : : ABI_FREE(gxyz)
1517 : :
1518 : : !----------------------------------------------------------------------
1519 : :
1520 : : else ! option switch
1521 : : message = 'Internal error, option should be 1 or 2!'
1522 : : ABI_ERROR(message)
1523 : : end if
1524 : :
1525 : : #else
1526 : 0 : BIGDFT_NOTENABLED_ERROR()
1527 : : if (.false.) write(std_out,*) option,natom,nfft,nspden,mpi_enreg%me,&
1528 : : & wvl_den%symObj,wvl_descr%h(1),rprimd(1,1),efield(1),grtn(1,1),vpsp(1),xcart(1,1)
1529 : : #endif
1530 : :
1531 : 0 : end subroutine mklocl_wavelets
1532 : : !!***
1533 : :
1534 : : !----------------------------------------------------------------------
1535 : :
1536 : : !!****f* mklocl_wavelets/local_forces_wvl
1537 : : !! NAME
1538 : : !! local_forces_wvl
1539 : : !!
1540 : : !! FUNCTION
1541 : : !!
1542 : : !! INPUTS
1543 : : !! hxh,hyh,hzh=wavelet grid spacings
1544 : : !! iproc=current MPI process number
1545 : : !! n1,n2,n3=number of wavelet points in each direction
1546 : : !! n1i,n2i=size of intermediate xy wvl grid
1547 : : !! n3pi=number of xy wvl planes handled by current MPI process
1548 : : !! i3s=starting index of local potential for current MPI process
1549 : : !! natom=number of atoms
1550 : : !! pot(*)=Hartree ionic potential
1551 : : !! rho(*)=electronic density
1552 : : !! rxyz(3,natom)=cartesian coordinates of atoms
1553 : : !! wvl=wavelet BigDFT object
1554 : : !!
1555 : : !! OUTPUT
1556 : : !! floc(3,natom)=local ionic potential contribution to forces
1557 : : !!
1558 : : !! SOURCE
1559 : :
1560 : : subroutine local_forces_wvl(iproc,natom,rxyz,hxh,hyh,hzh,n1,n2,n3,n3pi,i3s,n1i,n2i,&
1561 : : & rho,pot,floc,wvl)
1562 : :
1563 : : use defs_wvltypes
1564 : : #if defined HAVE_BIGDFT
1565 : : use BigDFT_API, only : PSPCODE_PAW,ind_positions
1566 : : #endif
1567 : :
1568 : : !Arguments -------------------------------
1569 : : !scalars
1570 : : integer,intent(in) :: i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom
1571 : : real(dp),intent(in) :: hxh,hyh,hzh
1572 : : type(wvl_internal_type),intent(in) :: wvl
1573 : : !arrays
1574 : : real(dp),intent(in) :: rxyz(3,natom)
1575 : : real(dp),dimension(*),intent(in) :: rho,pot
1576 : : real(dp),intent(out) :: floc(3,natom)
1577 : :
1578 : : !Local variables -------------------------
1579 : : #if defined HAVE_BIGDFT
1580 : : !scalars
1581 : : integer :: i1,i2,i3,iat,iex,iey,iez,iloc,ind,isx,isy,isz,ityp
1582 : : integer :: j1,j2,j3,msz=0
1583 : : integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,nloc
1584 : : logical :: perx,pery,perz,gox,goy,goz,USE_PAW
1585 : : real(dp) :: arg,charge,cutoff,dvhgh,forceloc,fxerf,fyerf,fzerf,fxgau,fygau,fzgau
1586 : : real(dp) :: forceleaked,prefactor,r2,rhoel,rloc,rloc2,rx,ry,rz,rzero,tt,vel,x,xp,y,z
1587 : : !arrays
1588 : : real(dp) :: cprime(4),dvpawdr(1),rr(1)
1589 : : real(dp),pointer :: psppar(:,:),rad(:),vloc(:),d2vloc(:)
1590 : : #endif
1591 : :
1592 : : ! *********************************************************************
1593 : :
1594 : : #if defined HAVE_BIGDFT
1595 : :
1596 : : if (iproc==0) write(std_out,'(a)',advance='no') ' Calculate local forces...'
1597 : :
1598 : : !PAW or NCPP ?
1599 : : USE_PAW=any(wvl%atoms%npspcode==PSPCODE_PAW)
1600 : :
1601 : : !Conditions for periodicity in the three directions
1602 : : perx=(wvl%atoms%astruct%geocode /= 'F')
1603 : : pery=(wvl%atoms%astruct%geocode == 'P')
1604 : : perz=(wvl%atoms%astruct%geocode /= 'F')
1605 : : call ext_buffers(perx,nbl1,nbr1)
1606 : : call ext_buffers(pery,nbl2,nbr2)
1607 : : call ext_buffers(perz,nbl3,nbr3)
1608 : :
1609 : : forceleaked=zero
1610 : :
1611 : : do iat=1,natom
1612 : : ityp=wvl%atoms%astruct%iatype(iat)
1613 : :
1614 : : ! Coordinates of the center
1615 : : rx=rxyz(1,iat)
1616 : : ry=rxyz(2,iat)
1617 : : rz=rxyz(3,iat)
1618 : :
1619 : : ! Initialization of the forces
1620 : : ! ion-electron term, error function part
1621 : : fxerf=zero
1622 : : fyerf=zero
1623 : : fzerf=zero
1624 : : ! ion-electron term, gaussian part
1625 : : fxgau=zero
1626 : : fygau=zero
1627 : : fzgau=zero
1628 : :
1629 : : ! Building array of coefficients of the derivative of the gaussian part
1630 : : psppar => wvl%atoms%psppar(:,:,ityp)
1631 : : cprime(1)=2._dp*wvl%atoms%psppar(0,2,ityp)-wvl%atoms%psppar(0,1,ityp)
1632 : : cprime(2)=4._dp*wvl%atoms%psppar(0,3,ityp)-wvl%atoms%psppar(0,2,ityp)
1633 : : cprime(3)=6._dp*wvl%atoms%psppar(0,4,ityp)-wvl%atoms%psppar(0,3,ityp)
1634 : : cprime(4)=-wvl%atoms%psppar(0,4,ityp)
1635 : :
1636 : : ! Determine number of local terms (HGH pot)
1637 : : nloc=0
1638 : : do iloc=1,4
1639 : : if (wvl%atoms%psppar(0,iloc,ityp).ne.zero) nloc=iloc
1640 : : end do
1641 : :
1642 : : ! Some constants depending on the atom type
1643 : : rloc=wvl%atoms%psppar(0,0,ityp) ; rloc2=rloc**2
1644 : : charge=real(wvl%atoms%nelpsp(ityp),kind=dp)
1645 : : prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1646 : :
1647 : : ! PAW specifics
1648 : : if (USE_PAW) then
1649 : : msz=wvl%rholoc%msz(ityp)
1650 : : rad => wvl%rholoc%rad(1:msz,ityp)
1651 : : vloc => wvl%rholoc%d(1:msz,3,ityp)
1652 : : d2vloc => wvl%rholoc%d(1:msz,4,ityp)
1653 : : rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2)
1654 : : end if
1655 : :
1656 : : ! Maximum extension of the gaussian
1657 : : cutoff=10._dp*rloc
1658 : : isx=floor((rx-cutoff)/hxh)
1659 : : isy=floor((ry-cutoff)/hyh)
1660 : : isz=floor((rz-cutoff)/hzh)
1661 : : iex=ceiling((rx+cutoff)/hxh)
1662 : : iey=ceiling((ry+cutoff)/hyh)
1663 : : iez=ceiling((rz+cutoff)/hzh)
1664 : :
1665 : : ! Calculate the forces near the atom due to the gaussian
1666 : : ! and error function parts of the potential
1667 : : if (n3pi>0) then
1668 : : do i3=isz,iez
1669 : : z=real(i3,kind=dp)*hzh-rz
1670 : : call ind_positions(perz,i3,n3,j3,goz)
1671 : : j3=j3+nbl3+1
1672 : : do i2=isy,iey
1673 : : y=real(i2,kind=dp)*hyh-ry
1674 : : call ind_positions(pery,i2,n2,j2,goy)
1675 : : do i1=isx,iex
1676 : : x=real(i1,kind=dp)*hxh-rx
1677 : : call ind_positions(perx,i1,n1,j1,gox)
1678 : :
1679 : : r2=x**2+y**2+z**2
1680 : : xp=exp(-0.5_dp*r2/rloc2)
1681 : :
1682 : : if ((j3>=i3s.and.j3<=i3s+n3pi-1) .and. (gox.and.goy)) then
1683 : : ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s+1-1)*n1i*n2i
1684 : :
1685 : : ! Short range part
1686 : : rhoel=rho(ind)
1687 : : ! HGH: V_S^prime=gaussian
1688 : : if (.not.USE_PAW) then
1689 : : if (nloc/=0) then
1690 : : arg=r2/rloc2
1691 : : tt=cprime(nloc)
1692 : : do iloc=nloc-1,1,-1
1693 : : tt=arg*tt+cprime(iloc)
1694 : : end do
1695 : : forceloc=xp*tt*rhoel
1696 : : else
1697 : : forceloc=zero
1698 : : end if
1699 : : ! PAW: V_PAW^prime-V_L^prime
1700 : : else
1701 : : rr(1)=sqrt(r2)
1702 : :
1703 : : if (rr(1)>=rzero) then
1704 : : call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
1705 : : call calcdVloc_wvl(dvhgh,rr(1),rloc,charge)
1706 : : forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1707 : : else
1708 : : forceloc=rhoel*rloc2*dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc)
1709 : : end if
1710 : : end if
1711 : :
1712 : : fxgau=fxgau+forceloc*x
1713 : : fygau=fygau+forceloc*y
1714 : : fzgau=fzgau+forceloc*z
1715 : :
1716 : : ! Long range part: error function
1717 : : vel=pot(ind)
1718 : : fxerf=fxerf+xp*vel*x
1719 : : fyerf=fyerf+xp*vel*y
1720 : : fzerf=fzerf+xp*vel*z
1721 : :
1722 : : else if ((.not.goz).and.(nloc>0)) then
1723 : : arg=r2/rloc2
1724 : : tt=cprime(nloc)
1725 : : do iloc=nloc-1,1,-1
1726 : : tt=arg*tt+cprime(iloc)
1727 : : end do
1728 : : forceleaked=forceleaked+prefactor*xp*tt*rho(1)
1729 : : end if
1730 : :
1731 : : end do ! i1
1732 : : end do ! i2
1733 : : end do ! i3
1734 : : end if ! n3pi>0
1735 : :
1736 : : ! Final result of the forces
1737 : : floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc2)*fxgau
1738 : : floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc2)*fygau
1739 : : floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc2)*fzgau
1740 : :
1741 : : end do ! iat
1742 : :
1743 : : forceleaked=forceleaked*hxh*hyh*hzh
1744 : : if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1745 : :
1746 : : #else
1747 : : BIGDFT_NOTENABLED_ERROR()
1748 : : if (.false.) write(std_out,*) i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom,hxh,hyh,hzh,&
1749 : : & rxyz(1,1),floc(1,1),rho(1),pot(1),wvl%h(1)
1750 : : #endif
1751 : :
1752 : : CONTAINS
1753 : : !!***
1754 : :
1755 : : !----------------------------------------------------------------------
1756 : :
1757 : : !!****f* mklocl_wavelets/calcdVloc_wvl
1758 : : !! NAME
1759 : : !! calcdVloc_wvl
1760 : : !!
1761 : : !! FUNCTION
1762 : : !! Compute 1st-derivative of long-range HGH local ionic potential (derf)
1763 : : !!
1764 : : !! INPUTS
1765 : : !!
1766 : : !! OUTPUT
1767 : : !!
1768 : : !! SOURCE
1769 : :
1770 : : subroutine calcdVloc_wvl(yy,xx,rloc,Z)
1771 : :
1772 : : !Arguments ------------------------------------
1773 : : !scalars
1774 : : real(dp),intent(in) :: xx,rloc,Z
1775 : : real(dp),intent(out) :: yy
1776 : :
1777 : : !Local variables-------------------------------
1778 : : !scalars
1779 : : real(dp):: arg,tt
1780 : :
1781 : : ! *************************************************************************
1782 : :
1783 : : arg=xx/(sqrt(2._dp)*rloc)
1784 : : call derf_ab(tt,arg)
1785 : : yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1786 : :
1787 : : end subroutine calcdVloc_wvl
1788 : : !!***
1789 : :
1790 : : !----------------------------------------------------------------------
1791 : :
1792 : : !!****f* mklocl_wavelets/dvloc_zero_wvl
1793 : : !! NAME
1794 : : !! dvloc_zero_wvl
1795 : : !!
1796 : : !! FUNCTION
1797 : : !! Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0
1798 : : !!
1799 : : !! INPUTS
1800 : : !!
1801 : : !! OUTPUT
1802 : : !!
1803 : : !! SOURCE
1804 : :
1805 : : function dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc)
1806 : :
1807 : :
1808 : : !Arguments ------------------------------------
1809 : : !scalars
1810 : : integer,intent(in) :: msz
1811 : : real(dp) :: dvloc_zero_wvl
1812 : : real(dp),intent(in) :: charge,rloc
1813 : : !arrays
1814 : : real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
1815 : :
1816 : : !Local variables-------------------------------
1817 : : !scalars
1818 : : real(dp) :: y1,y2,y3,zz=0._dp
1819 : : !arrays
1820 : : real(dp) :: ll(3),xx(3),yy(3)
1821 : :
1822 : : ! *************************************************************************
1823 : :
1824 : : !Select 3 points x1,x2,x3 near 0
1825 : : if (rad(1)>1.d-10) then
1826 : : xx(1:3)=rad(1:3)
1827 : : else
1828 : : xx(1:3)=rad(2:4)
1829 : : end if
1830 : :
1831 : : !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
1832 : : call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1833 : : call calcdVloc_wvl(y1,xx(1),rloc,charge)
1834 : : call calcdVloc_wvl(y2,xx(2),rloc,charge)
1835 : : call calcdVloc_wvl(y3,xx(3),rloc,charge)
1836 : : yy(1)=(yy(1)-y1)/xx(1)
1837 : : yy(2)=(yy(2)-y2)/xx(2)
1838 : : yy(3)=(yy(3)-y3)/xx(3)
1839 : :
1840 : : !Find a polynomial of the form (z=0):
1841 : : !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1842 : :
1843 : : !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1844 : : ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1845 : : !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1846 : : ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1847 : : !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1848 : : ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1849 : :
1850 : : dvloc_zero_wvl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1851 : :
1852 : : end function dvloc_zero_wvl
1853 : : !!***
1854 : :
1855 : : end subroutine local_forces_wvl
1856 : : !!***
1857 : :
1858 : : end module m_mklocl_realspace
1859 : : !!***
|
