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-2018 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

PARENTS

CHILDREN

SOURCE

25 #if defined HAVE_CONFIG_H
26 #include "config.h"
27 #endif
28 
29 #include "abi_common.h"
30 
31 module m_mklocl_realspace
32 
33  use defs_basis
34  use defs_datatypes
35  use defs_abitypes
36  use defs_wvltypes
37  use m_xmpi
38  use m_abicore
39  use m_errors
40 
41  use m_time,        only : timab
42  use m_geometry,    only : xred2xcart
43  use m_fft_mesh,    only : mkgrid_fft
44  use m_mpinfo,      only : ptabs_fourdp
45  use m_pawtab,      only : pawtab_type
46  use m_paw_numeric, only : paw_splint, paw_splint_der
47  use m_psolver,     only : psolver_hartree, psolver_kernel
48  use m_abi2big,     only : wvl_rhov_abi2big
49  use m_wvl_wfs,     only : derf_ab
50  use m_fft,         only : fourdp
51 
52 
53  implicit none
54 
55  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.

PARENTS

      mklocl

CHILDREN

SOURCE

110 subroutine mklocl_realspace(grtn,icoulomb,mpi_enreg,natom,nattyp,nfft,ngfft,nscforder, &
111 &                           nspden,ntypat,option,pawtab,psps,rhog,rhor,rprimd,typat,&
112 &                           ucvol,usewvl,vpsp,xred)
113 
114 #if defined HAVE_BIGDFT
115  use BigDFT_API,    only : coulomb_operator,deallocate_coulomb_operator
116  use defs_PSolver
117 #else
118  use defs_wvltypes, only : coulomb_operator
119 #endif
120 
121 !This section has been created automatically by the script Abilint (TD).
122 !Do not modify the following lines by hand.
123 #undef ABI_FUNC
124 #define ABI_FUNC 'mklocl_realspace'
125 !End of the abilint section
126 
127  implicit none
128 
129 !Arguments ------------------------------------
130 !scalars
131  integer,intent(in) :: natom,nfft,nspden,ntypat,option
132  real(dp),intent(in) :: ucvol
133  type(MPI_type),intent(in) :: mpi_enreg
134  type(pseudopotential_type),intent(in) :: psps
135  type(pawtab_type),intent(in)  :: pawtab(ntypat*psps%usepaw)
136 !arrays
137  integer,intent(in)  :: icoulomb,nscforder,usewvl
138  integer,intent(in)  :: nattyp(ntypat),ngfft(18),typat(natom)
139  real(dp),intent(in) :: rhog(2,nfft)
140  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3)
141  real(dp),intent(in) :: xred(3,natom)
142  real(dp),intent(out) :: grtn(3,natom),vpsp(nfft)
143 
144 !Local variables-------------------------------
145  character(len=1) :: geocode
146   !testing variables
147  !scalars
148  integer,parameter :: nStep=2
149  integer :: comm_fft,countParSeconde,i1,i2,i3
150  integer :: ia,ia1,ia2,igeo,ii,ind,itypat,ix,iy,iz,jj
151  integer :: kk,me_fft,n1,n2,n3,n3d,n_interpol
152  integer :: nproc_fft,tpsStart,tpsStop
153  real(dp),parameter :: min_rho_value=1.0d-12
154  real(dp) :: aa,bb,cc,dd,delta,deltaV,dr,dr2div6,invdr,r,vol_interpol,x,y,z,hgx,hgy,hgz,entmp
155  logical,parameter :: customRho=.false.,finiteDiff=.false.,testing=.false.
156  logical :: doIt
157  character(len=500) :: message
158 !arrays
159  integer :: ngfft_interpol(18)
160  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
161  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
162  real(dp) :: coord(3),coordXYZ(3),refValue(3),tsec(2)
163  real(dp),allocatable :: coordCart_interpol(:,:),coordRed_interpol(:,:)
164  real(dp),allocatable :: gridcart(:,:)
165  real(dp),allocatable :: grtn_cart_interpol(:,:),grtn_diff(:,:)
166  real(dp),allocatable :: rhog_interpol(:,:),rhog_testing(:,:),rhor_interpol(:)
167  real(dp),allocatable :: rhor_testing(:),rhor_work(:),xcart(:,:),vhartr(:),gxyz(:,:)
168  type(coulomb_operator):: kernel
169 
170 ! *************************************************************************
171 
172 !Keep track of total time spent here
173  if (option==2) then
174    call timab(72,1,tsec)
175  end if
176 
177 !Several constants (FFT sizes and parallelism)
178  n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
179  nproc_fft = ngfft(10) ;  me_fft = ngfft(11)
180  n3d = ngfft(13)          !for parallel runs
181  if (nproc_fft==1) n3d=n3  !for serial runs
182  comm_fft=mpi_enreg%comm_fft
183  if(me_fft /= mpi_enreg%me_fft .or. nproc_fft /= mpi_enreg%nproc_fft) then
184    MSG_BUG("mpi_enreg%x_fft not equal to the corresponding values in ngfft")
185  end if
186 
187 !Conditions for periodicity in the three directions
188  geocode='P'
189  if (icoulomb==1) geocode='F'
190  if (icoulomb==2) geocode='S'
191 
192 !Get the distrib associated with this fft_grid
193  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
194 
195 !Store xcart for each atom
196  ABI_ALLOCATE(xcart,(3, natom))
197  call xred2xcart(natom, rprimd, xcart, xred)
198 !Store cartesian coordinates for each grid points
199  ABI_ALLOCATE(gridcart,(3, nfft))
200  call mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd)
201 
202 !Check whether all the PSP considered are of type GTH-HGH or PAW
203  doIt=.true.
204  do ii=1,psps%npsp
205    doIt=doIt .and.&
206    (psps%pspcod(ii)==2.or.psps%pspcod(ii)==3.or.psps%pspcod(ii)==10.or.psps%pspcod(ii)==7)
207  end do
208 
209 !HGH-GTH/PAW treatment presumably starts here
210  if (doIt) then
211 
212 !  Definition of the grid spacings as in the kernel routine
213    hgx = rprimd(1,1)/(n1)
214    hgy = rprimd(2,2)/(n2)
215    hgz = rprimd(3,3)/(n3)
216 
217 
218 !----------------------------------------------------------------------
219 ! ----- Option 1: compute local ionic potential                   -----
220 !----------------------------------------------------------------------
221    if (option==1) then
222 
223      call psolver_kernel( (/ hgx, hgy, hgz /), 2, icoulomb, me_fft, kernel, comm_fft, &
224 &     (/n1,n2,n3/), nproc_fft, nscforder)
225 
226      call createIonicPotential_new(fftn3_distrib,ffti3_local,&
227 &     geocode,me_fft, nproc_fft, natom, &
228 &     ntypat, typat, psps%gth_params%psppar, &
229 &     int(psps%ziontypat), xcart,gridcart, hgx,hgy,hgz, &
230 &     n1,n2,n3d,n3, kernel, vpsp, comm_fft,pawtab,psps%usepaw)
231 
232 !----------------------------------------------------------------------
233 ! ----- Option 2: compute forces induced by local ionic potential -----
234 !----------------------------------------------------------------------
235    else if (option == 2) then
236 
237 !    Compute Hartree potential from rhor
238      ABI_ALLOCATE(vhartr,(nfft))
239      call psolver_hartree(entmp, (/ hgx, hgy, hgz /), icoulomb, me_fft, comm_fft, nfft, &
240 &     (/n1,n2,n3/), nproc_fft, nscforder, nspden, rhor, vhartr, usewvl)
241 
242 !    Allocate temporary array for forces
243      ABI_ALLOCATE(gxyz,(3, natom))
244 
245 !    Calculate local part of the forces grtn (inspired from BigDFT routine)
246      call local_forces_new(fftn3_distrib,ffti3_local,geocode,me_fft, ntypat, natom, &
247 &     typat, xcart, gridcart, psps%gth_params%psppar, int(psps%ziontypat), &
248 &     hgx,hgy,hgz, n1,n2,n3,n3d, rhor,vhartr, gxyz, pawtab,psps%usepaw)
249 
250 !    Forces should be in reduced coordinates.
251      do ia = 1, natom, 1
252        do igeo = 1, 3, 1
253          grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) &
254 &         - rprimd(2, igeo) * gxyz(2, ia) &
255 &         - rprimd(3, igeo) * gxyz(3, ia)
256        end do
257      end do
258 
259 !    Deallocate local variables
260      ABI_DEALLOCATE(vhartr)
261      ABI_DEALLOCATE(gxyz)
262    end if
263 
264 !----------------------------------------------------------------------
265 ! ----- Section for the non-HGH/GTH/PAW pseudopotentials (testing) ----
266 !----------------------------------------------------------------------
267  else
268 
269    if (testing) then
270      call system_clock(count_rate = countParSeconde)
271      call system_clock(tpsStart, count_rate = countParSeconde)
272    end if
273 
274 !  dr is the r step in the sampling psps%vlspl
275    dr = psps%qgrid_vl(2)
276    invdr = 1._dp / dr
277    dr2div6 = dr * dr / 6._dp
278 
279    if (option == 1) then
280 !    Set 0 in vpsp before summing
281      vpsp(:) = 0._dp
282    else if (option == 2) then
283 !    Allocate array to store cartesian gradient computed with
284 !    an interpolation of rhor
285      ABI_ALLOCATE(grtn_cart_interpol,(3, natom))
286      grtn_cart_interpol(:, :) = 0._dp
287 
288      n_interpol = nStep ** 3
289      ABI_ALLOCATE(coordRed_interpol,(3, nStep ** 3))
290      ABI_ALLOCATE(coordCart_interpol,(3, nStep ** 3))
291 
292      if (testing .and. customRho) then
293 !      Use a custom rho instead of the self-consistent one.
294        ABI_ALLOCATE(rhor_testing,(nfft))
295        ABI_ALLOCATE(rhog_testing,(2, nfft))
296      end if
297 
298      ABI_ALLOCATE(rhor_interpol,(nfft * n_interpol))
299      ABI_ALLOCATE(rhor_work,(nfft * n_interpol))
300      ABI_ALLOCATE(rhog_interpol,(2, nfft * n_interpol))
301 
302      if (testing .and. customRho) then
303 !      Testing only, changing rho with a centered gaussian
304        do ii = 1, nfft, 1
305 !        using the position of the first atom as center.
306          r = (gridcart(1, ii) - xcart(1, 1)) ** 2 + &
307 &         (gridcart(2, ii) - xcart(2, 1)) ** 2 + &
308 &         (gridcart(3, ii) - xcart(3, 1)) ** 2
309          rhor_testing(ii) = exp(-r/4._dp)
310        end do
311 !      Testing only, compute rhog_testing from rhor_testing
312        call fourdp(1,rhog_testing,rhor_testing,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
313      end if
314 
315 !    Compute the interpolation of rho, using a fourier transform
316      rhog_interpol(:, :) = 0._dp
317      ii = 0
318      do i3 = 1, n3, 1
319        if (i3 <= n3 / 2) then
320          iz = i3
321        else
322          iz = n3 * nStep - n3 + i3
323        end if
324        do i2 = 1, n2, 1
325          if (i2 <= n2 / 2) then
326            iy = i2
327          else
328            iy = n2 * nStep - n2 + i2
329          end if
330          do i1 = 1, n1, 1
331            ii = ii + 1
332            if (i1 <= n1 / 2) then
333              ix = i1
334            else
335              ix = n1 * nStep - n1 + i1
336            end if
337            jj = (iz - 1) * n2 * n1 * nStep ** 2 + (iy - 1) * n3 * nStep + ix
338            if (testing .and. customRho) then
339              rhog_interpol(:, jj) = rhog_testing(:, ii)
340            else
341              rhog_interpol(:, jj) = rhog(:, ii)
342            end if
343          end do
344        end do
345      end do
346 
347 !    Compute the interpolation of rho from the Fourier transformation
348      ngfft_interpol(:) = ngfft(:)
349      ngfft_interpol(1:3) = (/ n1 * nStep, n2 * nStep, n3 * nStep /)
350      ngfft_interpol(4:6) = (/ n1 * nStep + 1, n2 * nStep + 1, n3 * nStep /)
351      call fourdp(1,rhog_interpol,rhor_work,1,mpi_enreg,nfft*n_interpol,ngfft_interpol,mpi_enreg%paral_kgb,0)
352 
353 !    Reorder rhor_interpol to be able to read it linearly
354      jj = 0
355      do i3 = 1, n3, 1
356        do i2 = 1, n2, 1
357          do i1 = 1, n1, 1
358            do iz = 1, nStep, 1
359              do iy = 1, nStep, 1
360                do ix = 1, nStep, 1
361                  jj = jj + 1
362                  kk = ((i3 - 1) * nStep + iz - 1) ! z coordinate in the interpolated grid
363                  kk = kk * n1 * n2 * nStep ** 2
364                  kk = kk + ((i2 - 1) * nStep + iy - 1) * n1 * nStep ! adding y coordinate
365                  kk = kk + (i1 - 1) * nStep + ix ! adding x coordinate
366                  rhor_interpol(jj) = rhor_work(kk)
367                end do
368              end do
369            end do
370          end do
371        end do
372      end do
373      ABI_DEALLOCATE(rhor_work)
374 
375 !    Compute grid access in the interpolated volume
376      ii = 0
377      do iz = 1, nStep, 1
378        z = real(iz - 1, dp) / real(nStep, dp)
379        do iy = 1, nStep, 1
380          y = real(iy - 1, dp) / real(nStep, dp)
381          do ix = 1, nStep, 1
382            x = real(ix - 1, dp) / real(nStep, dp)
383            ii = ii + 1
384            coordRed_interpol(:, ii) = (/ x, y, z /)
385 !          Assuming orthogonal box (should be change later)
386            coordCart_interpol(:, ii) = (/ x * rprimd(1, 1) / real(n1, dp), &
387 &           y * rprimd(2, 2) / real(n2, dp), &
388 &           z * rprimd(3, 3) / real(n3, dp) /)
389          end do
390        end do
391      end do
392 
393      vol_interpol = 1._dp / real(nStep, dp) ** 3
394 !    Compute the coordinates (integer) of each atom and deduce
395 !    the max extens of the integral summation.
396 !    !!  do ia = 1, natom, 1
397 !    !!   coordAtom(1, ia) = int(xred(1, ia) * n1) + 1
398 !    !!   coordAtom(2, ia) = int(xred(2, ia) * n2) + 1
399 !    !!   coordAtom(3, ia) = int(xred(3, ia) * n3) + 1
400 !    !!  end do
401    end if
402 
403    if (testing .and. option == 2) then
404      call system_clock(tpsStop, count_rate = countParSeconde)
405      write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde)
406    end if
407 
408    ia1=1
409    do itypat = 1, ntypat, 1
410 !    ia1,ia2 sets range of loop over atoms:
411      ia2 = ia1 + nattyp(itypat) - 1
412 
413      do ii = 1, nfft, 1
414        do ia = ia1, ia2, 1
415          if (option == 1) then
416 !          Compute the potential
417 !          r is the distance between grid point and atom
418            r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + &
419 &           (gridcart(2, ii) - xcart(2, ia)) ** 2 + &
420 &           (gridcart(3, ii) - xcart(3, ia)) ** 2)
421 
422 !          Coefficients needed to compute the spline.
423            jj = int(r * invdr) + 1
424            if (jj > psps%mqgrid_vl - 2) then
425              write(message, '(3a,i0,a,i0,a,a)' )&
426 &             '  pseudo-potential local part sampling is not wide enough', ch10, &
427 &             '  want to access position ', jj, ' whereas mqgrid_vl = ', psps%mqgrid_vl, ch10, &
428 &             '  Action : no idea, contact developpers...'
429              MSG_ERROR(message)
430            end if
431            delta = r - psps%qgrid_vl(jj)
432            bb = delta * invdr
433            aa = 1._dp - bb
434            cc = aa * (aa ** 2 - 1._dp) * dr2div6
435            dd = bb * (bb ** 2 - 1._dp) * dr2div6
436 
437 !          compute V(r) from the spline, jj and jj + 1 is braketting r in
438 !          the sampling
439            deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + &
440 &           cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat)
441 !          Add on grid point ii the contribution of atom ia
442            vpsp(ii) = vpsp(ii) + deltaV
443          else if (option == 2) then
444 !          Compute the forces, as gradient of energy (V(r).rho(r))
445 
446 !          Testing only - reference points
447            if (.false.) then
448 !            r is the distance between grid point and atom
449              r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + &
450 &             (gridcart(2, ii) - xcart(2, ia)) ** 2 + &
451 &             (gridcart(3, ii) - xcart(3, ia)) ** 2)
452 
453 !            Coefficients needed to compute the spline.
454              jj = int(r * invdr) + 1
455              delta = r - psps%qgrid_vl(jj)
456              bb = delta * invdr
457              aa = 1._dp - bb
458              cc = aa * (aa ** 2 - 1._dp) * dr2div6
459              dd = bb * (bb ** 2 - 1._dp) * dr2div6
460 
461 !            When mesh position is on a node, forces are null.
462              if (r /= 0._dp) then
463 !              This value deltaV is the first derivative of V(r) taken at r.
464                deltaV = aa * psps%dvlspl(jj, 1, itypat) + bb * psps%dvlspl(jj + 1, 1, itypat) + &
465 &               cc * psps%dvlspl(jj, 2, itypat) + dd * psps%dvlspl(jj + 1, 2, itypat)
466 !              We multiply by rho(r) to have an energy.
467                deltaV = deltaV * rhor(ii, 1) / r
468                refValue(:) = - deltaV * (gridcart(:, ii) - xcart(:, ia))
469                grtn_cart_interpol(:, ia) = grtn_cart_interpol(:, ia) + refValue(:)
470              end if
471            end if
472 
473 !          Compute the interpolation for the point ii
474            ind = (ii - 1) * n_interpol
475            do kk = 1, n_interpol, 1
476              ind = ind + 1
477 
478              if (rhor_interpol(ind) > min_rho_value) then
479 !              Assume orthogonal box...
480                coordXYZ(1) = gridcart(1, ii) - xcart(1, ia) + coordCart_interpol(1, kk)
481                coordXYZ(2) = gridcart(2, ii) - xcart(2, ia) + coordCart_interpol(2, kk)
482                coordXYZ(3) = gridcart(3, ii) - xcart(3, ia) + coordCart_interpol(3, kk)
483                r = coordXYZ(1) ** 2 + coordXYZ(2) ** 2 + coordXYZ(3) ** 2
484 
485                if (r /= 0._dp) then
486                  r = sqrt(r)
487 !                Coefficients needed to compute the spline.
488                  jj = int(r * invdr) + 1
489                  delta = r - psps%qgrid_vl(jj)
490                  bb = delta * invdr
491                  aa = 1._dp - bb
492                  cc = aa * (aa ** 2 - 1._dp) * dr2div6
493                  dd = bb * (bb ** 2 - 1._dp) * dr2div6
494                  deltaV = aa * psps%dvlspl(jj, 1, itypat) + &
495 &                 bb * psps%dvlspl(jj + 1, 1, itypat) + &
496 &                 cc * psps%dvlspl(jj, 2, itypat) + &
497 &                 dd * psps%dvlspl(jj + 1, 2, itypat)
498                  deltaV = deltaV * rhor_interpol(ind) / r
499                  grtn_cart_interpol(1, ia) = grtn_cart_interpol(1, ia) - deltaV * coordXYZ(1)
500                  grtn_cart_interpol(2, ia) = grtn_cart_interpol(2, ia) - deltaV * coordXYZ(2)
501                  grtn_cart_interpol(3, ia) = grtn_cart_interpol(3, ia) - deltaV * coordXYZ(3)
502 !                do igeo = 1, 3, 1
503 !                grtn_cart_interpol(igeo, ia) = grtn_cart_interpol(igeo, ia) - deltaV * coordXYZ(igeo)
504 !                end do
505                end if
506              end if
507            end do
508 
509 !          =============
510 !          Testing only
511 !          =============
512 !          use of finite differences
513            if (finiteDiff) then
514              do igeo = 1, 3, 1
515                coord(:) = 0._dp
516                coord(igeo) = dr / 2.0_dp
517                r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + &
518 &               (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + &
519 &               (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2)
520 
521 !              Coefficients needed to compute the spline.
522                jj = int(r * invdr) + 1
523                delta = r - psps%qgrid_vl(jj)
524                bb = delta * invdr
525                aa = 1._dp - bb
526                cc = aa * (aa ** 2 - 1._dp) * dr2div6
527                dd = bb * (bb ** 2 - 1._dp) * dr2div6
528 
529                deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + &
530 &               cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat)
531 
532 
533                coord(:) = 0._dp
534                coord(igeo) = -dr / 2.0_dp
535                r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + &
536 &               (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + &
537 &               (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2)
538 
539 !              Coefficients needed to compute the spline.
540                jj = int(r * invdr) + 1
541                delta = r - psps%qgrid_vl(jj)
542                bb = delta * invdr
543                aa = 1._dp - bb
544                cc = aa * (aa ** 2 - 1._dp) * dr2div6
545                dd = bb * (bb ** 2 - 1._dp) * dr2div6
546 
547                deltaV = deltaV - (aa * psps%vlspl(jj, 1, itypat) + &
548 &               bb * psps%vlspl(jj + 1, 1, itypat) + &
549 &               cc * psps%vlspl(jj, 2, itypat) + &
550 &               dd * psps%vlspl(jj + 1, 2, itypat))
551                grtn_diff(igeo, ia) = grtn_diff(igeo, ia) - deltaV * rhor(ii, 1) / dr
552              end do
553            end if
554 !          =============
555 !          Testing only
556 !          =============
557 
558          end if
559        end do
560 !      End loop over atoms of type itypat
561      end do
562 !    End loop over real space grid points
563 
564      ia1 = ia2 + 1
565    end do
566 !  End loop over type of atoms
567 
568    if(option==2)then
569 !    multiply the forces by the volume of a single box mesh.
570      grtn_cart_interpol(:, :) = grtn_cart_interpol(:, :) * &
571 &     ucvol / real(n1 * n2 * n3, dp) * vol_interpol
572 !    Transform cartesian forces to reduce coordinates
573      do ia = 1, natom, 1
574        do igeo = 1, 3, 1
575          grtn(igeo, ia) = rprimd(1, igeo) * grtn_cart_interpol(1, ia) + &
576 &         rprimd(2, igeo) * grtn_cart_interpol(2, ia) + &
577 &         rprimd(3, igeo) * grtn_cart_interpol(3, ia)
578        end do
579      end do
580      ABI_DEALLOCATE(rhor_interpol)
581      ABI_DEALLOCATE(rhog_interpol)
582      ABI_DEALLOCATE(coordRed_interpol)
583      ABI_DEALLOCATE(coordCart_interpol)
584      if (testing .and. customRho) then
585        ABI_DEALLOCATE(rhor_testing)
586        ABI_DEALLOCATE(rhog_testing)
587      end if
588 
589      if (testing) then
590        call system_clock(tpsStop, count_rate = countParSeconde)
591        write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde)
592        write(std_out,*) grtn_cart_interpol
593        MSG_ERROR("Testing section!")
594      end if
595 
596    end if
597 
598 !-------------------------
599  end if ! GTH/HGH/PAW psps
600 
601 !Release temporary memory
602  ABI_DEALLOCATE(xcart)
603  ABI_DEALLOCATE(gridcart)
604 
605 !Close timing counters
606  if (option==2)then
607    call timab(72,2,tsec)
608  end if
609 
610 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=informations 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.

