TABLE OF CONTENTS


ABINIT/m_mklocl_realspace [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mklocl_realspace

FUNCTION

   Routines related to the local part of the pseudopotentials.
   Computation is done in real space (useful for isolated systems).

COPYRIGHT

  Copyright (C) 2013-2022 ABINIT group (TRangel, MT, DC)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

TODO

  This module could be merged with m_mklocl

SOURCE

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

ABINIT/mklocl_realspace [ Functions ]

[ Top ] [ Functions ]

NAME

  mklocl_realspace

FUNCTION

 This method is equivalent to mklocl_recipspace except that
 it uses real space pseudo-potentials. It is useful for isolated
 systems. Then the option 3 and 4 are not available for this implementation.

 Optionally compute:
  option=1 : local ionic potential throughout unit cell
  option=2 : contribution of local ionic potential to E gradient wrt xred

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in unit cell.
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nspden=number of spin-density components
  ntypat=number of types of atoms.
  option= (see above)
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$)
  rhor(nfft,nspden)=electron density in electrons/bohr**3.
    (needed if option==2 or if option==4)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  ucvol=unit cell volume ($\textrm{Bohr}^3$).
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  (if option==1) vpsp(nfft)=local crystal pseudopotential in real space.
  (if option==2) grtn(3,natom)=grads of Etot wrt tn. These gradients are in
                 reduced coordinates. Multiply them by rprimd to get
                 gradients in cartesian coordinates.

