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 : :
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 : : !!***
|