PARENTS

      mklocl,wvl_wfsinp_scratch

CHILDREN

      calcdvloc_wvl,derf_ab,paw_splint_der

SOURCE

1447 subroutine mklocl_wavelets(efield, grtn, mpi_enreg, natom, nfft, &
1448      & nspden, option, rprimd, vpsp, wvl_den, wvl_descr, xcart)
1449 
1450  use defs_wvltypes
1451 #if defined HAVE_BIGDFT
1452  use BigDFT_API, only : ELECTRONIC_DENSITY,createIonicPotential,local_forces
1453  use poisson_solver, only : H_potential
1454 #endif
1455 
1456 !This section has been created automatically by the script Abilint (TD).
1457 !Do not modify the following lines by hand.
1458 #undef ABI_FUNC
1459 #define ABI_FUNC 'mklocl_wavelets'
1460 !End of the abilint section
1461 
1462  implicit none
1463 
1464 !Arguments ------------------------------------
1465 !scalars
1466  integer,intent(in) :: option, natom, nfft, nspden
1467  type(MPI_type),intent(in) :: mpi_enreg
1468  type(wvl_denspot_type), intent(inout) :: wvl_den
1469  type(wvl_internal_type), intent(in) :: wvl_descr
1470 !arrays
1471  real(dp),intent(in) :: rprimd(3,3),efield(3)
1472  real(dp),intent(inout) :: grtn(3,natom)
1473  real(dp), intent(inout) :: vpsp(nfft)
1474  real(dp),intent(inout) :: xcart(3,natom)
1475 
1476 !Local variables-------------------------------
1477 #if defined HAVE_BIGDFT
1478 !scalars
1479  integer :: i,i1,i2,i3,ia,ierr,igeo,me,nproc,shift,comm
1480  real(dp) :: energ
1481  character(len=500) :: message
1482 !arrays
1483  real(dp) :: epot(3)
1484  real(dp) :: elecfield(3)=(/zero,zero,zero/) ! Not used here
1485  real(dp),allocatable :: gxyz(:,:),vhartr(:),rhov(:,:)
1486 #endif
1487 
1488 ! *********************************************************************
1489 
1490 #if defined HAVE_BIGDFT
1491 
1492  elecfield=zero !not used here
1493 
1494 !Manage parallelism
1495  comm=mpi_enreg%comm_wvl
1496  nproc=xmpi_comm_size(comm)
1497  me=xmpi_comm_rank(comm)
1498 
1499 !----------------------------------------------------------------------
1500 ! ----- Option 1: compute local ionic potential                   -----
1501 !----------------------------------------------------------------------
1502  if (option == 1) then
1503 
1504 !  We get the kernel for the Poisson solver (used to go from the ionic
1505 !  charge to the potential).If the kernel is uncomputed, it does it now.
1506 !  call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 2, icoulomb, me, kernel, &
1507 !&                     comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder)
1508 !  if (.not.associated(kernel%co%kernel)) then
1509 !    call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 1, icoulomb, me, kernel, &
1510 !&                       comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder)
1511 !  end if
1512 
1513    message=ch10//' mklocl_wavelets: Create local potential from ions.'
1514    call wrtout(std_out,message,'COLL')
1515 
1516    shift = 1 + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1517 &   * wvl_den%denspot%dpbox%i3xcsh
1518 
1519 !  Call the BigDFT routine
1520    call createIonicPotential(wvl_descr%atoms%astruct%geocode, me, nproc, (me == 0), wvl_descr%atoms, &
1521 &   xcart, wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), &
1522 &   wvl_den%denspot%dpbox%hgrids(3), &
1523 &   elecfield, wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, &
1524 &   wvl_den%denspot%dpbox%n3pi, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, &
1525 &   wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), &
1526 &   wvl_den%denspot%dpbox%ndims(3), wvl_den%denspot%pkernel, vpsp(shift), 0.d0,wvl_descr%rholoc)
1527 
1528 !  Copy vpsp into bigdft object:
1529    call wvl_rhov_abi2big(1,vpsp,wvl_den%denspot%v_ext,shift=(shift-1))
1530 
1531 !  Eventually add the electric field
1532    if (maxval(efield) > tol12) then
1533      message=ch10//'mklocl_wavelets: Add the electric field.'
1534      call wrtout(std_out,message,'COLL')
1535 !    We add here the electric field since in BigDFT, the field must be on x...
1536      epot(:) = real(0.5, dp) * efield(:) * wvl_den%denspot%dpbox%hgrids(:)
1537      do i3 = 1, wvl_den%denspot%dpbox%n3pi, 1
1538        ia = (i3 - 1) * wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2)
1539        do i2 = -14, 2 * wvl_descr%Glr%d%n2 + 16, 1
1540          i = ia + (i2 + 14) * wvl_den%denspot%dpbox%ndims(1)
1541          do i1 = -14, 2 * wvl_descr%Glr%d%n1 + 16, 1
1542            i = i + 1
1543            vpsp(shift + i) = vpsp(shift + i) + &
1544 &           epot(1) * real(i1 - wvl_descr%Glr%d%n1, dp) + &
1545 &           epot(2) * real(i2 - wvl_descr%Glr%d%n2, dp) + &
1546 &           epot(3) * real(i3 - wvl_descr%Glr%d%n3, dp)
1547          end do
1548        end do
1549      end do
1550    end if
1551 
1552 !----------------------------------------------------------------------
1553 ! ----- Option 2: compute forces induced by local ionic potential -----
1554 !----------------------------------------------------------------------
1555  else if (option == 2) then
1556 
1557    message=ch10//' mklocl_wavelets: compute local forces.'
1558    call wrtout(std_out,message,'COLL')
1559 
1560    if (wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then
1561      message='denspot bigdft datstructure should contain rhor!'
1562      MSG_BUG(message)
1563    end if
1564 
1565 !  Extract density rhor from bigDFT datastructure
1566    ABI_ALLOCATE(rhov,(nfft, nspden))
1567    ABI_ALLOCATE(vhartr,(nfft))
1568    shift = wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1569 &   * wvl_den%denspot%dpbox%i3xcsh
1570    do i = 1, nfft
1571      rhov(i, 1) = wvl_den%denspot%rhov(i + shift)
1572      vhartr(i)  = wvl_den%denspot%rhov(i + shift)
1573    end do
1574    if (nspden == 2) then
1575      shift = shift + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) &
1576 &     * wvl_den%denspot%dpbox%n3d
1577      do i = 1, nfft
1578        rhov(i, 2) =             wvl_den%denspot%rhov(i + shift)
1579        vhartr(i)  = vhartr(i) + wvl_den%denspot%rhov(i + shift)
1580      end do
1581    end if
1582 
1583 !  Compute Hartree potential from rhor
1584    call H_potential('D',wvl_den%denspot%pkernel,vhartr,vhartr,energ,zero,.false.)
1585 
1586 !  Allocate temporary array for forces
1587    ABI_ALLOCATE(gxyz,(3, natom))
1588 
1589 !  Calculate local part of the forces grtn (modified BigDFT routine)
1590    call local_forces_wvl(me,natom,xcart,&
1591 &   wvl_den%denspot%dpbox%hgrids(1),&
1592 &   wvl_den%denspot%dpbox%hgrids(2),&
1593 &   wvl_den%denspot%dpbox%hgrids(3),&
1594 &   wvl_descr%Glr%d%n1,wvl_descr%Glr%d%n2,wvl_descr%Glr%d%n3,&
1595 &   wvl_den%denspot%dpbox%n3p,&
1596 &   wvl_den%denspot%dpbox%i3s+wvl_den%denspot%dpbox%i3xcsh,&
1597 &   wvl_den%denspot%dpbox%ndims(1),wvl_den%denspot%dpbox%ndims(2),&
1598 &   rhov,vhartr,gxyz,wvl_descr)
1599 !    call local_forces(me, wvl_descr%atoms, xcart, &
1600 ! &   wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), &
1601 ! &   wvl_den%denspot%dpbox%hgrids(3), &
1602 ! &   wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, &
1603 ! &   wvl_den%denspot%dpbox%n3p, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, &
1604 ! &   wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), &
1605 ! &   rhov, vhartr,gxyz,locstrten,charge)
1606 
1607 !    Pending: floc,locstrten and charge are not used here.
1608 !    Pending: put mpi_enreg%nscatterarr... in object denspot, initialize object, etc.
1609 
1610    if (nproc > 1) then
1611      call xmpi_sum(gxyz, comm, ierr)
1612    end if
1613 
1614 !  Forces should be in reduced coordinates.
1615    do ia = 1, natom, 1
1616      do igeo = 1, 3, 1
1617        grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) &
1618 &       - rprimd(2, igeo) * gxyz(2, ia) &
1619 &       - rprimd(3, igeo) * gxyz(3, ia)
1620      end do
1621    end do
1622 
1623 !  Deallocate local variables
1624    ABI_DEALLOCATE(vhartr)
1625    ABI_DEALLOCATE(rhov)
1626    ABI_DEALLOCATE(gxyz)
1627 
1628 !----------------------------------------------------------------------
1629 
1630  else ! option switch
1631    message = 'Internal error, option should be 1 or 2!'
1632    MSG_ERROR(message)
1633  end if
1634 
1635 #else
1636  BIGDFT_NOTENABLED_ERROR()
1637  if (.false.) write(std_out,*) option,natom,nfft,nspden,mpi_enreg%me,&
1638 & wvl_den%symObj,wvl_descr%h(1),rprimd(1,1),efield(1),grtn(1,1),vpsp(1),xcart(1,1)
1639 #endif
1640 
1641 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