SOURCE

 99 subroutine mklocl_realspace(grtn,icoulomb,mpi_enreg,natom,nattyp,nfft,ngfft,nscforder, &
100 &                           nspden,ntypat,option,pawtab,psps,rhog,rhor,rprimd,typat,&
101 &                           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  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
142  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
143  real(dp) :: coord(3),coordXYZ(3),refValue(3),tsec(2)
144  real(dp),allocatable :: coordCart_interpol(:,:),coordRed_interpol(:,:)
145  real(dp),allocatable :: gridcart(:,:)
146  real(dp),allocatable :: grtn_cart_interpol(:,:),grtn_diff(:,:)
147  real(dp),allocatable :: rhog_interpol(:,:),rhog_testing(:,:),rhor_interpol(:)
148  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  if (option==2) then
155    call timab(72,1,tsec)
156  end if
157 
158 !Several constants (FFT sizes and parallelism)
159  n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
160  nproc_fft = ngfft(10) ;  me_fft = ngfft(11)
161  n3d = ngfft(13)          !for parallel runs
162  if (nproc_fft==1) n3d=n3  !for serial runs
163  comm_fft=mpi_enreg%comm_fft
164  if(me_fft /= mpi_enreg%me_fft .or. nproc_fft /= mpi_enreg%nproc_fft) then
165    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  geocode='P'
170  if (icoulomb==1) geocode='F'
171  if (icoulomb==2) geocode='S'
172 
173 !Get the distrib associated with this fft_grid
174  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
175 
176 !Store xcart for each atom
177  ABI_MALLOC(xcart,(3, natom))
178  call xred2xcart(natom, rprimd, xcart, xred)
179 !Store cartesian coordinates for each grid points
180  ABI_MALLOC(gridcart,(3, nfft))
181  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  doIt=.true.
185  do ii=1,psps%npsp
186    doIt=doIt .and.&
187    (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  if (doIt) then
192 
193 !  Definition of the grid spacings as in the kernel routine
194    hgx = rprimd(1,1)/(n1)
195    hgy = rprimd(2,2)/(n2)
196    hgz = rprimd(3,3)/(n3)
197 
198 
199 !----------------------------------------------------------------------
200 ! ----- Option 1: compute local ionic potential                   -----
201 !----------------------------------------------------------------------
202    if (option==1) then
203 
204      call psolver_kernel( (/ hgx, hgy, hgz /), 2, icoulomb, me_fft, kernel, comm_fft, &
205 &     (/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 &     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    else if (option == 2) then
217 
218 !    Compute Hartree potential from rhor
219      ABI_MALLOC(vhartr,(nfft))
220      call psolver_hartree(entmp, (/ hgx, hgy, hgz /), icoulomb, me_fft, comm_fft, nfft, &
221 &     (/n1,n2,n3/), nproc_fft, nscforder, nspden, rhor, vhartr, usewvl)
222 
223 !    Allocate temporary array for forces
224      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 &     hgx,hgy,hgz, n1,n2,n3,n3d, rhor,vhartr, gxyz, pawtab,psps%usepaw)
230 
231 !    Forces should be in reduced coordinates.
232      do ia = 1, natom, 1
233        do igeo = 1, 3, 1
234          grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) &
235 &         - rprimd(2, igeo) * gxyz(2, ia) &
236 &         - rprimd(3, igeo) * gxyz(3, ia)
237        end do
238      end do
239 
240 !    Deallocate local variables
241      ABI_FREE(vhartr)
242      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    dr = psps%qgrid_vl(2)
257    invdr = 1._dp / dr
258    dr2div6 = dr * dr / 6._dp
259 
260    if (option == 1) then
261 !    Set 0 in vpsp before summing
262      vpsp(:) = 0._dp
263    else if (option == 2) then
264 !    Allocate array to store cartesian gradient computed with
265 !    an interpolation of rhor
266      ABI_MALLOC(grtn_cart_interpol,(3, natom))
267      grtn_cart_interpol(:, :) = 0._dp
268 
269      n_interpol = nStep ** 3
270      ABI_MALLOC(coordRed_interpol,(3, nStep ** 3))
271      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      ABI_MALLOC(rhor_interpol,(nfft * n_interpol))
280      ABI_MALLOC(rhor_work,(nfft * n_interpol))
281      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      rhog_interpol(:, :) = 0._dp
298      ii = 0
299      do i3 = 1, n3, 1
300        if (i3 <= n3 / 2) then
301          iz = i3
302        else
303          iz = n3 * nStep - n3 + i3
304        end if
305        do i2 = 1, n2, 1
306          if (i2 <= n2 / 2) then
307            iy = i2
308          else
309            iy = n2 * nStep - n2 + i2
310          end if
311          do i1 = 1, n1, 1
312            ii = ii + 1
313            if (i1 <= n1 / 2) then
314              ix = i1
315            else
316              ix = n1 * nStep - n1 + i1
317            end if
318            jj = (iz - 1) * n2 * n1 * nStep ** 2 + (iy - 1) * n3 * nStep + ix
319            if (testing .and. customRho) then
320              rhog_interpol(:, jj) = rhog_testing(:, ii)
321            else
322              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      ngfft_interpol(:) = ngfft(:)
330      ngfft_interpol(1:3) = (/ n1 * nStep, n2 * nStep, n3 * nStep /)
331      ngfft_interpol(4:6) = (/ n1 * nStep + 1, n2 * nStep + 1, n3 * nStep /)
332      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      jj = 0
336      do i3 = 1, n3, 1
337        do i2 = 1, n2, 1
338          do i1 = 1, n1, 1
339            do iz = 1, nStep, 1
340              do iy = 1, nStep, 1
341                do ix = 1, nStep, 1
342                  jj = jj + 1
343                  kk = ((i3 - 1) * nStep + iz - 1) ! z coordinate in the interpolated grid
344                  kk = kk * n1 * n2 * nStep ** 2
345                  kk = kk + ((i2 - 1) * nStep + iy - 1) * n1 * nStep ! adding y coordinate
346                  kk = kk + (i1 - 1) * nStep + ix ! adding x coordinate
347                  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      ABI_FREE(rhor_work)
355 
356 !    Compute grid access in the interpolated volume
357      ii = 0
358      do iz = 1, nStep, 1
359        z = real(iz - 1, dp) / real(nStep, dp)
360        do iy = 1, nStep, 1
361          y = real(iy - 1, dp) / real(nStep, dp)
362          do ix = 1, nStep, 1
363            x = real(ix - 1, dp) / real(nStep, dp)
364            ii = ii + 1
365            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 &           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    ia1=1
390    do itypat = 1, ntypat, 1
391 !    ia1,ia2 sets range of loop over atoms:
392      ia2 = ia1 + nattyp(itypat) - 1
393 
394      do ii = 1, nfft, 1
395        do ia = ia1, ia2, 1
396          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 &           (gridcart(3, ii) - xcart(3, ia)) ** 2)
402 
403 !          Coefficients needed to compute the spline.
404            jj = int(r * invdr) + 1
405            if (jj > psps%mqgrid_vl - 2) then
406              write(message, '(3a,i0,a,i0,a,a)' )&
407 &             '  pseudo-potential local part sampling is not wide enough', ch10, &
408 &             '  want to access position ', jj, ' whereas mqgrid_vl = ', psps%mqgrid_vl, ch10, &
409 &             '  Action : no idea, contact developpers...'
410              ABI_ERROR(message)
411            end if
412            delta = r - psps%qgrid_vl(jj)
413            bb = delta * invdr
414            aa = 1._dp - bb
415            cc = aa * (aa ** 2 - 1._dp) * dr2div6
416            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 &           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            vpsp(ii) = vpsp(ii) + deltaV
424          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            ind = (ii - 1) * n_interpol
456            do kk = 1, n_interpol, 1
457              ind = ind + 1
458 
459              if (rhor_interpol(ind) > min_rho_value) then
460 !              Assume orthogonal box...
461                coordXYZ(1) = gridcart(1, ii) - xcart(1, ia) + coordCart_interpol(1, kk)
462                coordXYZ(2) = gridcart(2, ii) - xcart(2, ia) + coordCart_interpol(2, kk)
463                coordXYZ(3) = gridcart(3, ii) - xcart(3, ia) + coordCart_interpol(3, kk)
464                r = coordXYZ(1) ** 2 + coordXYZ(2) ** 2 + coordXYZ(3) ** 2
465 
466                if (r /= 0._dp) then
467                  r = sqrt(r)
468 !                Coefficients needed to compute the spline.
469                  jj = int(r * invdr) + 1
470                  delta = r - psps%qgrid_vl(jj)
471                  bb = delta * invdr
472                  aa = 1._dp - bb
473                  cc = aa * (aa ** 2 - 1._dp) * dr2div6
474                  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 &                 dd * psps%dvlspl(jj + 1, 2, itypat)
479                  deltaV = deltaV * rhor_interpol(ind) / r
480                  grtn_cart_interpol(1, ia) = grtn_cart_interpol(1, ia) - deltaV * coordXYZ(1)
481                  grtn_cart_interpol(2, ia) = grtn_cart_interpol(2, ia) - deltaV * coordXYZ(2)
482                  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      ia1 = ia2 + 1
546    end do
547 !  End loop over type of atoms
548 
549    if(option==2)then
550 !    multiply the forces by the volume of a single box mesh.
551      grtn_cart_interpol(:, :) = grtn_cart_interpol(:, :) * &
552 &     ucvol / real(n1 * n2 * n3, dp) * vol_interpol
553 !    Transform cartesian forces to reduce coordinates
554      do ia = 1, natom, 1
555        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 &         rprimd(3, igeo) * grtn_cart_interpol(3, ia)
559        end do
560      end do
561      ABI_FREE(rhor_interpol)
562      ABI_FREE(rhog_interpol)
563      ABI_FREE(coordRed_interpol)
564      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  ABI_FREE(xcart)
584  ABI_FREE(gridcart)
585 
586 !Close timing counters
587  if (option==2)then
588    call timab(72,2,tsec)
589  end if
590 
591 end subroutine mklocl_realspace

ABINIT/mklocl_wavelets [ Functions ]

[ Top ] [ Functions ]

NAME

 mklocl_wavelets

FUNCTION

 Compute the ionic local potential when the pseudo-potentials are GTH, using
 the special decomposition of these pseudo. The resulting potential is computed with
 free boundary conditions. It gives the same result than mklocl_realspace for the
 GTH pseudo only with a different way to compute the potential.

 Optionally compute :
  option=1 : local ionic potential throughout unit cell
  option=2 : contribution of local ionic potential to E gradient wrt xred

INPUTS

  efield (3)=external electric field
  mpi_enreg=information about MPI parallelization
  natom=number of atoms
  nfft=size of vpsp (local potential)
  nspden=number of spin-density components
  option=type of calculation (potential, forces, ...)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  xcart(3,natom)=cartesian atomic coordinates.
  wvl_den=density-potential BigDFT object
  wvl_descr=wavelet BigDFT object

OUTPUT

  (if option==1) vpsp(nfft)=the potential resulting from the ionic
                 density of charge.
  (if option==2) grtn(3,natom)=grads of Etot wrt tn. These gradients are in
                 reduced coordinates. Multiply them by rprimd to get
                 gradients in cartesian coordinates.

SOURCE

1337 subroutine mklocl_wavelets(efield, grtn, mpi_enreg, natom, nfft, &
1338      & 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  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 end subroutine mklocl_wavelets

mklocl_realspace/calcdVloc_mklocl [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  calcdVloc_mklocl

FUNCTION

  Compute 1st-derivative of long-range HGH local ionic potential (derf)

INPUTS

OUTPUT

SOURCE

1175 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      arg=xx/(sqrt(2._dp)*rloc)
1189      call derf_ab(tt,arg)
1190      yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1191 
1192    end subroutine calcdVloc_mklocl

mklocl_realspace/calcVloc_mklocl [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  calcVloc_mklocl

FUNCTION

INPUTS

OUTPUT

SOURCE

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

mklocl_realspace/createIonicPotential_new [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  createIonicPotential_new

FUNCTION

INPUTS

OUTPUT

SOURCE

608 subroutine createIonicPotential_new(fftn3_distrib,ffti3_local,geocode,iproc,&
609 &  nproc,nat,ntypes,iatype,psppar,nelpsp,rxyz,gridcart,&
610 &  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  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

mklocl_realspace/ind_positions_mklocl [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  ind_positions_mklocl

FUNCTION

 determine the index in which the potential must be inserted, following the BC
 determine also whether the index is inside or outside the box for free BC

INPUTS

OUTPUT

SOURCE

1279 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      if (periodic) then
1291        go=.true.
1292        j=modulo(i-1,n)+1
1293      else
1294        j=i
1295        go=(i >= 1 .and. i <= n)
1296      end if
1297 
1298    end subroutine ind_positions_mklocl

mklocl_realspace/local_forces_new [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  local_forces_new

FUNCTION

INPUTS

OUTPUT

SOURCE

 979 subroutine local_forces_new(fftn3_distrib,ffti3_local,&
 980      geocode,iproc,ntypes,nat,iatype,rxyz,gridcart,psppar,nelpsp,hxh,hyh,hzh,&
 981      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  real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
1010 
1011 ! *********************************************************************
1012 
1013    if (iproc == 0) write(std_out,'(1x,a)',advance='no')'Calculate local forces...'
1014 
1015 !Conditions for periodicity in the three directions
1016    perx=(geocode /= 'F')
1017    pery=(geocode == 'P')
1018    perz=(geocode /= 'F')
1019 
1020    forceleaked=zero
1021 
1022    do iat=1,nat
1023      ityp=iatype(iat)
1024 !  Coordinates of the center
1025      rx=rxyz(1,iat)
1026      ry=rxyz(2,iat)
1027      rz=rxyz(3,iat)
1028 
1029 !  Initialization of the forces
1030 !  ion-electron term, error function part
1031      fxerf=zero
1032      fyerf=zero
1033      fzerf=zero
1034 !  ion-electron term, gaussian part
1035      fxgau=zero
1036      fygau=zero
1037      fzgau=zero
1038 
1039 !  Building array of coefficients of the derivative of the gaussian part
1040      cprime(1)=2._dp*psppar(0,2,ityp)-psppar(0,1,ityp)
1041      cprime(2)=4._dp*psppar(0,3,ityp)-psppar(0,2,ityp)
1042      cprime(3)=6._dp*psppar(0,4,ityp)-psppar(0,3,ityp)
1043      cprime(4)=-psppar(0,4,ityp)
1044 
1045 !  Determine number of local terms (HGH pot)
1046      nloc=0
1047      do iloc=1,4
1048        if (psppar(0,iloc,ityp).ne.zero) nloc=iloc
1049      end do
1050 
1051 !  Some constants depending on the atom type
1052      rloc=psppar(0,0,ityp) ; rloc2=rloc**2
1053      charge=real(nelpsp(ityp),kind=dp)
1054      prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1055 
1056 !  PAW specifics
1057      if (usepaw==1) then
1058        msz=pawtab(ityp)%wvl%rholoc%msz
1059        rad    => pawtab(ityp)%wvl%rholoc%rad(1:msz)
1060        vloc   => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
1061        d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
1062        rzero=rad(1);if (rad(1)<=1.d-10) rzero=rad(2)
1063      end if
1064 
1065 !  Maximum extension of the gaussian
1066      cutoff=10._dp*rloc
1067      isx=floor((rx-cutoff)/hxh)
1068      isy=floor((ry-cutoff)/hyh)
1069      isz=floor((rz-cutoff)/hzh)
1070      iex=ceiling((rx+cutoff)/hxh)
1071      iey=ceiling((ry+cutoff)/hyh)
1072      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      do i3=isz,iez
1077        z=real(i3,kind=dp)*hzh-rz
1078        call ind_positions_mklocl(perz,i3,n3,j3,goz)
1079        if(fftn3_distrib(j3)==iproc) then
1080          i3loc=ffti3_local(j3)
1081          do i2=isy,iey
1082            y=real(i2,kind=dp)*hyh-ry
1083            call ind_positions_mklocl(pery,i2,n2,j2,goy)
1084            do i1=isx,iex
1085              x=real(i1,kind=dp)*hxh-rx
1086              call ind_positions_mklocl(perx,i1,n1,j1,gox)
1087 
1088              if (goz.and.goy.and.gox) then
1089                ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2
1090                x=(gridcart(1,ind)-rx)
1091                y=(gridcart(2,ind)-ry)
1092                z=(gridcart(3,ind)-rz)
1093                r2=x**2+y**2+z**2
1094                xp=exp(-0.5_dp*r2/rloc**2)
1095 
1096 !            Short range part
1097                rhoel=rho(ind)
1098 !            HGH: V_S^prime=gaussian
1099                if (usepaw==0) then
1100                  if (nloc/=0) then
1101                    arg=r2/rloc**2
1102                    tt=cprime(nloc)
1103                    do iloc=nloc-1,1,-1
1104                      tt=arg*tt+cprime(iloc)
1105                    end do
1106                    forceloc=xp*tt*rhoel
1107                  else
1108                    forceloc=zero
1109                  end if
1110 !            PAW: V_PAW^prime-V_L^prime
1111                else
1112                  rr(1)=sqrt(r2)
1113                  if (rr(1)>=rzero) then
1114                    call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
1115                    call calcdVloc_mklocl(dvhgh,rr(1),rloc,charge)
1116                    forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1117                  else
1118                    forceloc=rhoel*rloc2*dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
1119                  end if
1120                end if
1121 
1122                fxgau=fxgau+forceloc*x
1123                fygau=fygau+forceloc*y
1124                fzgau=fzgau+forceloc*z
1125 
1126 !            Long range part: error function
1127                Vel=pot(ind)
1128                fxerf=fxerf+xp*Vel*x
1129                fyerf=fyerf+xp*Vel*y
1130                fzerf=fzerf+xp*Vel*z
1131 
1132              else if (nloc>0) then
1133                r2=x**2+y**2+z**2
1134                arg=r2/rloc**2
1135                xp=exp(-0.5_dp*arg)
1136                tt=cprime(nloc)
1137                do iloc=nloc-1,1,-1
1138                  tt=arg*tt+cprime(iloc)
1139                end do
1140                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      floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc**2)*fxgau
1149      floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc**2)*fygau
1150      floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc**2)*fzgau
1151 
1152    end do
1153 
1154    forceleaked=forceleaked*prefactor*hxh*hyh*hzh
1155    if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1156 
1157    CONTAINS

mklocl_realspace/vloc_zero_mklocl [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  vloc_zero_mklocl

FUNCTION

  Use a quadratic interpolation to get limit of Vloc(x) at x->0

INPUTS

OUTPUT

SOURCE

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

mklocl_wavelets/calcdVloc_wvl [ Functions ]

[ Top ] [ mklocl_wavelets ] [ Functions ]

NAME

  calcdVloc_wvl

FUNCTION

  Compute 1st-derivative of long-range HGH local ionic potential (derf)

INPUTS

OUTPUT

SOURCE

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

mklocl_wavelets/dvloc_zero_mklocl [ Functions ]

[ Top ] [ mklocl_wavelets ] [ Functions ]

NAME

  dvloc_zero_mklocl

FUNCTION

  Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0

INPUTS

OUTPUT

SOURCE

1210 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      if (rad(1)>1.d-10) then
1231        xx(1:3)=rad(1:3)
1232      else
1233        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      call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1238      call calcdVloc_mklocl(y1,xx(1),rloc,charge)
1239      call calcdVloc_mklocl(y2,xx(2),rloc,charge)
1240      call calcdVloc_mklocl(y3,xx(3),rloc,charge)
1241      yy(1)=(yy(1)-y1)/xx(1)
1242      yy(2)=(yy(2)-y2)/xx(2)
1243      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      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      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      ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1254 
1255      dvloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1256 
1257    end function dvloc_zero_mklocl

mklocl_wavelets/dvloc_zero_wvl [ Functions ]

[ Top ] [ mklocl_wavelets ] [ Functions ]

NAME

  dvloc_zero_wvl

FUNCTION

  Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0

INPUTS

OUTPUT

SOURCE

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

mklocl_wavelets/local_forces_wvl [ Functions ]

[ Top ] [ mklocl_wavelets ] [ Functions ]

NAME

  local_forces_wvl

FUNCTION

INPUTS

  hxh,hyh,hzh=wavelet grid spacings
  iproc=current MPI process number
  n1,n2,n3=number of wavelet points in each direction
  n1i,n2i=size of intermediate xy wvl grid
  n3pi=number of xy wvl planes handled by current MPI process
  i3s=starting index of local potential for current MPI process
  natom=number of atoms
  pot(*)=Hartree ionic potential
  rho(*)=electronic density
  rxyz(3,natom)=cartesian coordinates of atoms
  wvl=wavelet BigDFT object

OUTPUT

  floc(3,natom)=local ionic potential contribution to forces

SOURCE

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