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, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
142 0 : integer, ABI_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 :
614 : !Arguments -------------------------------
615 : !scalars
616 : integer, intent(in) :: iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw
617 : real(dp), intent(in) :: hxh,hyh,hzh
618 : character(len=1), intent(in) :: geocode
619 : type(coulomb_operator), intent(in) :: kernel
620 : !arrays
621 : integer, dimension(nat), intent(in) :: iatype
622 : integer, dimension(ntypes), intent(in) :: nelpsp
623 : integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local
624 : real(dp), dimension(3,n1i*n2i*n3d), intent(in) :: gridcart
625 : real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar
626 : real(dp), dimension(3,nat), intent(in) :: rxyz
627 : real(dp), dimension(*), intent(inout) :: pot_ion
628 : type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw)
629 :
630 : !Local variables -------------------------
631 : #if defined HAVE_BIGDFT
632 : !scalars
633 : integer :: iat,i1,i2,i3,j1,j2,j3,isx,isy,isz,iex,iey,iez,ierr,ityp
634 : integer :: ind,nloc,iloc,i3loc,msz
635 : logical :: gox,goy,goz,perx,pery,perz
636 : real(dp) :: arg,charge,cutoff,ehart,eexcu,rholeaked,rholeaked_tot,rloc
637 : real(dp) :: rx,ry,rz,rzero,r2,tt,tt_tot,vexcu,vhgh,x,xp,y,z
638 : !arrays
639 : real(dp) :: rr(1),vpaw(1)
640 : real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
641 : #endif
642 :
643 : ! *********************************************************************
644 :
645 : #if defined HAVE_BIGDFT
646 :
647 : if(nproc<0)then
648 : ABI_ERROR('nproc should not be negative')
649 : end if
650 :
651 : !Ionic charge (must be calculated for the PS active processes)
652 : rholeaked=0._dp
653 : !Ionic energy (can be calculated for all the processors)
654 :
655 : !here we should insert the calculation of the ewald energy for the periodic BC case
656 : !!! eion=0._dp
657 : !!! do iat=1,nat
658 : !!! ityp=iatype(iat)
659 : !!! rx=rxyz(1,iat)
660 : !!! ry=rxyz(2,iat)
661 : !!! rz=rxyz(3,iat)
662 : !!! ! ion-ion interaction
663 : !!! do jat=1,iat-1
664 : !!! dist=sqrt( (rx-rxyz(1,jat))**2+(ry-rxyz(2,jat))**2+(rz-rxyz(3,jat))**2 )
665 : !!! jtyp=iatype(jat)
666 : !!! eion=eion+real(nelpsp(jtyp)*nelpsp(ityp),kind=dp)/dist
667 : !!! enddo
668 : !!! end do
669 : !!! if (iproc.eq.0) write(std_out,'(1x,a,1pe22.14)') 'ion-ion interaction energy',eion
670 :
671 : !Creates charge density arising from the ionic PSP cores
672 : !the n3pi dimension indicates the number of planes trated by each processor in the FFT parallelisation
673 : !for a plane wave treatment this value depends on whether the direct space is divided in planes or not
674 : !I don't know this variable, which in the future must be inserted at the place of n3pi (LG)
675 : !if n3pi=0 this means that the processors doesn't calculate anything
676 : !if (n3pi >0 ) then
677 :
678 : !conditions for periodicity in the three directions
679 : perx=(geocode /= 'F')
680 : pery=(geocode == 'P')
681 : perz=(geocode /= 'F')
682 :
683 : !this initialise the array to zero, it will work only if bigdft library is enabled
684 : pot_ion(1:n1i*n2i*n3d)=zero
685 :
686 : do iat=1,nat
687 : ityp=iatype(iat)
688 : rx=rxyz(1,iat)
689 : ry=rxyz(2,iat)
690 : rz=rxyz(3,iat)
691 :
692 : rloc=psppar(0,0,ityp)
693 : charge=real(nelpsp(ityp),kind=dp)/(2._dp*pi*sqrt(2._dp*pi)*rloc**3)
694 : cutoff=10._dp*rloc
695 :
696 : isx=floor((rx-cutoff)/hxh)
697 : isy=floor((ry-cutoff)/hyh)
698 : isz=floor((rz-cutoff)/hzh)
699 :
700 : iex=ceiling((rx+cutoff)/hxh)
701 : iey=ceiling((ry+cutoff)/hyh)
702 : iez=ceiling((rz+cutoff)/hzh)
703 :
704 : ! Calculate Ionic Density
705 : ! using HGH parameters.
706 : ! Eq. 1.104, T. Deutsch and L. Genovese, JDN. 12, 2011
707 : do i3=isz,iez
708 : z=real(i3,kind=dp)*hzh-rz
709 : call ind_positions_mklocl(perz,i3,n3i,j3,goz)
710 : if(fftn3_distrib(j3)==iproc) then
711 : i3loc=ffti3_local(j3)
712 : do i2=isy,iey
713 : y=real(i2,kind=dp)*hyh-ry
714 : call ind_positions_mklocl(pery,i2,n2i,j2,goy)
715 : do i1=isx,iex
716 : x=real(i1,kind=dp)*hxh-rx
717 : call ind_positions_mklocl(perx,i1,n1i,j1,gox)
718 : r2=x**2+y**2+z**2
719 : if (goz .and. goy .and. gox ) then
720 : ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i
721 : r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2
722 : end if
723 : arg=r2/rloc**2
724 : xp=exp(-.5d0*arg)
725 : if (goz .and. goy .and. gox ) then
726 : pot_ion(ind)=pot_ion(ind)-xp*charge
727 : else
728 : rholeaked=rholeaked+xp*charge
729 : end if
730 : end do
731 : end do
732 : end if
733 : end do
734 :
735 : end do
736 :
737 : !Check
738 : tt=0._dp
739 : do j3= 1,n3d
740 : do i2= 1,n2i
741 : do i1= 1,n1i
742 : ind=i1+(i2-1)*n1i+(j3-1)*n1i*n2i
743 : tt=tt+pot_ion(ind)
744 : end do
745 : end do
746 : end do
747 :
748 : tt=tt*hxh*hyh*hzh
749 : rholeaked=rholeaked*hxh*hyh*hzh
750 :
751 : call xmpi_sum(tt,tt_tot,spaceworld,ierr)
752 : call xmpi_sum(rholeaked,rholeaked_tot,spaceworld,ierr)
753 :
754 : if (iproc.eq.0) then
755 : write(std_out,'(1x,a,f26.12,2x,1pe10.3)') &
756 : & 'total ionic charge, leaked charge ',tt_tot,rholeaked_tot
757 : end if
758 :
759 : !Here the value of the datacode must be kept fixed
760 : !there can be some problems when running this stuff in parallel,
761 : ! if the ionic potential distribution does not agree with the
762 : ! plane distribution which is supposed to hold for the Poisson Solver
763 : call psolver(geocode,'D',iproc,nproc,n1i,n2i,n3i,0,hxh,hyh,hzh,&
764 : & pot_ion,kernel%kernel,pot_ion,ehart,eexcu,vexcu,0._dp,.false.,1)
765 :
766 : !Add the remaining short-range local terms
767 : do iat=1,nat
768 : ityp=iatype(iat)
769 :
770 : rx=rxyz(1,iat)
771 : ry=rxyz(2,iat)
772 : rz=rxyz(3,iat)
773 :
774 : ! determine number of local terms
775 : rloc=psppar(0,0,ityp)
776 : cutoff=10._dp*rloc
777 : charge=real(nelpsp(ityp),kind=dp)
778 :
779 : ! determine number of local terms (HGH pot)
780 : nloc=0
781 : do iloc=1,4
782 : if (psppar(0,iloc,ityp).ne.0._dp) nloc=iloc
783 : end do
784 :
785 : ! PAW specifics
786 : if (usepaw==1) then
787 : msz=pawtab(ityp)%wvl%rholoc%msz
788 : rad => pawtab(ityp)%wvl%rholoc%rad(1:msz)
789 : vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
790 : d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
791 : rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2)
792 : end if
793 :
794 : isx=floor((rx-cutoff)/hxh)
795 : isy=floor((ry-cutoff)/hyh)
796 : isz=floor((rz-cutoff)/hzh)
797 :
798 : iex=ceiling((rx+cutoff)/hxh)
799 : iey=ceiling((ry+cutoff)/hyh)
800 : iez=ceiling((rz+cutoff)/hzh)
801 :
802 : do i3=isz,iez
803 : z=real(i3,kind=dp)*hzh-rz
804 : call ind_positions_mklocl(perz,i3,n3i,j3,goz)
805 : if(fftn3_distrib(j3) == iproc .and. goz) then !MPI
806 : i3loc=ffti3_local(j3)
807 : if (goz) then
808 : do i2=isy,iey
809 : y=real(i2,kind=dp)*hyh-ry
810 : call ind_positions_mklocl(pery,i2,n2i,j2,goy)
811 : if (goy) then
812 : do i1=isx,iex
813 : x=real(i1,kind=dp)*hxh-rx
814 : call ind_positions_mklocl(perx,i1,n1i,j1,gox)
815 : if (gox) then
816 : ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i
817 : r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2
818 : ! r2=x**2+y**2+z**2
819 :
820 : ! HGH: V_S=gaussian potential of Eq. (9) in JCP 129, 014109(2008)
821 : if (usepaw==0) then
822 : if (nloc /= 0) then
823 : arg=r2/rloc**2
824 : xp=exp(-.5d0*arg)
825 : tt=psppar(0,nloc,ityp)
826 : do iloc=nloc-1,1,-1
827 : tt=arg*tt+psppar(0,iloc,ityp)
828 : end do
829 : pot_ion(ind)=pot_ion(ind)+xp*tt
830 : end if
831 :
832 : ! PAW: V_PAW-V_L^HGH
833 : else
834 : rr(1)=sqrt(r2)
835 : if (rr(1)>=rzero) then
836 : call paw_splint(msz,rad,vloc,d2vloc,1,rr,vpaw)
837 : call calcVloc_mklocl(vhgh,rr(1),rloc,charge)
838 : pot_ion(ind)=pot_ion(ind)+vpaw(1)-vhgh
839 : else
840 : pot_ion(ind)=pot_ion(ind)+vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
841 : end if
842 : end if
843 :
844 : end if
845 : end do
846 : end if
847 : end do
848 : end if
849 : end if
850 : end do
851 :
852 : end do !iat
853 : #else
854 0 : BIGDFT_NOTENABLED_ERROR()
855 : if (.false.) write(std_out,*) geocode,iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw,&
856 : & hxh,hyh,hzh,iatype(1),nelpsp(1),fftn3_distrib(1),ffti3_local(1),gridcart(1,1),psppar(1,1,1),&
857 : & rxyz(1,1),pot_ion(1),pawtab(1)%mesh_size,kernel%co
858 : #endif
859 :
860 : CONTAINS
861 : !!***
862 :
863 : !----------------------------------------------------------------------
864 :
865 : !!****f* mklocl_realspace/calcVloc_mklocl
866 : !! NAME
867 : !! calcVloc_mklocl
868 : !!
869 : !! FUNCTION
870 : !!
871 : !! INPUTS
872 : !!
873 : !! OUTPUT
874 : !!
875 : !! SOURCE
876 :
877 : subroutine calcVloc_mklocl(yy,xx,rloc,Z)
878 :
879 : !Arguments ------------------------------------
880 : !scalars
881 : real(dp),intent(in) :: xx,rloc,Z
882 : real(dp),intent(out) :: yy
883 :
884 : !Local variables-------------------------------
885 : !scalars
886 : real(dp):: arg,tt
887 :
888 : ! *************************************************************************
889 :
890 : arg=xx/(sqrt(2.0)*rloc)
891 : call derf_ab(tt,arg)
892 : yy=-Z/xx*tt
893 :
894 : end subroutine calcVloc_mklocl
895 : !!***
896 :
897 : !----------------------------------------------------------------------
898 :
899 : !!****f* mklocl_realspace/vloc_zero_mklocl
900 : !! NAME
901 : !! vloc_zero_mklocl
902 : !!
903 : !! FUNCTION
904 : !! Use a quadratic interpolation to get limit of Vloc(x) at x->0
905 : !!
906 : !! INPUTS
907 : !!
908 : !! OUTPUT
909 : !!
910 : !! SOURCE
911 :
912 : function vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
913 :
914 :
915 : !Arguments ------------------------------------
916 : !scalars
917 : integer,intent(in) :: msz
918 : real(dp) :: vloc_zero_mklocl
919 : real(dp),intent(in) :: charge,rloc
920 : !arrays
921 : real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
922 :
923 : !Local variables-------------------------------
924 : !scalars
925 : real(dp) :: y1,y2,y3,zz=0._dp
926 : !arrays
927 : real(dp) :: ll(3),xx(3),yy(3)
928 :
929 : ! *************************************************************************
930 :
931 : !Select 3 points x1,x2,x3 near 0
932 : if (rad(1)>1.d-10) then
933 : xx(1:3)=rad(1:3)
934 : else
935 : xx(1:3)=rad(2:4)
936 : end if
937 :
938 : !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
939 : call paw_splint(msz,rad,vloc,d2vloc,3,xx,yy)
940 : call calcVloc_mklocl(y1,xx(1),rloc,charge)
941 : call calcVloc_mklocl(y2,xx(2),rloc,charge)
942 : call calcVloc_mklocl(y3,xx(3),rloc,charge)
943 : yy(1)= yy(1)-y1
944 : yy(2)= yy(2)-y2
945 : yy(3)= yy(3)-y3
946 :
947 : !Find a polynomial of the form (z=0):
948 : !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
949 :
950 : !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
951 : ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
952 : !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
953 : ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
954 : !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
955 : ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
956 :
957 : vloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
958 :
959 : end function vloc_zero_mklocl
960 : !!***
961 :
962 : end subroutine createIonicPotential_new
963 : !!***
964 :
965 : !----------------------------------------------------------------------
966 :
967 : !!****f* mklocl_realspace/local_forces_new
968 : !! NAME
969 : !! local_forces_new
970 : !!
971 : !! FUNCTION
972 : !!
973 : !! INPUTS
974 : !!
975 : !! OUTPUT
976 : !!
977 : !! SOURCE
978 :
979 0 : subroutine local_forces_new(fftn3_distrib,ffti3_local,&
980 0 : geocode,iproc,ntypes,nat,iatype,rxyz,gridcart,psppar,nelpsp,hxh,hyh,hzh,&
981 0 : n1,n2,n3,n3d,rho,pot,floc,pawtab,usepaw)
982 :
983 :
984 : !Arguments -------------------------------
985 : !scalars
986 : integer, intent(in) :: iproc,ntypes,nat,n1,n2,n3,n3d,usepaw
987 : character(len=1), intent(in) :: geocode
988 : real(dp), intent(in) :: hxh,hyh,hzh
989 : !arrays
990 : integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local
991 : integer, dimension(nat), intent(in) :: iatype
992 : integer, dimension(ntypes), intent(in) :: nelpsp
993 : real(dp), dimension(3,n1*n2*n3d), intent(in) :: gridcart
994 : real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar
995 : real(dp), dimension(3,nat), intent(in) :: rxyz
996 : real(dp), dimension(*), intent(in) :: rho,pot
997 : real(dp), dimension(3,nat), intent(out) :: floc
998 : type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw)
999 :
1000 : !Local variables -------------------------
1001 : !scalars
1002 : integer :: isx,isy,isz,iex,iey,iez,i1,i2,i3,j1,j2,j3,ind,iat,ityp,iloc,i3loc,msz,nloc
1003 : logical :: perx,pery,perz,gox,goy,goz
1004 : real(dp) :: arg,charge,cutoff,dvhgh,fxerf,fyerf,fzerf,fxgau,fygau,fzgau,forceleaked
1005 : real(dp) :: forceloc,prefactor,rloc,rloc2,rhoel,rx,ry,rz,rzero,r2,tt,x,xp,y,z,Vel
1006 : real(dp), dimension(4) :: cprime
1007 : !arrays
1008 : real(dp) :: dvpawdr(1),rr(1)
1009 0 : real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
1010 :
1011 : ! *********************************************************************
1012 :
1013 0 : if (iproc == 0) write(std_out,'(1x,a)',advance='no')'Calculate local forces...'
1014 :
1015 : !Conditions for periodicity in the three directions
1016 0 : perx=(geocode /= 'F')
1017 0 : pery=(geocode == 'P')
1018 0 : perz=(geocode /= 'F')
1019 :
1020 0 : forceleaked=zero
1021 :
1022 0 : do iat=1,nat
1023 0 : ityp=iatype(iat)
1024 : ! Coordinates of the center
1025 0 : rx=rxyz(1,iat)
1026 0 : ry=rxyz(2,iat)
1027 0 : rz=rxyz(3,iat)
1028 :
1029 : ! Initialization of the forces
1030 : ! ion-electron term, error function part
1031 0 : fxerf=zero
1032 0 : fyerf=zero
1033 0 : fzerf=zero
1034 : ! ion-electron term, gaussian part
1035 0 : fxgau=zero
1036 0 : fygau=zero
1037 0 : fzgau=zero
1038 :
1039 : ! Building array of coefficients of the derivative of the gaussian part
1040 0 : cprime(1)=2._dp*psppar(0,2,ityp)-psppar(0,1,ityp)
1041 0 : cprime(2)=4._dp*psppar(0,3,ityp)-psppar(0,2,ityp)
1042 0 : cprime(3)=6._dp*psppar(0,4,ityp)-psppar(0,3,ityp)
1043 0 : cprime(4)=-psppar(0,4,ityp)
1044 :
1045 : ! Determine number of local terms (HGH pot)
1046 0 : nloc=0
1047 0 : do iloc=1,4
1048 0 : if (psppar(0,iloc,ityp).ne.zero) nloc=iloc
1049 : end do
1050 :
1051 : ! Some constants depending on the atom type
1052 0 : rloc=psppar(0,0,ityp) ; rloc2=rloc**2
1053 0 : charge=real(nelpsp(ityp),kind=dp)
1054 0 : prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1055 :
1056 : ! PAW specifics
1057 0 : if (usepaw==1) then
1058 0 : msz=pawtab(ityp)%wvl%rholoc%msz
1059 0 : rad => pawtab(ityp)%wvl%rholoc%rad(1:msz)
1060 0 : vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
1061 0 : d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
1062 0 : rzero=rad(1);if (rad(1)<=1.d-10) rzero=rad(2)
1063 : end if
1064 :
1065 : ! Maximum extension of the gaussian
1066 0 : cutoff=10._dp*rloc
1067 0 : isx=floor((rx-cutoff)/hxh)
1068 0 : isy=floor((ry-cutoff)/hyh)
1069 0 : isz=floor((rz-cutoff)/hzh)
1070 0 : iex=ceiling((rx+cutoff)/hxh)
1071 0 : iey=ceiling((ry+cutoff)/hyh)
1072 0 : iez=ceiling((rz+cutoff)/hzh)
1073 :
1074 : ! Calculate the forces near the atom due to the gaussian
1075 : ! and error function parts of the potential
1076 0 : do i3=isz,iez
1077 0 : z=real(i3,kind=dp)*hzh-rz
1078 0 : call ind_positions_mklocl(perz,i3,n3,j3,goz)
1079 0 : if(fftn3_distrib(j3)==iproc) then
1080 0 : i3loc=ffti3_local(j3)
1081 0 : do i2=isy,iey
1082 0 : y=real(i2,kind=dp)*hyh-ry
1083 0 : call ind_positions_mklocl(pery,i2,n2,j2,goy)
1084 0 : do i1=isx,iex
1085 0 : x=real(i1,kind=dp)*hxh-rx
1086 0 : call ind_positions_mklocl(perx,i1,n1,j1,gox)
1087 :
1088 0 : if (goz.and.goy.and.gox) then
1089 0 : ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2
1090 0 : x=(gridcart(1,ind)-rx)
1091 0 : y=(gridcart(2,ind)-ry)
1092 0 : z=(gridcart(3,ind)-rz)
1093 0 : r2=x**2+y**2+z**2
1094 0 : xp=exp(-0.5_dp*r2/rloc**2)
1095 :
1096 : ! Short range part
1097 0 : rhoel=rho(ind)
1098 : ! HGH: V_S^prime=gaussian
1099 0 : if (usepaw==0) then
1100 0 : if (nloc/=0) then
1101 0 : arg=r2/rloc**2
1102 0 : tt=cprime(nloc)
1103 0 : do iloc=nloc-1,1,-1
1104 0 : tt=arg*tt+cprime(iloc)
1105 : end do
1106 0 : forceloc=xp*tt*rhoel
1107 : else
1108 : forceloc=zero
1109 : end if
1110 : ! PAW: V_PAW^prime-V_L^prime
1111 : else
1112 0 : rr(1)=sqrt(r2)
1113 0 : if (rr(1)>=rzero) then
1114 0 : call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
1115 0 : call calcdVloc_mklocl(dvhgh,rr(1),rloc,charge)
1116 0 : forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1117 : else
1118 0 : forceloc=rhoel*rloc2*dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
1119 : end if
1120 : end if
1121 :
1122 0 : fxgau=fxgau+forceloc*x
1123 0 : fygau=fygau+forceloc*y
1124 0 : fzgau=fzgau+forceloc*z
1125 :
1126 : ! Long range part: error function
1127 0 : Vel=pot(ind)
1128 0 : fxerf=fxerf+xp*Vel*x
1129 0 : fyerf=fyerf+xp*Vel*y
1130 0 : fzerf=fzerf+xp*Vel*z
1131 :
1132 0 : else if (nloc>0) then
1133 0 : r2=x**2+y**2+z**2
1134 0 : arg=r2/rloc**2
1135 0 : xp=exp(-0.5_dp*arg)
1136 0 : tt=cprime(nloc)
1137 0 : do iloc=nloc-1,1,-1
1138 0 : tt=arg*tt+cprime(iloc)
1139 : end do
1140 0 : forceleaked=forceleaked+xp*(1._dp+tt)
1141 : end if
1142 : end do
1143 : end do
1144 : end if
1145 : end do
1146 :
1147 : ! Final result of the forces
1148 0 : floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc**2)*fxgau
1149 0 : floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc**2)*fygau
1150 0 : floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc**2)*fzgau
1151 :
1152 : end do
1153 :
1154 0 : forceleaked=forceleaked*prefactor*hxh*hyh*hzh
1155 0 : if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1156 :
1157 : CONTAINS
1158 : !!***
1159 :
1160 : !----------------------------------------------------------------------
1161 :
1162 : !!****f* mklocl_realspace/calcdVloc_mklocl
1163 : !! NAME
1164 : !! calcdVloc_mklocl
1165 : !!
1166 : !! FUNCTION
1167 : !! Compute 1st-derivative of long-range HGH local ionic potential (derf)
1168 : !!
1169 : !! INPUTS
1170 : !!
1171 : !! OUTPUT
1172 : !!
1173 : !! SOURCE
1174 :
1175 0 : subroutine calcdVloc_mklocl(yy,xx,rloc,Z)
1176 :
1177 : !Arguments ------------------------------------
1178 : !scalars
1179 : real(dp),intent(in) :: xx,rloc,Z
1180 : real(dp),intent(out) :: yy
1181 :
1182 : !Local variables-------------------------------
1183 : !scalars
1184 : real(dp):: arg,tt
1185 :
1186 : ! *************************************************************************
1187 :
1188 0 : arg=xx/(sqrt(2._dp)*rloc)
1189 0 : call derf_ab(tt,arg)
1190 0 : yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1191 :
1192 0 : end subroutine calcdVloc_mklocl
1193 : !!***
1194 :
1195 : !----------------------------------------------------------------------
1196 :
1197 : !!****f* mklocl_wavelets/dvloc_zero_mklocl
1198 : !! NAME
1199 : !! dvloc_zero_mklocl
1200 : !!
1201 : !! FUNCTION
1202 : !! Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0
1203 : !!
1204 : !! INPUTS
1205 : !!
1206 : !! OUTPUT
1207 : !!
1208 : !! SOURCE
1209 :
1210 0 : function dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
1211 :
1212 :
1213 : !Arguments ------------------------------------
1214 : !scalars
1215 : integer,intent(in) :: msz
1216 : real(dp) :: dvloc_zero_mklocl
1217 : real(dp),intent(in) :: charge,rloc
1218 : !arrays
1219 : real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
1220 :
1221 : !Local variables-------------------------------
1222 : !scalars
1223 : real(dp) :: y1,y2,y3,zz=0._dp
1224 : !arrays
1225 : real(dp) :: ll(3),xx(3),yy(3)
1226 :
1227 : ! *************************************************************************
1228 :
1229 : !Select 3 points x1,x2,x3 near 0
1230 0 : if (rad(1)>1.d-10) then
1231 0 : xx(1:3)=rad(1:3)
1232 : else
1233 0 : xx(1:3)=rad(2:4)
1234 : end if
1235 :
1236 : !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
1237 0 : call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1238 0 : call calcdVloc_mklocl(y1,xx(1),rloc,charge)
1239 0 : call calcdVloc_mklocl(y2,xx(2),rloc,charge)
1240 0 : call calcdVloc_mklocl(y3,xx(3),rloc,charge)
1241 0 : yy(1)=(yy(1)-y1)/xx(1)
1242 0 : yy(2)=(yy(2)-y2)/xx(2)
1243 0 : yy(3)=(yy(3)-y3)/xx(3)
1244 :
1245 : !Find a polynomial of the form (z=0):
1246 : !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1247 :
1248 : !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1249 0 : ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1250 : !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1251 0 : ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1252 : !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1253 0 : ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1254 :
1255 0 : dvloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1256 :
1257 0 : end function dvloc_zero_mklocl
1258 : !!***
1259 :
1260 : end subroutine local_forces_new
1261 : !!***
1262 :
1263 : !----------------------------------------------------------------------
1264 :
1265 : !!****f* mklocl_realspace/ind_positions_mklocl
1266 : !! NAME
1267 : !! ind_positions_mklocl
1268 : !!
1269 : !! FUNCTION
1270 : !! determine the index in which the potential must be inserted, following the BC
1271 : !! determine also whether the index is inside or outside the box for free BC
1272 : !!
1273 : !! INPUTS
1274 : !!
1275 : !! OUTPUT
1276 : !!
1277 : !! SOURCE
1278 :
1279 0 : subroutine ind_positions_mklocl(periodic,i,n,j,go)
1280 :
1281 :
1282 : !Arguments -------------------------------
1283 : logical, intent(in) :: periodic
1284 : integer, intent(in) :: i,n
1285 : logical, intent(out) :: go
1286 : integer, intent(out) :: j
1287 :
1288 : ! *********************************************************************
1289 :
1290 0 : if (periodic) then
1291 0 : go=.true.
1292 0 : j=modulo(i-1,n)+1
1293 : else
1294 0 : j=i
1295 0 : go=(i >= 1 .and. i <= n)
1296 : end if
1297 :
1298 0 : end subroutine ind_positions_mklocl
1299 : !!***
1300 :
1301 : !!****f* ABINIT/mklocl_wavelets
1302 : !!
1303 : !! NAME
1304 : !! mklocl_wavelets
1305 : !!
1306 : !! FUNCTION
1307 : !! Compute the ionic local potential when the pseudo-potentials are GTH, using
1308 : !! the special decomposition of these pseudo. The resulting potential is computed with
1309 : !! free boundary conditions. It gives the same result than mklocl_realspace for the
1310 : !! GTH pseudo only with a different way to compute the potential.
1311 : !!
1312 : !! Optionally compute :
1313 : !! option=1 : local ionic potential throughout unit cell
1314 : !! option=2 : contribution of local ionic potential to E gradient wrt xred
1315 : !!
1316 : !! INPUTS
1317 : !! efield (3)=external electric field
1318 : !! mpi_enreg=information about MPI parallelization
1319 : !! natom=number of atoms
1320 : !! nfft=size of vpsp (local potential)
1321 : !! nspden=number of spin-density components
1322 : !! option=type of calculation (potential, forces, ...)
1323 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1324 : !! xcart(3,natom)=cartesian atomic coordinates.
1325 : !! wvl_den=density-potential BigDFT object
1326 : !! wvl_descr=wavelet BigDFT object
1327 : !!
1328 : !! OUTPUT
1329 : !! (if option==1) vpsp(nfft)=the potential resulting from the ionic
1330 : !! density of charge.
1331 : !! (if option==2) grtn(3,natom)=grads of Etot wrt tn. These gradients are in
1332 : !! reduced coordinates. Multiply them by rprimd to get
1333 : !! gradients in cartesian coordinates.
1334 : !!
1335 : !! SOURCE
1336 :
1337 0 : subroutine mklocl_wavelets(efield, grtn, mpi_enreg, natom, nfft, &
1338 0 : & nspden, option, rprimd, vpsp, wvl_den, wvl_descr, xcart)
1339 :
1340 : use defs_wvltypes
1341 : #if defined HAVE_BIGDFT
1342 : use BigDFT_API, only : ELECTRONIC_DENSITY,createIonicPotential,local_forces
1343 : use poisson_solver, only : H_potential
1344 : #endif
1345 :
1346 : !Arguments ------------------------------------
1347 : !scalars
1348 : integer,intent(in) :: option, natom, nfft, nspden
1349 : type(MPI_type),intent(in) :: mpi_enreg
1350 : type(wvl_denspot_type), intent(inout) :: wvl_den
1351 : type(wvl_internal_type), intent(in) :: wvl_descr
1352 : !arrays
1353 : real(dp),intent(in) :: rprimd(3,3),efield(3)
1354 : real(dp),intent(inout) :: grtn(3,natom)
1355 : real(dp), intent(inout) :: vpsp(nfft)
1356 : real(dp),intent(inout) :: xcart(3,natom)
1357 :
1358 : !Local variables-------------------------------
1359 : #if defined HAVE_BIGDFT
1360 : !scalars
1361 : integer :: i,i1,i2,i3,ia,ierr,igeo,me,nproc,shift,comm
1362 : real(dp) :: energ
1363 : character(len=500) :: message
1364 : !arrays
1365 : real(dp) :: epot(3)
1366 : real(dp) :: elecfield(3)=(/zero,zero,zero/) ! Not used here
1367 : real(dp),allocatable :: gxyz(:,:),vhartr(:),rhov(:,:)
1368 : #endif
1369 :
1370 : ! *********************************************************************
1371 :
1372 : #if defined HAVE_BIGDFT
1373 :
1374 : elecfield=zero !not used here
1375 :
1376 : !Manage parallelism
1377 : comm=mpi_enreg%comm_wvl
1378 : nproc=xmpi_comm_size(comm)
1379 : me=xmpi_comm_rank(comm)
1380 :
1381 : !----------------------------------------------------------------------
1382 : ! ----- Option 1: compute local ionic potential -----
1383 : !----------------------------------------------------------------------
1384 : if (option == 1) then
1385 :
1386 : ! We get the kernel for the Poisson solver (used to go from the ionic
1387 : ! charge to the potential).If the kernel is uncomputed, it does it now.
1388 : ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 2, icoulomb, me, kernel, &
1389 : !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder)
1390 : ! if (.not.associated(kernel%co%kernel)) then
1391 : ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 1, icoulomb, me, kernel, &
1392 : !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder)
1393 : ! end if
1394 :
1395 : message=ch10//' mklocl_wavelets: Create local potential from ions.'
1396 : call wrtout(std_out,message,'COLL')
1397 :
1398 : shift = 1 + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1399 : & * wvl_den%denspot%dpbox%i3xcsh
1400 :
1401 : ! Call the BigDFT routine
1402 : call createIonicPotential(wvl_descr%atoms%astruct%geocode, me, nproc, (me == 0), wvl_descr%atoms, &
1403 : & xcart, wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), &
1404 : & wvl_den%denspot%dpbox%hgrids(3), &
1405 : & elecfield, wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, &
1406 : & wvl_den%denspot%dpbox%n3pi, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, &
1407 : & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), &
1408 : & wvl_den%denspot%dpbox%ndims(3), wvl_den%denspot%pkernel, vpsp(shift), 0.d0,wvl_descr%rholoc)
1409 :
1410 : ! Copy vpsp into bigdft object:
1411 : call wvl_rhov_abi2big(1,vpsp,wvl_den%denspot%v_ext,shift=(shift-1))
1412 :
1413 : ! Eventually add the electric field
1414 : if (maxval(efield) > tol12) then
1415 : message=ch10//'mklocl_wavelets: Add the electric field.'
1416 : call wrtout(std_out,message,'COLL')
1417 : ! We add here the electric field since in BigDFT, the field must be on x...
1418 : epot(:) = real(0.5, dp) * efield(:) * wvl_den%denspot%dpbox%hgrids(:)
1419 : do i3 = 1, wvl_den%denspot%dpbox%n3pi, 1
1420 : ia = (i3 - 1) * wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2)
1421 : do i2 = -14, 2 * wvl_descr%Glr%d%n2 + 16, 1
1422 : i = ia + (i2 + 14) * wvl_den%denspot%dpbox%ndims(1)
1423 : do i1 = -14, 2 * wvl_descr%Glr%d%n1 + 16, 1
1424 : i = i + 1
1425 : vpsp(shift + i) = vpsp(shift + i) + &
1426 : & epot(1) * real(i1 - wvl_descr%Glr%d%n1, dp) + &
1427 : & epot(2) * real(i2 - wvl_descr%Glr%d%n2, dp) + &
1428 : & epot(3) * real(i3 - wvl_descr%Glr%d%n3, dp)
1429 : end do
1430 : end do
1431 : end do
1432 : end if
1433 :
1434 : !----------------------------------------------------------------------
1435 : ! ----- Option 2: compute forces induced by local ionic potential -----
1436 : !----------------------------------------------------------------------
1437 : else if (option == 2) then
1438 :
1439 : message=ch10//' mklocl_wavelets: compute local forces.'
1440 : call wrtout(std_out,message,'COLL')
1441 :
1442 : if (wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then
1443 : message='denspot bigdft datstructure should contain rhor!'
1444 : ABI_BUG(message)
1445 : end if
1446 :
1447 : ! Extract density rhor from bigDFT datastructure
1448 : ABI_MALLOC(rhov,(nfft, nspden))
1449 : ABI_MALLOC(vhartr,(nfft))
1450 : rhov=zero ; vhartr=zero
1451 : shift = wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1452 : & * wvl_den%denspot%dpbox%i3xcsh
1453 : do i = 1, min(nfft,size(wvl_den%denspot%rhov)-shift)
1454 : rhov(i, 1) = wvl_den%denspot%rhov(i + shift)
1455 : vhartr(i) = wvl_den%denspot%rhov(i + shift)
1456 : end do
1457 : if (nspden == 2) then
1458 : shift = shift + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1459 : & * wvl_den%denspot%dpbox%n3d
1460 : do i = 1, min(nfft,size(wvl_den%denspot%rhov)-shift)
1461 : rhov(i, 2) = wvl_den%denspot%rhov(i + shift)
1462 : vhartr(i) = vhartr(i) + wvl_den%denspot%rhov(i + shift)
1463 : end do
1464 : end if
1465 :
1466 : ! Compute Hartree potential from rhor
1467 : call H_potential('D',wvl_den%denspot%pkernel,vhartr,vhartr,energ,zero,.false.)
1468 :
1469 : ! Allocate temporary array for forces
1470 : ABI_MALLOC(gxyz,(3, natom))
1471 :
1472 : ! Calculate local part of the forces grtn (modified BigDFT routine)
1473 : call local_forces_wvl(me,natom,xcart,&
1474 : & wvl_den%denspot%dpbox%hgrids(1),&
1475 : & wvl_den%denspot%dpbox%hgrids(2),&
1476 : & wvl_den%denspot%dpbox%hgrids(3),&
1477 : & wvl_descr%Glr%d%n1,wvl_descr%Glr%d%n2,wvl_descr%Glr%d%n3,&
1478 : & wvl_den%denspot%dpbox%n3p,&
1479 : & wvl_den%denspot%dpbox%i3s+wvl_den%denspot%dpbox%i3xcsh,&
1480 : & wvl_den%denspot%dpbox%ndims(1),wvl_den%denspot%dpbox%ndims(2),&
1481 : & rhov,vhartr,gxyz,wvl_descr)
1482 : ! call local_forces(me, wvl_descr%atoms, xcart, &
1483 : ! & wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), &
1484 : ! & wvl_den%denspot%dpbox%hgrids(3), &
1485 : ! & wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, &
1486 : ! & wvl_den%denspot%dpbox%n3p, 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,locstrten,charge)
1489 :
1490 : ! Pending: floc,locstrten and charge are not used here.
1491 : ! Pending: put mpi_enreg%nscatterarr... in object denspot, initialize object, etc.
1492 :
1493 : if (nproc > 1) then
1494 : call xmpi_sum(gxyz, comm, ierr)
1495 : end if
1496 :
1497 : ! Forces should be in reduced coordinates.
1498 : do ia = 1, natom, 1
1499 : do igeo = 1, 3, 1
1500 : grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) &
1501 : & - rprimd(2, igeo) * gxyz(2, ia) &
1502 : & - rprimd(3, igeo) * gxyz(3, ia)
1503 : end do
1504 : end do
1505 :
1506 : ! Deallocate local variables
1507 : ABI_FREE(vhartr)
1508 : ABI_FREE(rhov)
1509 : ABI_FREE(gxyz)
1510 :
1511 : !----------------------------------------------------------------------
1512 :
1513 : else ! option switch
1514 : message = 'Internal error, option should be 1 or 2!'
1515 : ABI_ERROR(message)
1516 : end if
1517 :
1518 : #else
1519 0 : BIGDFT_NOTENABLED_ERROR()
1520 : if (.false.) write(std_out,*) option,natom,nfft,nspden,mpi_enreg%me,&
1521 : & wvl_den%symObj,wvl_descr%h(1),rprimd(1,1),efield(1),grtn(1,1),vpsp(1),xcart(1,1)
1522 : #endif
1523 :
1524 0 : end subroutine mklocl_wavelets
1525 : !!***
1526 :
1527 : !----------------------------------------------------------------------
1528 :
1529 : !!****f* mklocl_wavelets/local_forces_wvl
1530 : !! NAME
1531 : !! local_forces_wvl
1532 : !!
1533 : !! FUNCTION
1534 : !!
1535 : !! INPUTS
1536 : !! hxh,hyh,hzh=wavelet grid spacings
1537 : !! iproc=current MPI process number
1538 : !! n1,n2,n3=number of wavelet points in each direction
1539 : !! n1i,n2i=size of intermediate xy wvl grid
1540 : !! n3pi=number of xy wvl planes handled by current MPI process
1541 : !! i3s=starting index of local potential for current MPI process
1542 : !! natom=number of atoms
1543 : !! pot(*)=Hartree ionic potential
1544 : !! rho(*)=electronic density
1545 : !! rxyz(3,natom)=cartesian coordinates of atoms
1546 : !! wvl=wavelet BigDFT object
1547 : !!
1548 : !! OUTPUT
1549 : !! floc(3,natom)=local ionic potential contribution to forces
1550 : !!
1551 : !! SOURCE
1552 :
1553 : subroutine local_forces_wvl(iproc,natom,rxyz,hxh,hyh,hzh,n1,n2,n3,n3pi,i3s,n1i,n2i,&
1554 : & rho,pot,floc,wvl)
1555 :
1556 : use defs_wvltypes
1557 : #if defined HAVE_BIGDFT
1558 : use BigDFT_API, only : PSPCODE_PAW,ind_positions
1559 : #endif
1560 :
1561 : !Arguments -------------------------------
1562 : !scalars
1563 : integer,intent(in) :: i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom
1564 : real(dp),intent(in) :: hxh,hyh,hzh
1565 : type(wvl_internal_type),intent(in) :: wvl
1566 : !arrays
1567 : real(dp),intent(in) :: rxyz(3,natom)
1568 : real(dp),dimension(*),intent(in) :: rho,pot
1569 : real(dp),intent(out) :: floc(3,natom)
1570 :
1571 : !Local variables -------------------------
1572 : #if defined HAVE_BIGDFT
1573 : !scalars
1574 : integer :: i1,i2,i3,iat,iex,iey,iez,iloc,ind,isx,isy,isz,ityp
1575 : integer :: j1,j2,j3,msz=0
1576 : integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,nloc
1577 : logical :: perx,pery,perz,gox,goy,goz,USE_PAW
1578 : real(dp) :: arg,charge,cutoff,dvhgh,forceloc,fxerf,fyerf,fzerf,fxgau,fygau,fzgau
1579 : real(dp) :: forceleaked,prefactor,r2,rhoel,rloc,rloc2,rx,ry,rz,rzero,tt,vel,x,xp,y,z
1580 : !arrays
1581 : real(dp) :: cprime(4),dvpawdr(1),rr(1)
1582 : real(dp),pointer :: psppar(:,:),rad(:),vloc(:),d2vloc(:)
1583 : #endif
1584 :
1585 : ! *********************************************************************
1586 :
1587 : #if defined HAVE_BIGDFT
1588 :
1589 : if (iproc==0) write(std_out,'(a)',advance='no') ' Calculate local forces...'
1590 :
1591 : !PAW or NCPP ?
1592 : USE_PAW=any(wvl%atoms%npspcode==PSPCODE_PAW)
1593 :
1594 : !Conditions for periodicity in the three directions
1595 : perx=(wvl%atoms%astruct%geocode /= 'F')
1596 : pery=(wvl%atoms%astruct%geocode == 'P')
1597 : perz=(wvl%atoms%astruct%geocode /= 'F')
1598 : call ext_buffers(perx,nbl1,nbr1)
1599 : call ext_buffers(pery,nbl2,nbr2)
1600 : call ext_buffers(perz,nbl3,nbr3)
1601 :
1602 : forceleaked=zero
1603 :
1604 : do iat=1,natom
1605 : ityp=wvl%atoms%astruct%iatype(iat)
1606 :
1607 : ! Coordinates of the center
1608 : rx=rxyz(1,iat)
1609 : ry=rxyz(2,iat)
1610 : rz=rxyz(3,iat)
1611 :
1612 : ! Initialization of the forces
1613 : ! ion-electron term, error function part
1614 : fxerf=zero
1615 : fyerf=zero
1616 : fzerf=zero
1617 : ! ion-electron term, gaussian part
1618 : fxgau=zero
1619 : fygau=zero
1620 : fzgau=zero
1621 :
1622 : ! Building array of coefficients of the derivative of the gaussian part
1623 : psppar => wvl%atoms%psppar(:,:,ityp)
1624 : cprime(1)=2._dp*wvl%atoms%psppar(0,2,ityp)-wvl%atoms%psppar(0,1,ityp)
1625 : cprime(2)=4._dp*wvl%atoms%psppar(0,3,ityp)-wvl%atoms%psppar(0,2,ityp)
1626 : cprime(3)=6._dp*wvl%atoms%psppar(0,4,ityp)-wvl%atoms%psppar(0,3,ityp)
1627 : cprime(4)=-wvl%atoms%psppar(0,4,ityp)
1628 :
1629 : ! Determine number of local terms (HGH pot)
1630 : nloc=0
1631 : do iloc=1,4
1632 : if (wvl%atoms%psppar(0,iloc,ityp).ne.zero) nloc=iloc
1633 : end do
1634 :
1635 : ! Some constants depending on the atom type
1636 : rloc=wvl%atoms%psppar(0,0,ityp) ; rloc2=rloc**2
1637 : charge=real(wvl%atoms%nelpsp(ityp),kind=dp)
1638 : prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1639 :
1640 : ! PAW specifics
1641 : if (USE_PAW) then
1642 : msz=wvl%rholoc%msz(ityp)
1643 : rad => wvl%rholoc%rad(1:msz,ityp)
1644 : vloc => wvl%rholoc%d(1:msz,3,ityp)
1645 : d2vloc => wvl%rholoc%d(1:msz,4,ityp)
1646 : rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2)
1647 : end if
1648 :
1649 : ! Maximum extension of the gaussian
1650 : cutoff=10._dp*rloc
1651 : isx=floor((rx-cutoff)/hxh)
1652 : isy=floor((ry-cutoff)/hyh)
1653 : isz=floor((rz-cutoff)/hzh)
1654 : iex=ceiling((rx+cutoff)/hxh)
1655 : iey=ceiling((ry+cutoff)/hyh)
1656 : iez=ceiling((rz+cutoff)/hzh)
1657 :
1658 : ! Calculate the forces near the atom due to the gaussian
1659 : ! and error function parts of the potential
1660 : if (n3pi>0) then
1661 : do i3=isz,iez
1662 : z=real(i3,kind=dp)*hzh-rz
1663 : call ind_positions(perz,i3,n3,j3,goz)
1664 : j3=j3+nbl3+1
1665 : do i2=isy,iey
1666 : y=real(i2,kind=dp)*hyh-ry
1667 : call ind_positions(pery,i2,n2,j2,goy)
1668 : do i1=isx,iex
1669 : x=real(i1,kind=dp)*hxh-rx
1670 : call ind_positions(perx,i1,n1,j1,gox)
1671 :
1672 : r2=x**2+y**2+z**2
1673 : xp=exp(-0.5_dp*r2/rloc2)
1674 :
1675 : if ((j3>=i3s.and.j3<=i3s+n3pi-1) .and. (gox.and.goy)) then
1676 : ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s+1-1)*n1i*n2i
1677 :
1678 : ! Short range part
1679 : rhoel=rho(ind)
1680 : ! HGH: V_S^prime=gaussian
1681 : if (.not.USE_PAW) then
1682 : if (nloc/=0) then
1683 : arg=r2/rloc2
1684 : tt=cprime(nloc)
1685 : do iloc=nloc-1,1,-1
1686 : tt=arg*tt+cprime(iloc)
1687 : end do
1688 : forceloc=xp*tt*rhoel
1689 : else
1690 : forceloc=zero
1691 : end if
1692 : ! PAW: V_PAW^prime-V_L^prime
1693 : else
1694 : rr(1)=sqrt(r2)
1695 :
1696 : if (rr(1)>=rzero) then
1697 : call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
1698 : call calcdVloc_wvl(dvhgh,rr(1),rloc,charge)
1699 : forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1700 : else
1701 : forceloc=rhoel*rloc2*dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc)
1702 : end if
1703 : end if
1704 :
1705 : fxgau=fxgau+forceloc*x
1706 : fygau=fygau+forceloc*y
1707 : fzgau=fzgau+forceloc*z
1708 :
1709 : ! Long range part: error function
1710 : vel=pot(ind)
1711 : fxerf=fxerf+xp*vel*x
1712 : fyerf=fyerf+xp*vel*y
1713 : fzerf=fzerf+xp*vel*z
1714 :
1715 : else if ((.not.goz).and.(nloc>0)) then
1716 : arg=r2/rloc2
1717 : tt=cprime(nloc)
1718 : do iloc=nloc-1,1,-1
1719 : tt=arg*tt+cprime(iloc)
1720 : end do
1721 : forceleaked=forceleaked+prefactor*xp*tt*rho(1)
1722 : end if
1723 :
1724 : end do ! i1
1725 : end do ! i2
1726 : end do ! i3
1727 : end if ! n3pi>0
1728 :
1729 : ! Final result of the forces
1730 : floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc2)*fxgau
1731 : floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc2)*fygau
1732 : floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc2)*fzgau
1733 :
1734 : end do ! iat
1735 :
1736 : forceleaked=forceleaked*hxh*hyh*hzh
1737 : if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1738 :
1739 : #else
1740 : BIGDFT_NOTENABLED_ERROR()
1741 : if (.false.) write(std_out,*) i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom,hxh,hyh,hzh,&
1742 : & rxyz(1,1),floc(1,1),rho(1),pot(1),wvl%h(1)
1743 : #endif
1744 :
1745 : CONTAINS
1746 : !!***
1747 :
1748 : !----------------------------------------------------------------------
1749 :
1750 : !!****f* mklocl_wavelets/calcdVloc_wvl
1751 : !! NAME
1752 : !! calcdVloc_wvl
1753 : !!
1754 : !! FUNCTION
1755 : !! Compute 1st-derivative of long-range HGH local ionic potential (derf)
1756 : !!
1757 : !! INPUTS
1758 : !!
1759 : !! OUTPUT
1760 : !!
1761 : !! SOURCE
1762 :
1763 : subroutine calcdVloc_wvl(yy,xx,rloc,Z)
1764 :
1765 : !Arguments ------------------------------------
1766 : !scalars
1767 : real(dp),intent(in) :: xx,rloc,Z
1768 : real(dp),intent(out) :: yy
1769 :
1770 : !Local variables-------------------------------
1771 : !scalars
1772 : real(dp):: arg,tt
1773 :
1774 : ! *************************************************************************
1775 :
1776 : arg=xx/(sqrt(2._dp)*rloc)
1777 : call derf_ab(tt,arg)
1778 : yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1779 :
1780 : end subroutine calcdVloc_wvl
1781 : !!***
1782 :
1783 : !----------------------------------------------------------------------
1784 :
1785 : !!****f* mklocl_wavelets/dvloc_zero_wvl
1786 : !! NAME
1787 : !! dvloc_zero_wvl
1788 : !!
1789 : !! FUNCTION
1790 : !! Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0
1791 : !!
1792 : !! INPUTS
1793 : !!
1794 : !! OUTPUT
1795 : !!
1796 : !! SOURCE
1797 :
1798 : function dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc)
1799 :
1800 :
1801 : !Arguments ------------------------------------
1802 : !scalars
1803 : integer,intent(in) :: msz
1804 : real(dp) :: dvloc_zero_wvl
1805 : real(dp),intent(in) :: charge,rloc
1806 : !arrays
1807 : real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
1808 :
1809 : !Local variables-------------------------------
1810 : !scalars
1811 : real(dp) :: y1,y2,y3,zz=0._dp
1812 : !arrays
1813 : real(dp) :: ll(3),xx(3),yy(3)
1814 :
1815 : ! *************************************************************************
1816 :
1817 : !Select 3 points x1,x2,x3 near 0
1818 : if (rad(1)>1.d-10) then
1819 : xx(1:3)=rad(1:3)
1820 : else
1821 : xx(1:3)=rad(2:4)
1822 : end if
1823 :
1824 : !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
1825 : call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1826 : call calcdVloc_wvl(y1,xx(1),rloc,charge)
1827 : call calcdVloc_wvl(y2,xx(2),rloc,charge)
1828 : call calcdVloc_wvl(y3,xx(3),rloc,charge)
1829 : yy(1)=(yy(1)-y1)/xx(1)
1830 : yy(2)=(yy(2)-y2)/xx(2)
1831 : yy(3)=(yy(3)-y3)/xx(3)
1832 :
1833 : !Find a polynomial of the form (z=0):
1834 : !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1835 :
1836 : !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1837 : ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1838 : !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1839 : ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1840 : !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1841 : ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1842 :
1843 : dvloc_zero_wvl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1844 :
1845 : end function dvloc_zero_wvl
1846 : !!***
1847 :
1848 : end subroutine local_forces_wvl
1849 : !!***
1850 :
1851 : end module m_mklocl_realspace
1852 : !!***
|