PARENTS

      mklocl_realspace

CHILDREN

SOURCE

1248 subroutine calcdVloc_mklocl(yy,xx,rloc,Z)
1249 
1250 
1251 !This section has been created automatically by the script Abilint (TD).
1252 !Do not modify the following lines by hand.
1253 #undef ABI_FUNC
1254 #define ABI_FUNC 'calcdVloc_mklocl'
1255 !End of the abilint section
1256 
1257  implicit none
1258 !Arguments ------------------------------------
1259 !scalars
1260  real(dp),intent(in)  :: xx,rloc,Z
1261  real(dp),intent(out) :: yy
1262 
1263 !Local variables-------------------------------
1264  !scalars
1265  real(dp):: arg,tt
1266 
1267 ! *************************************************************************
1268 
1269      arg=xx/(sqrt(2._dp)*rloc)
1270      call derf_ab(tt,arg)
1271      yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1272 
1273    end subroutine calcdVloc_mklocl

mklocl_realspace/calcVloc_mklocl [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  calcVloc_mklocl

FUNCTION

INPUTS

OUTPUT

PARENTS

      mklocl_realspace

CHILDREN

SOURCE

914 subroutine calcVloc_mklocl(yy,xx,rloc,Z)
915 
916 
917 !This section has been created automatically by the script Abilint (TD).
918 !Do not modify the following lines by hand.
919 #undef ABI_FUNC
920 #define ABI_FUNC 'calcVloc_mklocl'
921 !End of the abilint section
922 
923  implicit none
924 !Arguments ------------------------------------
925 !scalars
926  real(dp),intent(in)  :: xx,rloc,Z
927  real(dp),intent(out) :: yy
928 
929 !Local variables-------------------------------
930  !scalars
931  real(dp):: arg,tt
932 
933 ! *************************************************************************
934 
935    arg=xx/(sqrt(2.0)*rloc)
936    call derf_ab(tt,arg)
937    yy=-Z/xx*tt
938 
939  end subroutine calcVloc_mklocl

mklocl_realspace/createIonicPotential_new [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  createIonicPotential_new

FUNCTION

INPUTS

OUTPUT

PARENTS

      mklocl_realspace

CHILDREN

SOURCE

632 subroutine createIonicPotential_new(fftn3_distrib,ffti3_local,geocode,iproc,&
633 &  nproc,nat,ntypes,iatype,psppar,nelpsp,rxyz,gridcart,&
634 &  hxh,hyh,hzh,n1i,n2i,n3d,n3i,kernel,pot_ion,spaceworld,pawtab,usepaw)
635 
636  use defs_wvltypes, only : coulomb_operator
637 
638 !This section has been created automatically by the script Abilint (TD).
639 !Do not modify the following lines by hand.
640 #undef ABI_FUNC
641 #define ABI_FUNC 'createIonicPotential_new'
642 !End of the abilint section
643 
644 implicit none
645 
646 !Arguments -------------------------------
647 !scalars
648  integer, intent(in) :: iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw
649  real(dp), intent(in) :: hxh,hyh,hzh
650  character(len=1), intent(in) :: geocode
651  type(coulomb_operator), intent(in) :: kernel
652 !arrays
653  integer, dimension(nat), intent(in) :: iatype
654  integer, dimension(ntypes), intent(in) :: nelpsp
655  integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local
656  real(dp), dimension(3,n1i*n2i*n3d), intent(in) :: gridcart
657  real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar
658  real(dp), dimension(3,nat), intent(in) :: rxyz
659  real(dp), dimension(*), intent(inout) :: pot_ion
660  type(pawtab_type),intent(in)  :: pawtab(ntypes*usepaw)
661 
662 !Local variables -------------------------
663 #if defined HAVE_BIGDFT
664 !scalars
665  integer :: iat,i1,i2,i3,j1,j2,j3,isx,isy,isz,iex,iey,iez,ierr,ityp
666  integer :: ind,nloc,iloc,i3loc,msz
667  logical :: gox,goy,goz,perx,pery,perz
668  real(dp) :: arg,charge,cutoff,ehart,eexcu,rholeaked,rholeaked_tot,rloc
669  real(dp) :: rx,ry,rz,rzero,r2,tt,tt_tot,vexcu,vhgh,x,xp,y,z
670 !arrays
671  real(dp) :: rr(1),vpaw(1)
672  real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
673 #endif
674 
675 ! *********************************************************************
676 
677 #if defined HAVE_BIGDFT
678 
679  if(nproc<0)then
680    MSG_ERROR('nproc should not be negative')
681  end if
682 
683 !Ionic charge (must be calculated for the PS active processes)
684  rholeaked=0._dp
685 !Ionic energy (can be calculated for all the processors)
686 
687 !here we should insert the calculation of the ewald energy for the periodic BC case
688 !!!  eion=0._dp
689 !!!  do iat=1,nat
690 !!!     ityp=iatype(iat)
691 !!!     rx=rxyz(1,iat)
692 !!!     ry=rxyz(2,iat)
693 !!!     rz=rxyz(3,iat)
694 !!!     !    ion-ion interaction
695 !!!     do jat=1,iat-1
696 !!!        dist=sqrt( (rx-rxyz(1,jat))**2+(ry-rxyz(2,jat))**2+(rz-rxyz(3,jat))**2 )
697 !!!        jtyp=iatype(jat)
698 !!!        eion=eion+real(nelpsp(jtyp)*nelpsp(ityp),kind=dp)/dist
699 !!!     enddo
700 !!!  end do
701 !!!  if (iproc.eq.0) write(std_out,'(1x,a,1pe22.14)') 'ion-ion interaction energy',eion
702 
703 !Creates charge density arising from the ionic PSP cores
704 !the n3pi dimension indicates the number of planes trated by each processor in the FFT parallelisation
705 !for a plane wave treatment this value depends on whether the direct space is divided in planes or not
706 !I don't know this variable, which in the future must be inserted at the place of n3pi (LG)
707 !if n3pi=0 this means that the processors doesn't calculate anything
708 !if (n3pi >0 ) then
709 
710 !conditions for periodicity in the three directions
711  perx=(geocode /= 'F')
712  pery=(geocode == 'P')
713  perz=(geocode /= 'F')
714 
715 !this initialise the array to zero, it will work only if bigdft library is enabled
716  pot_ion(1:n1i*n2i*n3d)=zero
717 
718  do iat=1,nat
719    ityp=iatype(iat)
720    rx=rxyz(1,iat)
721    ry=rxyz(2,iat)
722    rz=rxyz(3,iat)
723 
724    rloc=psppar(0,0,ityp)
725    charge=real(nelpsp(ityp),kind=dp)/(2._dp*pi*sqrt(2._dp*pi)*rloc**3)
726    cutoff=10._dp*rloc
727 
728    isx=floor((rx-cutoff)/hxh)
729    isy=floor((ry-cutoff)/hyh)
730    isz=floor((rz-cutoff)/hzh)
731 
732    iex=ceiling((rx+cutoff)/hxh)
733    iey=ceiling((ry+cutoff)/hyh)
734    iez=ceiling((rz+cutoff)/hzh)
735 
736 !  Calculate Ionic Density
737 !  using HGH parameters.
738 !  Eq. 1.104, T. Deutsch and L. Genovese, JDN. 12, 2011
739    do i3=isz,iez
740      z=real(i3,kind=dp)*hzh-rz
741      call ind_positions_mklocl(perz,i3,n3i,j3,goz)
742      if(fftn3_distrib(j3)==iproc) then
743        i3loc=ffti3_local(j3)
744        do i2=isy,iey
745          y=real(i2,kind=dp)*hyh-ry
746          call ind_positions_mklocl(pery,i2,n2i,j2,goy)
747          do i1=isx,iex
748            x=real(i1,kind=dp)*hxh-rx
749            call ind_positions_mklocl(perx,i1,n1i,j1,gox)
750            r2=x**2+y**2+z**2
751            if (goz  .and. goy  .and. gox ) then
752              ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i
753              r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2
754            end if
755            arg=r2/rloc**2
756            xp=exp(-.5d0*arg)
757            if (goz  .and. goy  .and. gox ) then
758              pot_ion(ind)=pot_ion(ind)-xp*charge
759            else
760              rholeaked=rholeaked+xp*charge
761            end if
762          end do
763        end do
764      end if
765    end do
766 
767  end do
768 
769 !Check
770  tt=0._dp
771  do j3= 1,n3d
772    do i2= 1,n2i
773      do i1= 1,n1i
774        ind=i1+(i2-1)*n1i+(j3-1)*n1i*n2i
775        tt=tt+pot_ion(ind)
776      end do
777    end do
778  end do
779 
780  tt=tt*hxh*hyh*hzh
781  rholeaked=rholeaked*hxh*hyh*hzh
782 
783  call xmpi_sum(tt,tt_tot,spaceworld,ierr)
784  call xmpi_sum(rholeaked,rholeaked_tot,spaceworld,ierr)
785 
786  if (iproc.eq.0) then
787    write(std_out,'(1x,a,f26.12,2x,1pe10.3)') &
788 &   'total ionic charge, leaked charge ',tt_tot,rholeaked_tot
789  end if
790 
791 !Here the value of the datacode must be kept fixed
792 !there can be some problems when running this stuff in parallel,
793 !  if the ionic potential distribution does not agree with the
794 !  plane distribution which is supposed to hold for the Poisson Solver
795  call psolver(geocode,'D',iproc,nproc,n1i,n2i,n3i,0,hxh,hyh,hzh,&
796 & pot_ion,kernel%kernel,pot_ion,ehart,eexcu,vexcu,0._dp,.false.,1)
797 
798 !Add the remaining short-range local terms
799  do iat=1,nat
800    ityp=iatype(iat)
801 
802    rx=rxyz(1,iat)
803    ry=rxyz(2,iat)
804    rz=rxyz(3,iat)
805 
806 !  determine number of local terms
807    rloc=psppar(0,0,ityp)
808    cutoff=10._dp*rloc
809    charge=real(nelpsp(ityp),kind=dp)
810 
811 !  determine number of local terms (HGH pot)
812    nloc=0
813    do iloc=1,4
814      if (psppar(0,iloc,ityp).ne.0._dp) nloc=iloc
815    end do
816 
817 !  PAW specifics
818    if (usepaw==1) then
819      msz=pawtab(ityp)%wvl%rholoc%msz
820      rad    => pawtab(ityp)%wvl%rholoc%rad(1:msz)
821      vloc   => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
822      d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
823      rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2)
824    end if
825 
826    isx=floor((rx-cutoff)/hxh)
827    isy=floor((ry-cutoff)/hyh)
828    isz=floor((rz-cutoff)/hzh)
829 
830    iex=ceiling((rx+cutoff)/hxh)
831    iey=ceiling((ry+cutoff)/hyh)
832    iez=ceiling((rz+cutoff)/hzh)
833 
834    do i3=isz,iez
835      z=real(i3,kind=dp)*hzh-rz
836      call ind_positions_mklocl(perz,i3,n3i,j3,goz)
837      if(fftn3_distrib(j3) == iproc .and. goz) then !MPI
838        i3loc=ffti3_local(j3)
839        if (goz) then
840          do i2=isy,iey
841            y=real(i2,kind=dp)*hyh-ry
842            call ind_positions_mklocl(pery,i2,n2i,j2,goy)
843            if (goy) then
844              do i1=isx,iex
845                x=real(i1,kind=dp)*hxh-rx
846                call ind_positions_mklocl(perx,i1,n1i,j1,gox)
847                if (gox) then
848                  ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i
849                  r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2
850 !                r2=x**2+y**2+z**2
851 
852 !                HGH: V_S=gaussian potential of Eq. (9) in JCP 129, 014109(2008)
853                  if (usepaw==0) then
854                    if (nloc /= 0) then
855                      arg=r2/rloc**2
856                      xp=exp(-.5d0*arg)
857                      tt=psppar(0,nloc,ityp)
858                      do iloc=nloc-1,1,-1
859                        tt=arg*tt+psppar(0,iloc,ityp)
860                      end do
861                      pot_ion(ind)=pot_ion(ind)+xp*tt
862                    end if
863 
864 !                PAW: V_PAW-V_L^HGH
865                  else
866                    rr(1)=sqrt(r2)
867                    if (rr(1)>=rzero) then
868                      call paw_splint(msz,rad,vloc,d2vloc,1,rr,vpaw)
869                      call calcVloc_mklocl(vhgh,rr(1),rloc,charge)
870                      pot_ion(ind)=pot_ion(ind)+vpaw(1)-vhgh
871                    else
872                      pot_ion(ind)=pot_ion(ind)+vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
873                    end if
874                  end if
875 
876                end if
877              end do
878            end if
879          end do
880        end if
881      end if
882    end do
883 
884  end do  !iat
885 #else
886  BIGDFT_NOTENABLED_ERROR()
887  if (.false.) write(std_out,*) geocode,iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw,&
888 & hxh,hyh,hzh,iatype(1),nelpsp(1),fftn3_distrib(1),ffti3_local(1),gridcart(1,1),psppar(1,1,1),&
889 & rxyz(1,1),pot_ion(1),pawtab(1)%mesh_size,kernel%co
890 #endif
891 
892  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

PARENTS

      mklocl_realspace

CHILDREN

SOURCE

1375 subroutine ind_positions_mklocl(periodic,i,n,j,go)
1376 
1377 
1378 !This section has been created automatically by the script Abilint (TD).
1379 !Do not modify the following lines by hand.
1380 #undef ABI_FUNC
1381 #define ABI_FUNC 'ind_positions_mklocl'
1382 !End of the abilint section
1383 
1384  implicit none
1385 
1386 !Arguments -------------------------------
1387  logical, intent(in) :: periodic
1388  integer, intent(in) :: i,n
1389  logical, intent(out) :: go
1390  integer, intent(out) :: j
1391 
1392 ! *********************************************************************
1393 
1394      if (periodic) then
1395        go=.true.
1396        j=modulo(i-1,n)+1
1397      else
1398        j=i
1399        go=(i >= 1 .and. i <= n)
1400      end if
1401 
1402    end subroutine ind_positions_mklocl

mklocl_realspace/local_forces_new [ Functions ]

[ Top ] [ mklocl_realspace ] [ Functions ]

NAME

  local_forces_new

FUNCTION

INPUTS

OUTPUT

PARENTS

      mklocl_realspace

CHILDREN

SOURCE

1039 subroutine local_forces_new(fftn3_distrib,ffti3_local,&
1040      geocode,iproc,ntypes,nat,iatype,rxyz,gridcart,psppar,nelpsp,hxh,hyh,hzh,&
1041      n1,n2,n3,n3d,rho,pot,floc,pawtab,usepaw)
1042 
1043 
1044 !This section has been created automatically by the script Abilint (TD).
1045 !Do not modify the following lines by hand.
1046 #undef ABI_FUNC
1047 #define ABI_FUNC 'local_forces_new'
1048 !End of the abilint section
1049 
1050   implicit none
1051 
1052 !Arguments -------------------------------
1053 !scalars
1054  integer, intent(in) :: iproc,ntypes,nat,n1,n2,n3,n3d,usepaw
1055  character(len=1), intent(in) :: geocode
1056  real(dp), intent(in) :: hxh,hyh,hzh
1057 !arrays
1058  integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local
1059  integer, dimension(nat), intent(in) :: iatype
1060  integer, dimension(ntypes), intent(in) :: nelpsp
1061  real(dp), dimension(3,n1*n2*n3d), intent(in) :: gridcart
1062  real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar
1063  real(dp), dimension(3,nat), intent(in) :: rxyz
1064  real(dp), dimension(*), intent(in) :: rho,pot
1065  real(dp), dimension(3,nat), intent(out) :: floc
1066  type(pawtab_type),intent(in)  :: pawtab(ntypes*usepaw)
1067 
1068 !Local variables -------------------------
1069 !scalars
1070  integer :: isx,isy,isz,iex,iey,iez,i1,i2,i3,j1,j2,j3,ind,iat,ityp,iloc,i3loc,msz,nloc
1071  logical :: perx,pery,perz,gox,goy,goz
1072  real(dp) :: arg,charge,cutoff,dvhgh,fxerf,fyerf,fzerf,fxgau,fygau,fzgau,forceleaked
1073  real(dp) :: forceloc,prefactor,rloc,rloc2,rhoel,rx,ry,rz,rzero,r2,tt,x,xp,y,z,Vel
1074  real(dp), dimension(4) :: cprime
1075 !arrays
1076  real(dp) :: dvpawdr(1),rr(1)
1077  real(dp),pointer :: rad(:),vloc(:),d2vloc(:)
1078 
1079 ! *********************************************************************
1080 
1081    if (iproc == 0) write(std_out,'(1x,a)',advance='no')'Calculate local forces...'
1082 
1083 !Conditions for periodicity in the three directions
1084    perx=(geocode /= 'F')
1085    pery=(geocode == 'P')
1086    perz=(geocode /= 'F')
1087 
1088    forceleaked=zero
1089 
1090    do iat=1,nat
1091      ityp=iatype(iat)
1092 !  Coordinates of the center
1093      rx=rxyz(1,iat)
1094      ry=rxyz(2,iat)
1095      rz=rxyz(3,iat)
1096 
1097 !  Initialization of the forces
1098 !  ion-electron term, error function part
1099      fxerf=zero
1100      fyerf=zero
1101      fzerf=zero
1102 !  ion-electron term, gaussian part
1103      fxgau=zero
1104      fygau=zero
1105      fzgau=zero
1106 
1107 !  Building array of coefficients of the derivative of the gaussian part
1108      cprime(1)=2._dp*psppar(0,2,ityp)-psppar(0,1,ityp)
1109      cprime(2)=4._dp*psppar(0,3,ityp)-psppar(0,2,ityp)
1110      cprime(3)=6._dp*psppar(0,4,ityp)-psppar(0,3,ityp)
1111      cprime(4)=-psppar(0,4,ityp)
1112 
1113 !  Determine number of local terms (HGH pot)
1114      nloc=0
1115      do iloc=1,4
1116        if (psppar(0,iloc,ityp).ne.zero) nloc=iloc
1117      end do
1118 
1119 !  Some constants depending on the atom type
1120      rloc=psppar(0,0,ityp) ; rloc2=rloc**2
1121      charge=real(nelpsp(ityp),kind=dp)
1122      prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1123 
1124 !  PAW specifics
1125      if (usepaw==1) then
1126        msz=pawtab(ityp)%wvl%rholoc%msz
1127        rad    => pawtab(ityp)%wvl%rholoc%rad(1:msz)
1128        vloc   => pawtab(ityp)%wvl%rholoc%d(1:msz,3)
1129        d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4)
1130        rzero=rad(1);if (rad(1)<=1.d-10) rzero=rad(2)
1131      end if
1132 
1133 !  Maximum extension of the gaussian
1134      cutoff=10._dp*rloc
1135      isx=floor((rx-cutoff)/hxh)
1136      isy=floor((ry-cutoff)/hyh)
1137      isz=floor((rz-cutoff)/hzh)
1138      iex=ceiling((rx+cutoff)/hxh)
1139      iey=ceiling((ry+cutoff)/hyh)
1140      iez=ceiling((rz+cutoff)/hzh)
1141 
1142 !  Calculate the forces near the atom due to the gaussian
1143 !  and error function parts of the potential
1144      do i3=isz,iez
1145        z=real(i3,kind=dp)*hzh-rz
1146        call ind_positions_mklocl(perz,i3,n3,j3,goz)
1147        if(fftn3_distrib(j3)==iproc) then
1148          i3loc=ffti3_local(j3)
1149          do i2=isy,iey
1150            y=real(i2,kind=dp)*hyh-ry
1151            call ind_positions_mklocl(pery,i2,n2,j2,goy)
1152            do i1=isx,iex
1153              x=real(i1,kind=dp)*hxh-rx
1154              call ind_positions_mklocl(perx,i1,n1,j1,gox)
1155 
1156              if (goz.and.goy.and.gox) then
1157                ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2
1158                x=(gridcart(1,ind)-rx)
1159                y=(gridcart(2,ind)-ry)
1160                z=(gridcart(3,ind)-rz)
1161                r2=x**2+y**2+z**2
1162                xp=exp(-0.5_dp*r2/rloc**2)
1163 
1164 !            Short range part
1165                rhoel=rho(ind)
1166 !            HGH: V_S^prime=gaussian
1167                if (usepaw==0) then
1168                  if (nloc/=0) then
1169                    arg=r2/rloc**2
1170                    tt=cprime(nloc)
1171                    do iloc=nloc-1,1,-1
1172                      tt=arg*tt+cprime(iloc)
1173                    end do
1174                    forceloc=xp*tt*rhoel
1175                  else
1176                    forceloc=zero
1177                  end if
1178 !            PAW: V_PAW^prime-V_L^prime
1179                else
1180                  rr(1)=sqrt(r2)
1181                  if (rr(1)>=rzero) then
1182                    call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
1183                    call calcdVloc_mklocl(dvhgh,rr(1),rloc,charge)
1184                    forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1185                  else
1186                    forceloc=rhoel*rloc2*dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
1187                  end if
1188                end if
1189 
1190                fxgau=fxgau+forceloc*x
1191                fygau=fygau+forceloc*y
1192                fzgau=fzgau+forceloc*z
1193 
1194 !            Long range part: error function
1195                Vel=pot(ind)
1196                fxerf=fxerf+xp*Vel*x
1197                fyerf=fyerf+xp*Vel*y
1198                fzerf=fzerf+xp*Vel*z
1199 
1200              else if (nloc>0) then
1201                r2=x**2+y**2+z**2
1202                arg=r2/rloc**2
1203                xp=exp(-0.5_dp*arg)
1204                tt=cprime(nloc)
1205                do iloc=nloc-1,1,-1
1206                  tt=arg*tt+cprime(iloc)
1207                end do
1208                forceleaked=forceleaked+xp*(1._dp+tt)
1209              end if
1210            end do
1211          end do
1212        end if
1213      end do
1214 
1215 !  Final result of the forces
1216      floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc**2)*fxgau
1217      floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc**2)*fygau
1218      floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc**2)*fzgau
1219 
1220    end do
1221 
1222    forceleaked=forceleaked*prefactor*hxh*hyh*hzh
1223    if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1224 
1225    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

PARENTS

SOURCE

 959 function vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
 960 
 961 
 962 !This section has been created automatically by the script Abilint (TD).
 963 !Do not modify the following lines by hand.
 964 #undef ABI_FUNC
 965 #define ABI_FUNC 'vloc_zero_mklocl'
 966 !End of the abilint section
 967 
 968  implicit none
 969 
 970 !Arguments ------------------------------------
 971 !scalars
 972  integer,intent(in) :: msz
 973  real(dp) :: vloc_zero_mklocl
 974  real(dp),intent(in)  :: charge,rloc
 975 !arrays
 976  real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
 977 
 978 !Local variables-------------------------------
 979 !scalars
 980  real(dp) :: y1,y2,y3,zz=0._dp
 981 !arrays
 982  real(dp) :: ll(3),xx(3),yy(3)
 983 
 984 ! *************************************************************************
 985 
 986 !Select 3 points x1,x2,x3 near 0
 987    if (rad(1)>1.d-10) then
 988      xx(1:3)=rad(1:3)
 989    else
 990      xx(1:3)=rad(2:4)
 991    end if
 992 
 993 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
 994    call paw_splint(msz,rad,vloc,d2vloc,3,xx,yy)
 995    call calcVloc_mklocl(y1,xx(1),rloc,charge)
 996    call calcVloc_mklocl(y2,xx(2),rloc,charge)
 997    call calcVloc_mklocl(y3,xx(3),rloc,charge)
 998    yy(1)= yy(1)-y1
 999    yy(2)= yy(2)-y2
1000    yy(3)= yy(3)-y3
1001 
1002 !Find a polynomial of the form (z=0):
1003 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1004 
1005 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1006    ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1007 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1008    ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1009 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1010    ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1011 
1012    vloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1013 
1014  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

PARENTS

      mklocl_wavelets

CHILDREN

      calcdvloc_wvl,derf_ab,paw_splint_der

SOURCE

1900 subroutine calcdVloc_wvl(yy,xx,rloc,Z)
1901 
1902 
1903 !This section has been created automatically by the script Abilint (TD).
1904 !Do not modify the following lines by hand.
1905 #undef ABI_FUNC
1906 #define ABI_FUNC 'calcdVloc_wvl'
1907 !End of the abilint section
1908 
1909  implicit none
1910 !Arguments ------------------------------------
1911 !scalars
1912  real(dp),intent(in)  :: xx,rloc,Z
1913  real(dp),intent(out) :: yy
1914 
1915 !Local variables-------------------------------
1916  !scalars
1917  real(dp):: arg,tt
1918 
1919 ! *************************************************************************
1920 
1921    arg=xx/(sqrt(2._dp)*rloc)
1922    call derf_ab(tt,arg)
1923    yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) )
1924 
1925  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

PARENTS

SOURCE

1293 function dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc)
1294 
1295 
1296 !This section has been created automatically by the script Abilint (TD).
1297 !Do not modify the following lines by hand.
1298 #undef ABI_FUNC
1299 #define ABI_FUNC 'dvloc_zero_mklocl'
1300 !End of the abilint section
1301 
1302  implicit none
1303 
1304 !Arguments ------------------------------------
1305 !scalars
1306  integer,intent(in) :: msz
1307  real(dp) :: dvloc_zero_mklocl
1308  real(dp),intent(in)  :: charge,rloc
1309 !arrays
1310  real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
1311 
1312 !Local variables-------------------------------
1313 !scalars
1314  real(dp) :: y1,y2,y3,zz=0._dp
1315 !arrays
1316  real(dp) :: ll(3),xx(3),yy(3)
1317 
1318 ! *************************************************************************
1319 
1320 !Select 3 points x1,x2,x3 near 0
1321      if (rad(1)>1.d-10) then
1322        xx(1:3)=rad(1:3)
1323      else
1324        xx(1:3)=rad(2:4)
1325      end if
1326 
1327 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
1328      call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1329      call calcdVloc_mklocl(y1,xx(1),rloc,charge)
1330      call calcdVloc_mklocl(y2,xx(2),rloc,charge)
1331      call calcdVloc_mklocl(y3,xx(3),rloc,charge)
1332      yy(1)=(yy(1)-y1)/xx(1)
1333      yy(2)=(yy(2)-y2)/xx(2)
1334      yy(3)=(yy(3)-y3)/xx(3)
1335 
1336 !Find a polynomial of the form (z=0):
1337 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1338 
1339 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1340      ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1341 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1342      ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1343 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1344      ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1345 
1346      dvloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1347 
1348    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

PARENTS

SOURCE

1945 function dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc)
1946 
1947 
1948 !This section has been created automatically by the script Abilint (TD).
1949 !Do not modify the following lines by hand.
1950 #undef ABI_FUNC
1951 #define ABI_FUNC 'dvloc_zero_wvl'
1952 !End of the abilint section
1953 
1954  implicit none
1955 
1956 !Arguments ------------------------------------
1957 !scalars
1958  integer,intent(in) :: msz
1959  real(dp) :: dvloc_zero_wvl
1960  real(dp),intent(in)  :: charge,rloc
1961 !arrays
1962  real(dp) :: rad(msz),vloc(msz),d2vloc(msz)
1963 
1964 !Local variables-------------------------------
1965 !scalars
1966  real(dp) :: y1,y2,y3,zz=0._dp
1967 !arrays
1968  real(dp) :: ll(3),xx(3),yy(3)
1969 
1970 ! *************************************************************************
1971 
1972 !Select 3 points x1,x2,x3 near 0
1973    if (rad(1)>1.d-10) then
1974      xx(1:3)=rad(1:3)
1975    else
1976      xx(1:3)=rad(2:4)
1977    end if
1978 
1979 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x
1980    call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy)
1981    call calcdVloc_wvl(y1,xx(1),rloc,charge)
1982    call calcdVloc_wvl(y2,xx(2),rloc,charge)
1983    call calcdVloc_wvl(y3,xx(3),rloc,charge)
1984    yy(1)=(yy(1)-y1)/xx(1)
1985    yy(2)=(yy(2)-y2)/xx(2)
1986    yy(3)=(yy(3)-y3)/xx(3)
1987 
1988 !Find a polynomial of the form (z=0):
1989 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z)
1990 
1991 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3))
1992    ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3)))
1993 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3))
1994    ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3)))
1995 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2))
1996    ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2)))
1997 
1998    dvloc_zero_wvl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3)
1999 
2000  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

PARENTS

      mklocl_wavelets

CHILDREN

      calcdvloc_wvl,derf_ab,paw_splint_der

SOURCE

1676 subroutine local_forces_wvl(iproc,natom,rxyz,hxh,hyh,hzh,n1,n2,n3,n3pi,i3s,n1i,n2i,&
1677 &                           rho,pot,floc,wvl)
1678 
1679  use defs_wvltypes
1680 #if defined HAVE_BIGDFT
1681  use BigDFT_API, only : PSPCODE_PAW,ind_positions
1682 #endif
1683 
1684 !This section has been created automatically by the script Abilint (TD).
1685 !Do not modify the following lines by hand.
1686 #undef ABI_FUNC
1687 #define ABI_FUNC 'local_forces_wvl'
1688 !End of the abilint section
1689 
1690  implicit none
1691 
1692 !Arguments -------------------------------
1693 !scalars
1694  integer,intent(in) :: i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom
1695  real(dp),intent(in) :: hxh,hyh,hzh
1696  type(wvl_internal_type),intent(in) :: wvl
1697 !arrays
1698  real(dp),intent(in) :: rxyz(3,natom)
1699  real(dp),dimension(*),intent(in) :: rho,pot
1700  real(dp),intent(out) :: floc(3,natom)
1701 
1702 !Local variables -------------------------
1703 #if defined HAVE_BIGDFT
1704 !scalars
1705  integer :: i1,i2,i3,iat,iex,iey,iez,iloc,ind,isx,isy,isz,ityp
1706  integer :: j1,j2,j3,msz=0
1707  integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,nloc
1708  logical :: perx,pery,perz,gox,goy,goz,USE_PAW
1709  real(dp) :: arg,charge,cutoff,dvhgh,forceloc,fxerf,fyerf,fzerf,fxgau,fygau,fzgau
1710  real(dp) :: forceleaked,prefactor,r2,rhoel,rloc,rloc2,rx,ry,rz,rzero,tt,vel,x,xp,y,z
1711 !arrays
1712  real(dp) :: cprime(4),dvpawdr(1),rr(1)
1713  real(dp),pointer :: psppar(:,:),rad(:),vloc(:),d2vloc(:)
1714 #endif
1715 
1716 ! *********************************************************************
1717 
1718 #if defined HAVE_BIGDFT
1719 
1720  if (iproc==0) write(std_out,'(a)',advance='no') ' Calculate local forces...'
1721 
1722 !PAW or NCPP ?
1723  USE_PAW=any(wvl%atoms%npspcode==PSPCODE_PAW)
1724 
1725 !Conditions for periodicity in the three directions
1726  perx=(wvl%atoms%astruct%geocode /= 'F')
1727  pery=(wvl%atoms%astruct%geocode == 'P')
1728  perz=(wvl%atoms%astruct%geocode /= 'F')
1729  call ext_buffers(perx,nbl1,nbr1)
1730  call ext_buffers(pery,nbl2,nbr2)
1731  call ext_buffers(perz,nbl3,nbr3)
1732 
1733  forceleaked=zero
1734 
1735  do iat=1,natom
1736    ityp=wvl%atoms%astruct%iatype(iat)
1737 
1738 !  Coordinates of the center
1739    rx=rxyz(1,iat)
1740    ry=rxyz(2,iat)
1741    rz=rxyz(3,iat)
1742 
1743 !  Initialization of the forces
1744 !  ion-electron term, error function part
1745    fxerf=zero
1746    fyerf=zero
1747    fzerf=zero
1748 !  ion-electron term, gaussian part
1749    fxgau=zero
1750    fygau=zero
1751    fzgau=zero
1752 
1753 !  Building array of coefficients of the derivative of the gaussian part
1754    psppar => wvl%atoms%psppar(:,:,ityp)
1755    cprime(1)=2._dp*wvl%atoms%psppar(0,2,ityp)-wvl%atoms%psppar(0,1,ityp)
1756    cprime(2)=4._dp*wvl%atoms%psppar(0,3,ityp)-wvl%atoms%psppar(0,2,ityp)
1757    cprime(3)=6._dp*wvl%atoms%psppar(0,4,ityp)-wvl%atoms%psppar(0,3,ityp)
1758    cprime(4)=-wvl%atoms%psppar(0,4,ityp)
1759 
1760 !  Determine number of local terms (HGH pot)
1761    nloc=0
1762    do iloc=1,4
1763      if (wvl%atoms%psppar(0,iloc,ityp).ne.zero) nloc=iloc
1764    end do
1765 
1766 !  Some constants depending on the atom type
1767    rloc=wvl%atoms%psppar(0,0,ityp) ; rloc2=rloc**2
1768    charge=real(wvl%atoms%nelpsp(ityp),kind=dp)
1769    prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5)
1770 
1771 !  PAW specifics
1772    if (USE_PAW) then
1773      msz=wvl%rholoc%msz(ityp)
1774      rad    => wvl%rholoc%rad(1:msz,ityp)
1775      vloc   => wvl%rholoc%d(1:msz,3,ityp)
1776      d2vloc => wvl%rholoc%d(1:msz,4,ityp)
1777      rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2)
1778    end if
1779 
1780 !  Maximum extension of the gaussian
1781    cutoff=10._dp*rloc
1782    isx=floor((rx-cutoff)/hxh)
1783    isy=floor((ry-cutoff)/hyh)
1784    isz=floor((rz-cutoff)/hzh)
1785    iex=ceiling((rx+cutoff)/hxh)
1786    iey=ceiling((ry+cutoff)/hyh)
1787    iez=ceiling((rz+cutoff)/hzh)
1788 
1789 !  Calculate the forces near the atom due to the gaussian
1790 !  and error function parts of the potential
1791    if (n3pi>0) then
1792      do i3=isz,iez
1793        z=real(i3,kind=dp)*hzh-rz
1794        call ind_positions(perz,i3,n3,j3,goz)
1795        j3=j3+nbl3+1
1796        do i2=isy,iey
1797          y=real(i2,kind=dp)*hyh-ry
1798          call ind_positions(pery,i2,n2,j2,goy)
1799          do i1=isx,iex
1800            x=real(i1,kind=dp)*hxh-rx
1801            call ind_positions(perx,i1,n1,j1,gox)
1802 
1803            r2=x**2+y**2+z**2
1804            xp=exp(-0.5_dp*r2/rloc2)
1805 
1806            if ((j3>=i3s.and.j3<=i3s+n3pi-1) .and. (gox.and.goy)) then
1807              ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s+1-1)*n1i*n2i
1808 
1809 !            Short range part
1810              rhoel=rho(ind)
1811 !            HGH: V_S^prime=gaussian
1812              if (.not.USE_PAW) then
1813                if (nloc/=0) then
1814                  arg=r2/rloc2
1815                  tt=cprime(nloc)
1816                  do iloc=nloc-1,1,-1
1817                    tt=arg*tt+cprime(iloc)
1818                  end do
1819                  forceloc=xp*tt*rhoel
1820                else
1821                  forceloc=zero
1822                end if
1823 !            PAW: V_PAW^prime-V_L^prime
1824              else
1825                rr(1)=sqrt(r2)
1826 
1827                if (rr(1)>=rzero) then
1828                  call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr)
1829                  call calcdVloc_wvl(dvhgh,rr(1),rloc,charge)
1830                  forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1)
1831                else
1832                  forceloc=rhoel*rloc2*dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc)
1833                end if
1834              end if
1835 
1836              fxgau=fxgau+forceloc*x
1837              fygau=fygau+forceloc*y
1838              fzgau=fzgau+forceloc*z
1839 
1840 !            Long range part: error function
1841              vel=pot(ind)
1842              fxerf=fxerf+xp*vel*x
1843              fyerf=fyerf+xp*vel*y
1844              fzerf=fzerf+xp*vel*z
1845 
1846            else if ((.not.goz).and.(nloc>0)) then
1847              arg=r2/rloc2
1848              tt=cprime(nloc)
1849              do iloc=nloc-1,1,-1
1850                tt=arg*tt+cprime(iloc)
1851              end do
1852              forceleaked=forceleaked+prefactor*xp*tt*rho(1)
1853            end if
1854 
1855          end do ! i1
1856        end do ! i2
1857      end do ! i3
1858    end if ! n3pi>0
1859 
1860 !  Final result of the forces
1861    floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc2)*fxgau
1862    floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc2)*fygau
1863    floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc2)*fzgau
1864 
1865  end do ! iat
1866 
1867  forceleaked=forceleaked*hxh*hyh*hzh
1868  if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked
1869 
1870 #else
1871  BIGDFT_NOTENABLED_ERROR()
1872  if (.false.) write(std_out,*) i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom,hxh,hyh,hzh,&
1873 & rxyz(1,1),floc(1,1),rho(1),pot(1),wvl%h(1)
1874 #endif
1875 
1876  CONTAINS