TABLE OF CONTENTS


ABINIT/cc_derivatives [ Functions ]

[ Top ] [ Functions ]

NAME

 cc_derivatives

FUNCTION

 subroutine to spline the core charge and get derivatives
   extracted from previous version of psp6cc_drh
 input on log grid, and splined to regular grid between 0 and rchrg

INPUTS

  mmax=maximum number of points in real space grid in the psp file
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  rchrg=cut-off radius for the core density
  rad=radial grid points
  ff=core charge at points in rad
  ff1=first derivative of ff on log grid
  ff2=second derivative of ff on log grid

OUTPUT

  xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives

PARENTS

      psp6cc_drh,upf2abinit

CHILDREN

      spline,splint

NOTES

 Test version by DRH - requires very smooth model core charge

SOURCE

1675 subroutine cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d)
1676 
1677 
1678 !This section has been created automatically by the script Abilint (TD).
1679 !Do not modify the following lines by hand.
1680 #undef ABI_FUNC
1681 #define ABI_FUNC 'cc_derivatives'
1682 !End of the abilint section
1683 
1684  implicit none
1685 
1686 !Arguments ------------------------------------
1687 ! scalars
1688  integer,intent(in) :: mmax,n1xccc
1689  real(dp),intent(in) :: rchrg
1690 !arrays
1691  real(dp),intent(in) :: rad(mmax),ff(mmax),ff1(mmax),ff2(mmax)
1692  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
1693 
1694 !Local variables-------------------------------
1695 ! scalars
1696  integer :: i1xccc
1697  real(dp) :: der1,dern
1698 !arrays
1699  real(dp),allocatable :: ff3(:),ff4(:),gg(:),gg1(:),gg2(:)
1700  real(dp),allocatable :: gg3(:),gg4(:),work(:),xx(:)
1701 
1702 ! *************************************************************************
1703  ABI_ALLOCATE(ff3,(mmax))
1704  ABI_ALLOCATE(ff4,(mmax))
1705  ABI_ALLOCATE(gg,(n1xccc))
1706  ABI_ALLOCATE(gg1,(n1xccc))
1707  ABI_ALLOCATE(gg2,(n1xccc))
1708  ABI_ALLOCATE(gg3,(n1xccc))
1709  ABI_ALLOCATE(gg4,(n1xccc))
1710  ABI_ALLOCATE(work,(mmax))
1711  ABI_ALLOCATE(xx,(n1xccc))
1712 
1713  !write(std_out,*) 'cc_derivatives : enter'
1714 
1715 !calculate third derivative ff3 on logarithmic grid
1716  der1=ff2(1)
1717  dern=ff2(mmax)
1718  call spline(rad,ff1,mmax,der1,dern,ff3)
1719 
1720 !calculate fourth derivative ff4 on logarithmic grid
1721  der1=0.d0
1722  dern=0.d0
1723  call spline(rad,ff2,mmax,der1,dern,ff4)
1724 
1725 !generate uniform mesh xx in the box cut by rchrg:
1726 
1727  do i1xccc=1,n1xccc
1728    xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1)
1729  end do
1730 !
1731 !now interpolate core charge and derivatives on the uniform grid
1732 !
1733 !core charge, input=ff,  output=gg
1734  call splint(mmax,rad,ff,ff2,n1xccc,xx,gg)
1735 
1736 !first derivative input=ff1, output=gg1
1737  call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1)
1738 
1739 !normalize gg1
1740 !gg1(:)=gg1(:)*rchrg
1741 
1742 !second derivative input=ff2, output=gg2
1743  call splint(mmax,rad,ff2,ff4,n1xccc,xx,gg2)
1744 
1745 !normalize gg2
1746 !gg2(:)=gg2(:)*rchrg**2
1747 
1748 !reallocate work otherwise the calls to spline crash (n1xccc /= mmax)
1749  ABI_DEALLOCATE(work)
1750  ABI_ALLOCATE(work,(n1xccc))
1751 
1752 !recalculate 3rd derivative consistent with spline fit to first derivative
1753 !on linear grid
1754  der1=gg2(1)
1755  dern=gg2(n1xccc)
1756  call spline(xx,gg1,n1xccc,der1,dern,gg3)
1757 
1758 !calculate 4th derivative consistent with spline fit to second derivative
1759 !on linear grid
1760  der1=0.0d0
1761  dern=0.0d0
1762  call spline(xx,gg2,n1xccc,der1,dern,gg4)
1763 
1764 !now calculate second to fourth derivative by forward differences
1765 !to avoid numerical noise uses a smoothing function
1766 !
1767 !call smooth(gg1,n1xccc,10)
1768 
1769 !gg2(n1xccc)=0.0
1770 !do i1xccc=1,n1xccc-1
1771 !gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1)
1772 !end do
1773 
1774 !call smooth(gg2,n1xccc,10)
1775 
1776 !gg3(n1xccc)=0.0
1777 !do i1xccc=1,n1xccc-1
1778 !gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1)
1779 !end do
1780 
1781 !call smooth(gg3,n1xccc,10)
1782 
1783 !gg4(n1xccc)=0.0
1784 !do i1xccc=1,n1xccc-1
1785 !gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1)
1786 !end do
1787 
1788 !call smooth(gg4,n1xccc,10)
1789 
1790 !write on xcc1d
1791 !normalize to unit range usage later in program
1792  xccc1d(:,1)=gg(:)
1793  xccc1d(:,2)=gg1(:)*rchrg
1794  xccc1d(:,3)=gg2(:)*rchrg**2
1795  xccc1d(:,4)=gg3(:)*rchrg**3
1796  xccc1d(:,5)=gg4(:)*rchrg**4
1797 !***drh test
1798 !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc
1799 !***end drh test
1800 
1801 
1802 !DEBUG
1803 !note: the normalization condition is the following:
1804 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg
1805 !
1806 !norm=0.d0
1807 !do i1xccc=1,n1xccc
1808 !norm = norm + 4.d0*pi*rchrg/dble(n1xccc-1)*&
1809 !&             xx(i1xccc)**2*xccc1d(i1xccc,1)
1810 !end do
1811 !write(std_out,*) ' norm=',norm
1812 !
1813 !write(std_out,*)' psp6cc_drh : output of core charge density and derivatives '
1814 !write(std_out,*)'   xx          gg           gg1  '
1815 !do i1xccc=1,n1xccc
1816 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
1817 !end do
1818 !write(std_out,*)'   xx          gg2          gg3  '
1819 !do i1xccc=1,n1xccc
1820 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
1821 !end do
1822 !write(std_out,*)'   xx          gg4          gg5  '
1823 !do i1xccc=1,n1xccc
1824 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
1825 !end do
1826 !write(std_out,*)' psp1cc : debug done, stop '
1827 !stop
1828 !ENDDEBUG
1829 
1830  ABI_DEALLOCATE(ff3)
1831  ABI_DEALLOCATE(ff4)
1832  ABI_DEALLOCATE(gg)
1833  ABI_DEALLOCATE(gg1)
1834  ABI_DEALLOCATE(gg2)
1835  ABI_DEALLOCATE(gg3)
1836  ABI_DEALLOCATE(gg4)
1837  ABI_DEALLOCATE(work)
1838  ABI_DEALLOCATE(xx)
1839 
1840 end subroutine cc_derivatives

ABINIT/gg1cc [ Functions ]

[ Top ] [ Functions ]

NAME

 gg1cc

FUNCTION

 gg1cc_xx=$(\frac{\sin(2\pi xx)}{(2\pi xx)(1-4xx^2)(1-xx^2)})^2$

INPUTS

  xx= abscisse to which gg1cc_xx is calculated

OUTPUT

  gg1cc_xx= gg1cc_x(xx)

PARENTS

      psp1cc

CHILDREN

SOURCE

238 subroutine gg1cc(gg1cc_xx,xx)
239 
240 
241 !This section has been created automatically by the script Abilint (TD).
242 !Do not modify the following lines by hand.
243 #undef ABI_FUNC
244 #define ABI_FUNC 'gg1cc'
245 !End of the abilint section
246 
247  implicit none
248 
249 !Arguments ------------------------------------
250 !scalars
251  real(dp),intent(in) :: xx
252  real(dp),intent(out) :: gg1cc_xx
253 
254 !Local variables -------------------------------------------
255 !The c s are coefficients for Taylor expansion of the analytic form near xx=0, 1/2, and 1.
256 !scalars
257  real(dp) :: c21=4.d0/9.d0,c22=-40.d0/27.d0,c23=20.d0/3.d0-16.d0*pi**2/27.d0
258  real(dp) :: c24=-4160.d0/243.d0+160.d0*pi**2/81.d0,c31=1.d0/36.d0
259  real(dp) :: c32=-25.d0/108.d0,c33=485.d0/432.d0-pi**2/27.d0
260  real(dp) :: c34=-4055.d0/972.d0+25.d0*pi**2/81.d0
261 
262 ! *************************************************************************
263 
264 !Cut off beyond 3/gcut=xcccrc
265  if (xx>3.0d0) then
266    gg1cc_xx=0.0d0
267 !  Take care of difficult limits near x=0, 1/2, and 1
268  else if (abs(xx)<=1.d-09) then
269    gg1cc_xx=1.d0
270  else if (abs(xx-0.5d0)<=1.d-04) then
271 !  (this limit and next are more troublesome for numerical cancellation)
272    gg1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*c24))
273  else if (abs(xx-1.d0)<=1.d-04) then
274    gg1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*c34))
275  else
276 !  The following is the square of the Fourier transform of a
277 !  function built out of two spherical bessel functions in G
278 !  space and cut off absolutely beyond gcut
279    gg1cc_xx=(sin(2.0d0*pi*xx)/( (2.0d0*pi*xx) * &
280 &   (1.d0-4.0d0*xx**2)*(1.d0-xx**2) )  )**2
281  end if
282 
283 end subroutine gg1cc

ABINIT/gp1cc [ Functions ]

[ Top ] [ Functions ]

NAME

 gp1cc

FUNCTION

 Derivative of gg(xx) wrt xx.

INPUTS

  xx=abscisse to which gp1cc_xx is calculated

OUTPUT

  gp1cc_xx=derivative of gg(xx) wrt xx.

NOTES

 $ phi(x) = \frac{\sin(2\pi x)}{(2\pi x)(1-4x^2)(1-x^2)}$
 $ gg(x)= phi(x)^2$
 $ gp(x)= 2 * phi(x) * phi''(x)$
 $ phi''(x)=\frac{\cos(2\pi x)-(1-15x^2+20x^4) phi(x)}{x(1-4x^2)(1-x^2)}$

PARENTS

      psp1cc

CHILDREN

SOURCE

312 subroutine gp1cc(gp1cc_xx,xx)
313 
314 
315 !This section has been created automatically by the script Abilint (TD).
316 !Do not modify the following lines by hand.
317 #undef ABI_FUNC
318 #define ABI_FUNC 'gp1cc'
319 !End of the abilint section
320 
321  implicit none
322 
323 !Arguments ------------------------------------
324 !scalars
325  real(dp),intent(in) :: xx
326  real(dp),intent(out) :: gp1cc_xx
327 
328 !Local variables -------------------------------------------
329 !scalars
330  real(dp),parameter :: c11=20.d0-8.d0*pi**2/3.d0
331  real(dp),parameter :: c12=268.d0-160.d0/3.d0*pi**2+128.d0/45.d0*pi**4
332  real(dp),parameter :: c21=-40.d0/27.d0,c22=40.d0/3.d0-32.d0*pi**2/27.d0
333  real(dp),parameter :: c23=-4160.d0/81.d0+160.d0*pi**2/27.d0
334  real(dp),parameter :: c24=157712.d0/729.d0-320.d0*pi**2/9.d0+512.d0*pi**4/405.d0
335  real(dp),parameter :: c25=-452200.d0/729.d0+83200.d0*pi**2/729.d0-1280.d0*pi**4/243.d0
336  real(dp),parameter :: c31=-25.d0/108.d0,c32=485.d0/216.d0-2.d0*pi**2/27.d0
337  real(dp),parameter :: c33=-4055.d0/324.d0+25.d0*pi**2/27.d0
338  real(dp),parameter :: c34=616697.d0/11664.d0-485.d0*pi**2/81.d0+32.d0*pi**4/405.d0
339  real(dp),parameter :: c35=-2933875.d0/15552.d0+20275.d0*pi**2/729.d0-200.d0*pi**4/243.d0
340  real(dp),parameter :: two_pim1=1.0d0/two_pi
341  real(dp) :: denom,phi,phip
342 
343 ! *************************************************************************
344 
345 !Cut off beyond r=3*xcccrc is already done at the calling level
346  if (xx>1.001d0) then
347 !  The part that follows will be repeated later, but written in this way,
348 !  only one "if" condition is tested in most of the cases (1.001 < x < 3.0)
349    denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
350    phi=denom*sin(two_pi*xx)*two_pim1
351    phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
352    gp1cc_xx=2.d0*phi*phip
353 !  Handle limits where denominator vanishes
354  else if (abs(xx)<1.d-03) then
355    gp1cc_xx=xx*(c11+xx**2*c12)
356  else if (abs(xx-0.5d0)<=1.d-03) then
357    gp1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*(c24+(xx-0.5d0)*c25)))
358  else if (abs(xx-1.d0)<=1.d-03) then
359    gp1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*(c34+(xx-1.0d0)*c35)))
360  else
361 !  Here is the repeated part ...
362    denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
363    phi=denom*sin(two_pi*xx)*two_pim1
364    phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
365    gp1cc_xx=2.d0*phi*phip
366  end if
367 
368 end subroutine gp1cc

ABINIT/gpp1cc [ Functions ]

[ Top ] [ Functions ]

NAME

 gpp1cc

FUNCTION

 Second derivative of gg wrt xx.

INPUTS

  xx= abscisse to which gpp1cc_xx is calculated

OUTPUT

  gpp1cc_xx=second derivative of gg wrt xx.

PARENTS

      psp1cc

CHILDREN

SOURCE

392 subroutine gpp1cc(gpp1cc_xx,xx)
393 
394 
395 !This section has been created automatically by the script Abilint (TD).
396 !Do not modify the following lines by hand.
397 #undef ABI_FUNC
398 #define ABI_FUNC 'gpp1cc'
399 !End of the abilint section
400 
401  implicit none
402 
403 !Arguments ------------------------------------
404 !scalars
405  real(dp),intent(in) :: xx
406  real(dp),intent(out) :: gpp1cc_xx
407 
408 !Local variables -------------------------------------------
409 !scalars
410  real(dp),parameter :: c1=20.d0-8.d0*pi**2/3.d00
411  real(dp),parameter :: c2=40.d0/3.d0-32.d0*pi**2/27.d0
412  real(dp),parameter :: c3=-8320.d0/81.d0+320.d0*pi**2/27.d0
413  real(dp),parameter :: c4=157712.d0/243.d0-320.d0*pi**2/3.d0+512.d0*pi**4/135.d0
414  real(dp),parameter :: c5=-18088.d2/729.d0+3328.d2*pi**2/729.d0-5120.d0*pi**4/243.d0
415  real(dp),parameter :: c6=485.d0/216.d0-2.d0*pi**2/27.d0
416  real(dp),parameter :: c7=-4055.d0/162.d0+50.d0*pi**2/27.d0
417  real(dp),parameter :: c8=616697.d0/3888.d0-485.d0*pi**2/27.d0+32.d0*pi**4/135.d0
418  real(dp),parameter :: c9=-2933875.d0/3888.d0+81100.d0*pi**2/729.d0-800.d0*pi**4/243.d0
419  real(dp) :: t1,t10,t100,t11,t12,t120,t121,t122,t127,t138,t14,t140,t15,t152
420  real(dp) :: t157,t16,t160,t17,t174,t175,t18,t19,t2,t20,t21,t23,t24,t3,t31,t33
421  real(dp) :: t34,t4,t41,t42,t44,t45,t46,t5,t54,t55,t56,t57,t6,t62,t64,t65,t7
422  real(dp) :: t72,t78,t79,t8,t85,t9,t93
423 
424 ! *************************************************************************
425 
426  if (xx>3.0d0) then
427 !  Cut off beyond 3/gcut=3*xcccrc
428    gpp1cc_xx=0.0d0
429 !  Take care of difficult limits near xx=0, 1/2, and 1
430  else if (abs(xx)<=1.d-09) then
431    gpp1cc_xx=c1
432  else if (abs(xx-0.5d0)<=1.d-04) then
433 !  (this limit and next are more troublesome for numerical cancellation)
434    gpp1cc_xx=c2+(xx-0.5d0)*(c3+(xx-0.5d0)*(c4+(xx-0.5d0)*c5))
435  else if (abs(xx-1.d0)<=1.d-04) then
436    gpp1cc_xx=c6+(xx-1.0d0)*(c7+(xx-1.0d0)*(c8+(xx-1.0d0)*c9))
437  else
438 
439 !  Should fix up this Maple fortran later
440    t1 = xx**2
441    t2 = 1/t1
442    t3 = 1/Pi
443    t4 = 2*xx
444    t5 = t4-1
445    t6 = t5**2
446    t7 = 1/t6
447    t8 = t4+1
448    t9 = t8**2
449    t10 = 1/t9
450    t11 = xx-1
451    t12 = t11**2
452    t14 = 1/t12/t11
453    t15 = xx+1
454    t16 = t15**2
455    t17 = 1/t16
456    t18 = Pi*xx
457    t19 = sin(t18)
458    t20 = cos(t18)
459    t21 = t20**2
460    t23 = t19*t21*t20
461    t24 = t17*t23
462    t31 = t19**2
463    t33 = t31*t19*t20
464    t34 = t17*t33
465    t41 = Pi**2
466    t42 = 1/t41
467    t44 = 1/t16/t15
468    t45 = t31*t21
469    t46 = t44*t45
470    t54 = 1/t1/xx
471    t55 = 1/t12
472    t56 = t55*t46
473    t57 = t10*t56
474    t62 = t9**2
475    t64 = t17*t45
476    t65 = t55*t64
477    t72 = 1/t9/t8
478    t78 = t14*t64
479    t79 = t10*t78
480    t85 = t12**2
481    t93 = t21**2
482    t100 = t31**2
483    t120 = 1/t6/t5
484    t121 = t55*t34
485    t122 = t10*t121
486    t127 = t16**2
487    t138 = t6**2
488    t140 = t10*t65
489    t152 = t72*t65
490    t157 = t7*t140
491    t160 = t1**2
492    t174 = t55*t24
493    t175 = t10*t174
494    gpp1cc_xx = 8*t2*t3*t7*t10*t14*t34+8*t2*t42*t7*t10*t14*t46&
495 &   -8*t2*t3*t7*t10*t14*t24+8*t2*t3*t7*t10*t55*t44*t33+&
496 &   6*t2*t42*t7*t10*t55/t127*t45+24*t2*t42/t138*t140+&
497 &   16*t54*t42*t120*t140+16*t2*t3*t120*t122+16*t2&
498 &   *t42*t7*t72*t78-8*t2*t3*t7*t10*t55*t44*t23-8*t54*t3*t7*t175&
499 &   +2*t2*t7*t10*t55*t17*t100+2*t2*t7*t10*t55*t17*t93+&
500 &   8*t54*t42*t7*t79+16*t2*t42*t7*t72*t56+6*t2*t42*t7*t10/t85&
501 &   *t64+24*t2*t42*t7/t62*t65+8*t54*t42*t7*t57-&
502 &   16*t2*t3*t7*t72*t174+8*t54*t3*t7*t122-16*t2*t3*t120*t175&
503 &   +16*t2*t42*t120*t79+16*t2*t42*t120*t57+16*t54*t42*t7*t152+&
504 &   32*t2*t42*t120*t152+16*t2*t3*t7*t72*t121-12*t2*t157+&
505 &   6/t160*t42*t157
506  end if
507 
508 end subroutine gpp1cc

ABINIT/m_psptk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psptk

FUNCTION

  This module collects low-level procedures used by the other psp modules

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (XG, DCA, MM, DRH, FrD, GZ, AF)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_psptk
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_splines
33 
34  use m_numeric_tools,  only : ctrap
35  use m_special_funcs,  only : sbf8
36 
37  implicit none
38 
39  private

ABINIT/psp1cc [ Functions ]

[ Top ] [ Functions ]

NAME

 psp1cc

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for the format 1 and 5 of pseudopotentials.
 WARNING : the fifth derivate is actually set to zero

INPUTS

  fchrg=magnitude of the core charge correction
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used

OUTPUT

  xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives

NOTES

 This is a revised expression for core density (5 Nov 1992) :
 density(r)=fchrg*gg(xx)
 with
 $ gg(xx)=(\frac{\sin(2\pi xx)}{(2\pi xx)(1-4 xx^2)(1-xx^2)})^2 $
 and
 $ xx=\frac{r}{rchrg}=\frac{r}{xcccrc/3.0d0}=3*\frac{r}{xcccrc}=3*yy $

 Code for gg(xx), gp(xx), and gpp(xx) has been tested by numerical
 derivatives--looks ok. gpp(x) should still be rewritten.
 The argument of xccc1d is assumed to be normalized, and to vary
 from yy=0 to 1 (from r=0 to r=xcccrc, or from xx=0 to 3)
 Thus :
 xccc1d(yy)=fchrg*[\frac{\sin(2*\pi*(3yy))}
 {(6*\pi*(3yy))(1-4*(3yy)^2)(1-(3yy)^2)}]^2

WARNINGS

 Warning : the fifth derivative is not yet delivered.

PARENTS

      psp1in,psp5in

CHILDREN

      gg1cc,gp1cc,gpp1cc,spline

SOURCE

 99 subroutine psp1cc(fchrg,n1xccc,xccc1d)
100 
101 
102 !This section has been created automatically by the script Abilint (TD).
103 !Do not modify the following lines by hand.
104 #undef ABI_FUNC
105 #define ABI_FUNC 'psp1cc'
106 !End of the abilint section
107 
108  implicit none
109 
110 !Arguments ------------------------------------
111 !scalars
112  integer,intent(in) :: n1xccc
113  real(dp),intent(in) :: fchrg
114 !arrays
115  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
116 
117 !Local variables-------------------------------
118 !scalars
119  integer :: i1xccc,ider
120  real(dp) :: der1,dern,factor,gg1cc_xx,gp1cc_xx,gpp1cc_xx,xx
121  character(len=500) :: message
122 !arrays
123  real(dp),allocatable :: ff(:),ff2(:),work(:),yy(:)
124 
125 ! *************************************************************************
126 
127  ABI_ALLOCATE(ff,(n1xccc))
128  ABI_ALLOCATE(ff2,(n1xccc))
129  ABI_ALLOCATE(work,(n1xccc))
130  ABI_ALLOCATE(yy,(n1xccc))
131 
132  if(n1xccc > 1)then
133    factor=one/dble(n1xccc-1)
134    do i1xccc=1,n1xccc
135      yy(i1xccc)=(i1xccc-1)*factor
136    end do
137  else
138    write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc
139    MSG_BUG(message)
140  end if
141 
142 !Initialization, to avoid some problem with some compilers
143  xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
144 
145 !Take care of each derivative separately
146  do ider=0,2
147 
148    if(ider==0)then
149 !    Generate spline fitting for the function gg
150      do i1xccc=1,n1xccc
151        xx=three*yy(i1xccc)
152        call gg1cc(gg1cc_xx,xx)
153        ff(i1xccc)=fchrg*gg1cc_xx
154      end do
155 !    Complete with derivatives at end points
156      der1=zero
157      call gp1cc(gp1cc_xx,three)
158      dern=three*fchrg*gp1cc_xx
159    else if(ider==1)then
160 !    Generate spline fitting for the function gp
161      do i1xccc=1,n1xccc
162        xx=three*yy(i1xccc)
163        call gp1cc(gp1cc_xx,xx)
164        ff(i1xccc)=three*fchrg*gp1cc_xx
165      end do
166 !    Complete with derivatives at end points, already estimated
167      der1=xccc1d(1,ider+2)
168      dern=xccc1d(n1xccc,ider+2)
169    else if(ider==2)then
170 !    Generate spline fitting for the function gpp
171 !    (note : the function gpp has already been estimated, for the spline
172 !    fitting of the function gg, but it is replaced here by the more
173 !    accurate analytic derivative)
174      do i1xccc=1,n1xccc
175        xx=three*yy(i1xccc)
176        call gpp1cc(gpp1cc_xx,xx)
177        ff(i1xccc)=9.0_dp*fchrg*gpp1cc_xx
178      end do
179 !    Complete with derivatives of end points
180      der1=xccc1d(1,ider+2)
181      dern=xccc1d(n1xccc,ider+2)
182    end if
183 
184 !  Produce second derivative numerically, for use with splines
185    call spline(yy,ff,n1xccc,der1,dern,ff2)
186    xccc1d(:,ider+1)=ff(:)
187    xccc1d(:,ider+3)=ff2(:)
188  end do
189 
190  xccc1d(:,6)=zero
191 
192 !DEBUG
193 !write(std_out,*)' psp1cc : output of core charge density and derivatives '
194 !write(std_out,*)'   yy          gg           gp  '
195 !do i1xccc=1,n1xccc
196 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
197 !end do
198 !write(std_out,*)'   yy          gpp          gg2  '
199 !do i1xccc=1,n1xccc
200 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
201 !end do
202 !write(std_out,*)'   yy          gp2          gpp2  '
203 !do i1xccc=1,n1xccc
204 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
205 !end do
206 !write(std_out,*)' psp1cc : debug done, stop '
207 !stop
208 !ENDDEBUG
209 
210  ABI_DEALLOCATE(ff)
211  ABI_DEALLOCATE(ff2)
212  ABI_DEALLOCATE(work)
213  ABI_DEALLOCATE(yy)
214 
215 end subroutine psp1cc

ABINIT/psp5lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp5lo

FUNCTION

 Compute sine transform to transform from V(r) to q^2 V(q).
 Computes integrals on logarithmic grid using related uniform
 grid in exponent and corrected trapezoidal integration.

INPUTS

  al=spacing in exponent for radial atomic grid.
  mmax=number of radial r grid points (logarithmic atomic grid).
  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=q grid values (bohr**-1).
  rad(mmax)=r grid values (bohr).
  vloc(mmax)=V(r) on radial grid.
  zion=nominal valence charge of atom.

OUTPUT

  epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
  q2vq(mqgrid)
   =q^2 V(q)
   = -\frac{Zv}{\pi}
     + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
  yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax
   (needed for spline fitter).

PARENTS

      psp5in,psp6in

CHILDREN

      ctrap

SOURCE

548 subroutine psp5lo(al,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
549 &                  vloc,yp1,ypn,zion)
550 
551 
552 !This section has been created automatically by the script Abilint (TD).
553 !Do not modify the following lines by hand.
554 #undef ABI_FUNC
555 #define ABI_FUNC 'psp5lo'
556 !End of the abilint section
557 
558  implicit none
559 
560 !Arguments----------------------------------------------------------
561 !scalars
562  integer,intent(in) :: mmax,mqgrid
563  real(dp),intent(in) :: al,zion
564  real(dp),intent(out) :: epsatm,yp1,ypn
565 !arrays
566  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
567  real(dp),intent(out) :: q2vq(mqgrid)
568 
569 !Local variables-------------------------------
570 !scalars
571  integer :: iq,ir
572  real(dp),parameter :: scale=10.0d0
573  real(dp) :: arg,result,rmtoin,test,ztor1
574 !arrays
575  real(dp),allocatable :: work(:)
576 
577 ! *************************************************************************
578 
579  ABI_ALLOCATE(work,(mmax))
580 
581 !Do q=0 separately (compute epsatm)
582 !Do integral from 0 to r1
583  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
584 
585 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
586 !with extra factor of r to convert to uniform grid in exponent
587  do ir=1,mmax
588 !  First handle tail region
589    test=vloc(ir)+zion/rad(ir)
590 !  DEBUG
591 !  write(std_out,*)ir,rad(ir),test
592 !  ENDDEBUG
593 !  Ignore small contributions, or impose a cut-off in the case
594 !  the pseudopotential data are in single precision.
595 !  (it is indeed expected that vloc is very close to zero beyond 20,
596 !  so a value larger than 2.0d-8 is considered anomalous)
597    if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
598      work(ir)=zero
599    else
600      work(ir)=(rad(ir)*rad(ir))*(rad(ir)*vloc(ir)+zion)
601    end if
602  end do
603 !DEBUG
604 !write(std_out,*)' psp5lo : stop '
605 !stop
606 !ENDDEBUG
607 
608 !Do integral from r(1) to r(max)
609  call ctrap(mmax,work,al,result)
610 !Do integral from r(mmax) to infinity
611 !compute decay length lambda at r(mmax)
612 !$\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ &
613 !$(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$
614 !rmtoin=$(rad(mmax)*vloc(mmax)+zion)*(rad(mmax)+1.d0/\lambda)/\lambda$
615 !Due to inability to fit exponential decay to r*V(r)+Zv
616 !in tail, NO TAIL CORRECTION IS APPLIED
617 !(numerical trouble might be removed if atomic code is
618 !cleaned up in tail region)
619  rmtoin=0.0d0
620 
621  epsatm=4.d0*pi*(result+ztor1+rmtoin)
622 
623  q2vq(1)=-zion/pi
624 
625 !Loop over q values
626  do iq=2,mqgrid
627    arg=2.d0*pi*qgrid(iq)
628 !  ztor1=$ -Zv/\pi+2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
629    ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion)* &
630 &   cos(arg*rad(1)) )/pi
631 
632 !  set up integrand
633    do  ir=1,mmax
634      test=vloc(ir)+zion/rad(ir)
635 !    Ignore contributions within decade of machine precision
636      if ((scale+abs(test)).eq.scale) then
637        work(ir)=zero
638      else
639        work(ir)=rad(ir)*sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
640      end if
641    end do
642 !  do integral from r(1) to r(mmax)
643    call ctrap(mmax,work,al,result)
644 
645 !  do integral from r(mmax) to infinity
646 !  rmtoin=(r(mmax)*vr(mmax)+zion)*(lambda*sin(arg*r(mmax))+
647 !  arg*cos(arg*r(mmax)))/(arg**2+lambda**2)
648 !  See comment above; no tail correction
649    rmtoin=0.0d0
650 
651 !  store q^2 v(q)
652    q2vq(iq)=ztor1+2.d0*qgrid(iq)*(result+rmtoin)
653 
654  end do
655 
656 !Compute derivatives of q^2 v(q) at ends of interval
657  yp1=0.0d0
658 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
659 !integral from 0 to r1
660  arg=2.0d0*pi*qgrid(mqgrid)
661  ztor1=zion*rad(1)*sin(arg*rad(1))
662  ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + &
663 & (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1))
664 !integral from r(mmax) to infinity is overkill; ignore
665 !set up integrand
666  do ir=1,mmax
667    test=vloc(ir)+zion/rad(ir)
668 !  Ignore contributions within decade of machine precision
669    if ((scale+abs(test)).eq.scale) then
670      work(ir)=0.0d0
671    else
672      work(ir)=rad(ir)*(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * &
673 &     (rad(ir)*vloc(ir)+zion)
674    end if
675  end do
676  call ctrap(mmax,work,al,result)
677  ypn=2.0d0 * (ztor1 + result)
678 
679  ABI_DEALLOCATE(work)
680 
681 end subroutine psp5lo

ABINIT/psp5nl [ Functions ]

[ Top ] [ Functions ]

NAME

 psp5nl

FUNCTION

 Make Kleinman-Bylander form factors f_l(q) for each l from
 0 to lmax; Vloc is assumed local potential.

INPUTS

  al=grid spacing in exponent for radial grid
  lmax=maximum ang momentum for which nonlocal form factor is desired.
   Usually lmax=1, sometimes = 0 (e.g. for oxygen); lmax <= 2 allowed.
  mmax=number of radial grid points for atomic grid
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=number of grid points for q grid
  qgrid(mqgrid)=values at which form factors are returned
  rad(mmax)=radial grid values
  vloc(mmax)=local pseudopotential on radial grid
  vpspll(mmax,3)=nonlocal pseudopotentials for each l on radial grid
  wfll(mmax,3)=reference state wavefunctions on radial grid
                mmax and mqgrid

OUTPUT

  ekb(mpsang)=Kleinman-Bylander energy,
             {{\\ \begin{equation}
               \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
               {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))   dr]}
               \end{equation} }}
                for each l
  ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum

NOTES

 u_l(r) is reference state wavefunction (input as wf);
 j_l(q) is a spherical Bessel function;
 dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
 f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
 where dvms = $\int_0^{rmax} [(u_l(r) dV_l(r))^2 dr]$ is the mean
 square value of the nonlocal correction for angular momentum l.
 Xavier Gonze s E_KB = $ dvms/\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]$.
 This is the eigenvalue of the Kleinman-Bylander operator and sets
 the energy scale of the nonlocal psp corrections.

PARENTS

      psp5in,psp6in

CHILDREN

      ctrap,spline

SOURCE

 735 subroutine psp5nl(al,ekb,ffspl,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll)
 736 
 737 
 738 !This section has been created automatically by the script Abilint (TD).
 739 !Do not modify the following lines by hand.
 740 #undef ABI_FUNC
 741 #define ABI_FUNC 'psp5nl'
 742 !End of the abilint section
 743 
 744  implicit none
 745 
 746 !Arguments ------------------------------------
 747 !scalars
 748  integer,intent(in) :: lmax,mmax,mpsang,mqgrid
 749  real(dp),intent(in) :: al
 750 !arrays
 751  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax),vpspll(mmax,mpsang)
 752  real(dp),intent(in) :: wfll(mmax,mpsang)
 753  real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang)
 754 
 755 !Local variables-------------------------------
 756 !scalars
 757  integer,parameter :: dpsang=5
 758  integer :: iq,ir,lp1
 759  real(dp) :: arg,bessel,dvwf,qr,result,yp1,ypn,ztor1
 760  character(len=500) :: message
 761 !arrays
 762  real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang)
 763  real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:)
 764 
 765 !*************************************************************************
 766 
 767 !l=0,1,2 and 3 spherical Bessel functions
 768 !The accuracy of the bes1, bes2, bes3 functions for small arguments
 769 !may be insufficient. In the present version
 770 !of the routines, some care is taken with the value of the argument.
 771 !If smaller than 1.d-3, a two terms
 772 !Taylor series expansion is prefered.
 773 ! bes0(arg)=sin(arg)/arg
 774 ! bes1(arg)=(sin(arg)-arg*cos(arg))/arg**2
 775 ! bes2(arg)=( (3.0d0-arg**2)*sin(arg)-&
 776 !& 3.0d0*arg*cos(arg) )      /arg**3
 777 
 778 ! bes3(arg)=(15.d0*sin(arg)-15.d0*arg*cos(arg) &
 779 !& -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4
 780 
 781 !Zero out Kleinman-Bylander energies ekb
 782  ekb(:)=0.0d0
 783 
 784  ABI_ALLOCATE(work1,(mmax))
 785  ABI_ALLOCATE(work2,(mmax))
 786  ABI_ALLOCATE(work3,(mmax))
 787  ABI_ALLOCATE(work4,(mmax))
 788 
 789 !Allow for no nonlocal correction (lmax=-1)
 790  if (lmax/=-1) then
 791 
 792 !  Check that lmax is within allowed range
 793    if (lmax<0.or.lmax>3) then
 794      write(message, '(a,i12,a,a,a,a,a,a,a)' )&
 795 &     'lmax=',lmax,' is not an allowed value.',ch10,&
 796 &     'Allowed values are -1 for no nonlocal correction or else',ch10,&
 797 &     '0, 1,2 or 3 for maximum l nonlocal correction.',ch10,&
 798 &     'Action: check the input atomic psp data file for lmax.'
 799      MSG_ERROR(message)
 800    end if
 801 
 802 !  Compute normalizing integrals eta=<dV> and mean square
 803 !  nonlocal psp correction dvms=<dV^2>
 804 !  "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction
 805    do lp1=1,lmax+1
 806 
 807 !    integral from 0 to r1
 808      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)
 809      ztor1=(wfll(1,lp1)*dvwf)*rad(1)/dble(2*(lp1-1)+3)
 810 !    integrand for r1 to r(mmax) (incl extra factor of r)
 811      do ir=1,mmax
 812        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
 813        work1(ir)=rad(ir)*(wfll(ir,lp1)*dvwf)
 814      end do
 815 !    do integral by corrected trapezoidal integration
 816      call ctrap(mmax,work1,al,result)
 817      eta(lp1)=ztor1+result
 818 
 819      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)
 820      ztor1=dvwf**2*rad(1)/dble(2*(lp1-1)+3)
 821      do ir=1,mmax
 822        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
 823        work1(ir)=rad(ir)*(dvwf**2)
 824      end do
 825      call ctrap(mmax,work1,al,result)
 826      dvms(lp1)=ztor1+result
 827 
 828 !    DEBUG
 829 !    Compute the norm of wfll
 830 !    wf=wfll(1,lp1)
 831 !    ztor1=wf**2*rad(1)/dble(2*(lp1-1)+3)
 832 !    do ir=1,mmax
 833 !    wf=wfll(ir,lp1)
 834 !    work1(ir)=rad(ir)*(wf**2)
 835 !    end do
 836 !    call ctrap(mmax,work1,al,result)
 837 !    norm=ztor1+result
 838 !    write(std_out,*)' lp1, norm',lp1,norm
 839 !    ENDDEBUG
 840 
 841 !    If dvms is not 0 for any given angular momentum l,
 842 !    compute Xavier Gonze s definition of the Kleinman-Bylander
 843 !    energy E_KB = dvms/eta.  In this case also renormalize
 844 !    the projection operator to u_KB(r)=$u_l(r)*dV(r)/\sqrt{dvms}$.
 845 !    This means dvwf gets multiplied by the normalization factor
 846 !    "renorm"=$1/\sqrt{dvms}$ as seen below.
 847      if (dvms(lp1)/=0.0d0) then
 848        ekb(lp1)=dvms(lp1)/eta(lp1)
 849        renorm(lp1)=1.0d0/sqrt(dvms(lp1))
 850 !      ckb is Kleinman-Bylander "cosine" (Xavier Gonze)
 851        ckb(lp1)=eta(lp1)/sqrt(dvms(lp1))
 852      else
 853        ekb(lp1)=0.0d0
 854      end if
 855 
 856    end do
 857 
 858 !  l=0 form factor if ekb(1) not 0 (lmax always at least 0)
 859    if (ekb(1)/=0.0d0) then
 860 
 861 !    do q=0 separately
 862      lp1=1
 863 !    0 to r1 integral
 864      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 865      ztor1=(rad(1)*dvwf)*rad(1)/3.0d0
 866 !    integrand
 867      do ir=1,mmax
 868        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 869        work1(ir)=rad(ir)*(rad(ir)*dvwf)
 870      end do
 871      call ctrap(mmax,work1,al,result)
 872      ffspl(1,1,1)=ztor1+result
 873 
 874 !    do rest of q points
 875      do iq=2,mqgrid
 876        arg=two_pi*qgrid(iq)
 877 !      0 to r1 integral
 878        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 879        ztor1=(bes0_psp5(arg*rad(1))*rad(1)*dvwf)*rad(1)/3.0d0
 880 !      integrand
 881        do ir=1,mmax
 882          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 883          work1(ir)=rad(ir)*(rad(ir)*bes0_psp5(arg*rad(ir))*dvwf)
 884        end do
 885        call ctrap(mmax,work1,al,result)
 886        ffspl(iq,1,1)=ztor1+result
 887      end do
 888 
 889 !    Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
 890 !    yp1=0 for l=0
 891      yp1=0.0d0
 892 !    ypn=$ \int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$
 893      arg=two_pi*qgrid(mqgrid)
 894      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 895      qr=arg*rad(1)
 896      if(qr<1.d-3)then
 897        bessel=(10.d0-qr*qr)*qr/30.0d0
 898      else
 899        bessel=bes1_psp5(qr)
 900      end if
 901 !    ztor1=(-bes1(arg*rad(1))*two_pi*rad(1)*r(1)*dvwf)*rad(1)/5.0d0
 902      ztor1=(-bessel*two_pi*rad(1)*rad(1)*dvwf)*rad(1)/5.0d0
 903      do ir=1,mmax
 904        qr=arg*rad(ir)
 905        if(qr<1.d-3)then
 906          bessel=(10.d0-qr*qr)*qr/30.0d0
 907        else
 908          bessel=bes1_psp5(qr)
 909        end if
 910        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 911 !      work(ir)=rad(ir)*(-bes1(arg*rad(ir))*two_pi*rad(ir)*rad(ir)*dvwf)
 912        work1(ir)=rad(ir)*(-bessel*two_pi*rad(ir)*rad(ir)*dvwf)
 913      end do
 914      call ctrap(mmax,work1,al,result)
 915      ypn=ztor1+result
 916 
 917 !    Fit spline to get second derivatives by spline fit
 918      call spline(qgrid,ffspl(1,1,1),mqgrid,yp1,ypn,ffspl(1,2,1))
 919 
 920    else
 921 !    or else put nonlocal correction at l=0 to 0
 922      ffspl(:,:,1)=0.0d0
 923    end if
 924 
 925 !  Finished if lmax=0 (highest nonlocal correction)
 926 !  Do l=1 form factor if ekb(2) not 0 and lmax>=1
 927    if (lmax>0)then
 928      if(ekb(2)/=0.0d0) then
 929 
 930        lp1=2
 931 !      do q=0 separately: f_1(q=0) vanishes !
 932        ffspl(1,1,2)=0.0d0
 933 
 934 !      do rest of q points
 935        do iq=2,mqgrid
 936          arg=two_pi*qgrid(iq)
 937          dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 938          qr=arg*rad(1)
 939          if(qr<1.d-3)then
 940            bessel=(10.d0-qr*qr)*qr/30.0d0
 941          else
 942            bessel=bes1_psp5(qr)
 943          end if
 944 !        ztor1=(bes1(arg*rad(1))*rad(1)*dvwf)*rad(1)/5.0d0
 945          ztor1=(bessel*rad(1)*dvwf)*rad(1)/5.0d0
 946 
 947          do ir=1,mmax
 948            qr=arg*rad(ir)
 949            if(qr<1.d-3)then
 950              bessel=(10.d0-qr*qr)*qr/30.0d0
 951            else
 952              bessel=bes1_psp5(qr)
 953            end if
 954            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 955            work2(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
 956          end do
 957 
 958          call ctrap(mmax,work2,al,result)
 959          ffspl(iq,1,2)=ztor1+result
 960        end do
 961 
 962 !      Compute yp1,ypn for l=1
 963 !      yp1=$\displaystyle \int [2\pi r^2 wf(r) dV(r)]/3$
 964        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 965        ztor1=((two_pi*rad(1)**2)*dvwf)*rad(1)/(3.0d0*5.0d0)
 966        do ir=1,mmax
 967          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 968          work2(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf/3.0d0)
 969        end do
 970        call ctrap(mmax,work2,al,result)
 971        yp1=ztor1+result
 972 !      ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$
 973 !      where x=2 Pi qgrid(mqgrid) r
 974        arg=two_pi*qgrid(mqgrid)
 975        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 976        qr=arg*rad(1)
 977        if(qr<1.d-3)then
 978          bessel=(10.d0-3.0d0*qr*qr)/30.0d0
 979        else
 980          bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr
 981        end if
 982 !      ztor1=( (two_pi*rad(1)**2)*dvwf* (bes0(arg*rad(1))-
 983 !      2.0d0*bes1(arg*rad(1))/(arg*rad(1))) ) * rad(1)/5.0d0
 984        ztor1=( (two_pi*rad(1)**2)*dvwf*bessel)*  rad(1)/5.0d0
 985 
 986        do ir=1,mmax
 987          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 988          qr=arg*rad(ir)
 989          if(qr<1.d-3)then
 990            bessel=(10.d0-3.0d0*qr*qr)/30.0d0
 991          else
 992            bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr
 993          end if
 994 !        work(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
 995 !        (bes0(arg*rad(ir))-2.d0*bes1(arg*rad(ir))/(arg*rad(ir))) )
 996          work2(ir)=rad(ir)*(two_pi*rad(ir)**2)*dvwf*bessel
 997        end do
 998        call ctrap(mmax,work2,al,result)
 999        ypn=ztor1+result
1000 
1001 !      Fit spline for l=1 Kleinman-Bylander form factor
1002        call spline(qgrid,ffspl(1,1,2),mqgrid,yp1,ypn,ffspl(1,2,2))
1003 
1004      else
1005 !      or else put form factor to 0 for l=1
1006        ffspl(:,:,2)=0.0d0
1007      end if
1008 !    Endif condition of lmax>0
1009    end if
1010 
1011 !  Finished if lmax=1 (highest nonlocal correction)
1012 !  Do l=2 nonlocal form factor if eta(3) not 0 and lmax>=2
1013    if (lmax>1)then
1014      if(ekb(3)/=0.0d0) then
1015 
1016        lp1=3
1017 !      do q=0 separately; f_2(q=0) vanishes
1018        ffspl(1,1,3)=0.0d0
1019 
1020 !      do rest of q points
1021        do iq=2,mqgrid
1022          arg=two_pi*qgrid(iq)
1023          dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1024          qr=arg*rad(1)
1025          if(qr<1.d-3)then
1026            bessel=qr*qr/15.0d0-qr**4/210.0d0
1027          else
1028            bessel=bes2_psp5(qr)
1029          end if
1030 !        ztor1=(bes2(arg*rad(1))*rad(1)*dvwf)*rad(1)/7.0d0
1031          ztor1=(bessel*rad(1)*dvwf)*rad(1)/7.0d0
1032          do ir=1,mmax
1033            qr=arg*rad(ir)
1034            if(qr<1.d-3)then
1035              bessel=qr*qr/15.0d0-qr**4/210.0d0
1036            else
1037              bessel=bes2_psp5(qr)
1038            end if
1039            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1040 !          work(ir)=rad(ir)*(r(ir)*bes2(arg*rad(ir))*dvwf)
1041            work3(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
1042          end do
1043          call ctrap(mmax,work3,al,result)
1044          ffspl(iq,1,3)=ztor1+result
1045        end do
1046 
1047 !      Compute yp1,ypn for l=2
1048 !      yp1=0 for l=2
1049        yp1=0.0d0
1050 !      ypn=$\int [2 \pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$
1051 !      where x=2 Pi qgrid(mqgrid) r
1052        arg=two_pi*qgrid(mqgrid)
1053        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1054        qr=arg*rad(1)
1055        if(qr<1.d-3)then
1056          bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0
1057        else
1058          bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr
1059        end if
1060 !      ztor1=( (two_pi*rad(1)**2)*dvwf* (bes1(arg*rad(1))-
1061 !      3.0d0*bes2(arg*rad(1))/(arg*rad(1))) ) * rad(1)/7.0d0
1062        ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/7.0d0
1063        do ir=1,mmax
1064          qr=arg*rad(ir)
1065          if(qr<1.d-3)then
1066            bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0
1067          else
1068            bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr
1069          end if
1070          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1071 !        work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
1072 !        (bes1(arg*rad(ir))-3.d0*bes2(arg*rad(ir))/(arg*rad(ir))) )
1073          work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel)
1074        end do
1075        call ctrap(mmax,work3,al,result)
1076        ypn=ztor1+result
1077 
1078 !      Fit spline for l=2 Kleinman-Bylander form factor
1079        call spline(qgrid,ffspl(1,1,3),mqgrid,yp1,ypn,ffspl(1,2,3))
1080 
1081      else
1082 !      or else put form factor to 0 for l=1
1083        ffspl(:,:,3)=0.0d0
1084      end if
1085 !    Endif condition of lmax>1
1086    end if
1087 
1088 !  Finished if lmax=2 (highest nonlocal correction)
1089 !  Do l=3 nonlocal form factor if eta(4) not 0 and lmax>=3
1090    if (lmax>2)then
1091      if(ekb(4)/=0.0d0) then
1092 
1093        lp1=4
1094 !      do q=0 separately; f_3(q=0) vanishes
1095        ffspl(1,1,4)=0.0d0
1096 
1097 !      do rest of q points
1098        do iq=2,mqgrid
1099          arg=two_pi*qgrid(iq)
1100          dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1101          qr=arg*rad(1)
1102          if(qr<1.d-3)then
1103            bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0
1104          else
1105            bessel=bes3_psp5(qr)
1106          end if
1107 !        ztor1=(bes3(arg*rad(1))*rad(1)*dvwf)*rad(1)/9.0d0
1108          ztor1=(bessel*rad(1)*dvwf)*rad(1)/9.0d0
1109          do ir=1,mmax
1110            qr=arg*rad(ir)
1111            if(qr<1.d-3)then
1112              bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0
1113            else
1114              bessel=bes3_psp5(qr)
1115            end if
1116            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1117 !          work(ir)=rad(ir)*(rad(ir)*bes3(arg*rad(ir))*dvwf)
1118            work4(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
1119          end do
1120          call ctrap(mmax,work4,al,result)
1121          ffspl(iq,1,4)=ztor1+result
1122        end do
1123 
1124 !      Compute yp1,ypn for l=3
1125 !      yp1=0 for l=3
1126        yp1=0.0d0
1127 !      ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$
1128 !      where x=2 Pi qgrid(mqgrid) r
1129        arg=two_pi*qgrid(mqgrid)
1130        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1131        qr=arg*rad(1)
1132        if(qr<1.d-3)then
1133          bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0
1134        else
1135          bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr
1136        end if
1137 !      ztor1=( (two_pi*rad(1)**2)*dvwf* (bes2(arg*rad(1))-
1138 !      3.0d0*bes3(arg*rad(1))/(arg*rad(1))) ) * rad(1)/9.0d0
1139        ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/9.0d0
1140        do ir=1,mmax
1141          qr=arg*rad(ir)
1142          if(qr<1.d-3)then
1143            bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0
1144          else
1145            bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr
1146          end if
1147          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1148 !        work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
1149 !        (bes2(arg*rad(ir))-4.d0*bes3(arg*rad(ir))/(arg*rad(ir))) )
1150          work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel)
1151        end do
1152        call ctrap(mmax,work4,al,result)
1153        ypn=ztor1+result
1154 
1155 !      Fit spline for l=3 Kleinman-Bylander form factor
1156        call spline(qgrid,ffspl(1,1,4),mqgrid,yp1,ypn,ffspl(1,2,4))
1157 
1158      else
1159 !      or else put form factor to 0 for l=3
1160        ffspl(:,:,4)=0.0d0
1161      end if
1162 !    Endif condition of lmax>2
1163    end if
1164 
1165 !  Endif condition lmax/=-1
1166  end if
1167 
1168 !DEBUG
1169 !write(std_out,*) 'EKB=',(ekb(iq),iq=1,3)
1170 !write(std_out,*) 'COSKB=',(ckb(iq),iq=1,3)
1171 !ENDDEBUG
1172 
1173  ABI_DEALLOCATE(work1)
1174  ABI_DEALLOCATE(work2)
1175  ABI_DEALLOCATE(work3)
1176  ABI_DEALLOCATE(work4)
1177 
1178  contains
1179 
1180    function  bes0_psp5(arg)
1181 
1182 
1183 !This section has been created automatically by the script Abilint (TD).
1184 !Do not modify the following lines by hand.
1185 #undef ABI_FUNC
1186 #define ABI_FUNC 'bes0_psp5'
1187 !End of the abilint section
1188 
1189    real(dp) :: bes0_psp5,arg
1190    bes0_psp5=sin(arg)/arg
1191  end function bes0_psp5
1192 
1193    function bes1_psp5(arg)
1194 
1195 
1196 !This section has been created automatically by the script Abilint (TD).
1197 !Do not modify the following lines by hand.
1198 #undef ABI_FUNC
1199 #define ABI_FUNC 'bes1_psp5'
1200 !End of the abilint section
1201 
1202    real(dp) :: bes1_psp5,arg
1203    bes1_psp5=(sin(arg)-arg*cos(arg))/arg**2
1204  end function bes1_psp5
1205 
1206    function bes2_psp5(arg)
1207 
1208 
1209 !This section has been created automatically by the script Abilint (TD).
1210 !Do not modify the following lines by hand.
1211 #undef ABI_FUNC
1212 #define ABI_FUNC 'bes2_psp5'
1213 !End of the abilint section
1214 
1215    real(dp) :: bes2_psp5,arg
1216    bes2_psp5=( (3.0d0-arg**2)*sin(arg)-&
1217 &   3.0d0*arg*cos(arg) )      /arg**3
1218  end function bes2_psp5
1219 
1220    function bes3_psp5(arg)
1221 
1222 
1223 !This section has been created automatically by the script Abilint (TD).
1224 !Do not modify the following lines by hand.
1225 #undef ABI_FUNC
1226 #define ABI_FUNC 'bes3_psp5'
1227 !End of the abilint section
1228 
1229    real(dp) :: bes3_psp5, arg
1230    bes3_psp5=(15.d0*sin(arg)-15.d0*arg*cos(arg) &
1231 &   -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4
1232  end function bes3_psp5
1233 
1234 end subroutine psp5nl

ABINIT/psp8lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp8lo

FUNCTION

 Compute sine transform to transform from V(r) to q^2 V(q).
 Computes integrals on linear grid interpolated from the linear input
 grid with a spacing adjusted to ensure convergence at the maximum
 wavevector using corrected trapezoidal integration.

INPUTS

  amesh=spacing for linear radial atomic grid.
  mmax=number of radial r grid points
  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=q grid values (bohr**-1).
  rad(mmax)=r grid values (bohr).
  vloc(mmax)=V(r) on radial grid.
  zion=nominal valence charge of atom.

OUTPUT

  epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
  q2vq(mqgrid)
   =q^2 V(q)
   = -\frac{Zv}{\pi}
     + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
  yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax
   (needed for spline fitter).

PARENTS

      psp8in,psp9in

CHILDREN

      ctrap,splfit,spline

SOURCE

1275 subroutine psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion)
1276 
1277 
1278 !This section has been created automatically by the script Abilint (TD).
1279 !Do not modify the following lines by hand.
1280 #undef ABI_FUNC
1281 #define ABI_FUNC 'psp8lo'
1282 !End of the abilint section
1283 
1284  implicit none
1285 
1286 !Arguments----------------------------------------------------------
1287 !scalars
1288  integer,intent(in) :: mmax,mqgrid
1289  real(dp),intent(in) :: amesh,zion
1290  real(dp),intent(out) :: epsatm,yp1,ypn
1291 !arrays
1292  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
1293  real(dp),intent(out) :: q2vq(mqgrid)
1294 
1295 !Local variables-------------------------------
1296 !Following parameter controls accuracy of Fourier transform based on qmax
1297 !and represents the minimun number of integration points in one period.
1298 !scalars
1299  integer,parameter :: NPT_IN_2PI=200
1300  integer :: ider,iq,ir,irmu,irn,mesh_mult,mmax_new
1301  real(dp) :: amesh_new,arg,fp1,fpn,qmesh,result,ztor1
1302 !arrays
1303  real(dp),allocatable :: rad_new(:),rvlpz(:),rvlpz_new(:),sprvlpz(:,:),work(:)
1304 
1305 ! *************************************************************************
1306 
1307  ABI_ALLOCATE(work,(mmax))
1308  ABI_ALLOCATE(rvlpz,(mmax))
1309 
1310 !Do q=0 separately (compute epsatm)
1311  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
1312 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
1313  do ir=1,mmax
1314    rvlpz(ir)=rad(ir)*vloc(ir)+zion
1315    work(ir)=rad(ir)*rvlpz(ir)
1316  end do
1317 
1318 !Do integral from zero to r(max)
1319  call ctrap(mmax,work,amesh,result)
1320 
1321  epsatm=4.d0*pi*result
1322  q2vq(1)=-zion/pi
1323 
1324 !Find r mesh spacing necessary for accurate integration at qmax
1325  amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1326 
1327 !Choose submultiple of input mesh
1328  mesh_mult=int(amesh/amesh_new) + 1
1329  mmax_new=mesh_mult*(mmax-1)+1
1330  amesh_new=amesh/dble(mesh_mult)
1331 
1332  ABI_ALLOCATE(rad_new,(mmax_new))
1333  ABI_ALLOCATE(rvlpz_new,(mmax_new))
1334 
1335  if(mesh_mult==1) then
1336    rad_new(:)=rad(:)
1337    rvlpz_new(:)=rvlpz(:)
1338  else
1339 !  Set up spline and interpolate to finer mesh.
1340 !  First, compute derivatives at end points
1341    fp1=(-50.d0*rvlpz(1)+96.d0*rvlpz(2)-72.d0*rvlpz(3)+32.d0*rvlpz(4)&
1342 &   -6.d0*rvlpz(5))/(24.d0*amesh)
1343    fpn=(6.d0*rvlpz(mmax-4)-32.d0*rvlpz(mmax-3)+72.d0*rvlpz(mmax-2)&
1344 &   -96.d0*rvlpz(mmax-1)+50.d0*rvlpz(mmax))/(24.d0*amesh)
1345    ABI_ALLOCATE(sprvlpz,(mmax,2))
1346    work(:)=zero
1347 
1348 !  Spline fit
1349    call spline(rad, rvlpz,mmax,fp1,fpn,sprvlpz(:,2))
1350    sprvlpz(:,1)=rvlpz(:)
1351 
1352 !  Set up new radial mesh
1353    irn=1
1354    do ir=1,mmax-1
1355      do irmu=0,mesh_mult-1
1356        rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
1357        irn=irn+1
1358      end do
1359    end do
1360    rad_new(mmax_new)=rad(mmax)
1361 
1362    ider=0
1363    call splfit(rad,work,sprvlpz,ider,rad_new,rvlpz_new,mmax,mmax_new)
1364 
1365    ABI_DEALLOCATE(sprvlpz)
1366    ABI_DEALLOCATE(work)
1367    ABI_ALLOCATE(work,(mmax_new))
1368  end if
1369 
1370 !Loop over q values
1371  do iq=2,mqgrid
1372    arg=2.d0*pi*qgrid(iq)
1373 
1374 !  Set up integrand
1375    do  ir=1,mmax_new
1376      work(ir)=sin(arg*rad_new(ir))*rvlpz_new(ir)
1377    end do
1378 
1379 !  Do integral from zero to rad(mmax)
1380    call ctrap(mmax_new,work,amesh_new,result)
1381 
1382 !  Store q^2 v(q)
1383    q2vq(iq)=q2vq(1)+2.d0*qgrid(iq)*result
1384 
1385  end do
1386 
1387 !Compute derivatives of q^2 v(q) at ends of interval
1388  qmesh=qgrid(2)-qgrid(1)
1389  yp1=(-50.d0*q2vq(1)+96.d0*q2vq(2)-72.d0*q2vq(3)+32.d0*q2vq(4)&
1390 & -6.d0*q2vq(5))/(24.d0*qmesh)
1391  ypn=(6.d0*q2vq(mqgrid-4)-32.d0*q2vq(mqgrid-3)+72.d0*q2vq(mqgrid-2)&
1392 & -96.d0*q2vq(mqgrid-1)+50.d0*q2vq(mqgrid))/(24.d0*qmesh)
1393 
1394  ABI_DEALLOCATE(work)
1395  ABI_DEALLOCATE(rad_new)
1396  ABI_DEALLOCATE(rvlpz_new)
1397  ABI_DEALLOCATE(rvlpz)
1398 
1399 end subroutine psp8lo

ABINIT/psp8nl [ Functions ]

[ Top ] [ Functions ]

NAME

 psp8nl

FUNCTION

 Make Kleinman-Bylander/Bloechl form factors f_ln(q) for each
  projector n for each angular momentum l excepting an l corresponding
  to the local potential.
 Note that an arbitrary local potential can be used, so all l from
  0 to lmax may be represented.

INPUTS

  amesh=grid spacing for uniform (linear) radial grid
  indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln  (if useylm=0)
                                           or i=lmn (if useylm=1)
  lmax=maximum ang momentum for which nonlocal form factor is desired.
    lmax <= 2 allowed.
  lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
        =if useylm=0, max number of (l,n)   comp. over all type of psps
  lnmax=max. number of (l,n) components over all type of psps
  mmax=number of radial grid points for atomic grid
  mqgrid=number of grid points for q grid
  pspso=spin-orbit characteristics, govern the content of ffspl and ekb
   if =0 : this input requires NO spin-orbit characteristics of the psp
   if =2 : this input requires HGH or psp8 characteristics of the psp
   if =3 : this input requires HFN characteristics of the psp
  qgrid(mqgrid)=values at which form factors are returned
  rad(mmax)=radial grid values
  vpspll(mmax,lnmax)=nonlocal projectors for each (l,n) on linear
   radial grid.  Here, these are the  product of the reference
   wave functions and (v(l,n)-vloc), calculated in the psp generation
   program and normalized so that integral(0,rc(l)) vpsll^2 dr = 1,
   which leads to the the usual convention for the energies ekb(l,n)
   also calculated in the psp generation program.

OUTPUT

  ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_ln(q) and
   second derivative from spline fit for each (l,n).

NOTES

 u_l(r) is reference state wavefunction (input as wf);
 j_l(q) is a spherical Bessel function;
 dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
 f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
 where dvms = $\int_0^{rmax} [(u_l(r) dV_l(r))^2 dr]$ is the mean
 square value of the nonlocal correction for angular momentum l.
 Xavier Gonze s E_KB = $ dvms/\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]$.
 This is the eigenvalue of the Kleinman-Bylander operator and sets
 the energy scale of the nonlocal psp corrections.

PARENTS

      psp8in,psp9in

CHILDREN

      ctrap,sbf8,spline

SOURCE

1460 subroutine psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,mqgrid,qgrid,rad,vpspll)
1461 
1462 
1463 !This section has been created automatically by the script Abilint (TD).
1464 !Do not modify the following lines by hand.
1465 #undef ABI_FUNC
1466 #define ABI_FUNC 'psp8nl'
1467 !End of the abilint section
1468 
1469  implicit none
1470 
1471 !Arguments----------------------------------------------------------
1472 !scalars
1473  integer,intent(in) :: lmax,lmnmax,lnmax,mmax,mqgrid
1474  real(dp),intent(in) :: amesh
1475 !arrays
1476  integer,intent(in) :: indlmn(6,lmnmax)
1477  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vpspll(mmax,lnmax)
1478  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
1479 
1480 !Local variables-------------------------------
1481 !Following parameter controls accuracy of Fourier transform based on qmax
1482 !and represents the minimun number of integration points in one period.
1483 !scalars
1484  integer,parameter :: NPT_IN_2PI=200
1485  integer :: iln,iln0,ilmn,iq,ir,irmu,irn,ll,mesh_mult,mmax_new,mvpspll
1486  real(dp) :: amesh_new,arg,c1,c2,c3,c4,dri,qmesh,result,tv,xp,xpm1,xpm2,xpp1
1487  real(dp) :: yp1,ypn
1488 !arrays
1489  real(dp) :: sb_out(4)
1490  real(dp),allocatable :: rad_new(:),vpspll_new(:,:),work(:,:),work2(:)
1491 
1492 ! *************************************************************************
1493 
1494 !Find r mesh spacing necessary for accurate integration at qmax
1495  amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1496 
1497 !Choose submultiple of input mesh
1498  mesh_mult=int(amesh/amesh_new) + 1
1499  mmax_new=mesh_mult*(mmax-1)+1
1500  amesh_new=amesh/dble(mesh_mult)
1501 
1502  ABI_ALLOCATE(rad_new,(mmax_new))
1503  ABI_ALLOCATE(vpspll_new,(mmax_new,lnmax))
1504 
1505  if(mesh_mult==1) then
1506    rad_new(:)=rad(:)
1507  else
1508 !  Set up new radial mesh
1509    irn=1
1510    do ir=1,mmax-1
1511      do irmu=0,mesh_mult-1
1512        rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
1513        irn=irn+1
1514      end do
1515    end do
1516    rad_new(mmax_new)=rad(mmax)
1517  end if
1518 
1519 !Interpolate projectors onto new grid if called for
1520 !Cubic polynomial interpolation is used which is consistent
1521 !with the original interpolation of these functions from
1522 !a log grid to the input linear grid.
1523  dri = one/amesh
1524  do irn=1,mmax_new
1525 !  index to find bracketing input mesh points
1526    if(mesh_mult>1) then
1527      ir = irn/mesh_mult + 1
1528      ir = max(ir,2)
1529      ir = min(ir,mmax-2)
1530 !    interpolation coefficients
1531      xp = dri * (rad_new(irn) - rad(ir))
1532      xpp1 = xp + one
1533      xpm1 = xp - one
1534      xpm2 = xp - two
1535      c1 = -xp * xpm1 * xpm2 * sixth
1536      c2 = xpp1 * xpm1 * xpm2 * half
1537      c3 = - xp * xpp1 * xpm2 * half
1538      c4 = xp * xpp1 * xpm1 * sixth
1539 !    Now do the interpolation on all projectors for this grid point
1540 
1541      iln0=0
1542      do ilmn=1,lmnmax
1543        iln=indlmn(5,ilmn)
1544        if (iln>iln0) then
1545          iln0=iln
1546          tv =  c1 * vpspll(ir - 1, iln) &
1547 &         + c2 * vpspll(ir    , iln) &
1548 &         + c3 * vpspll(ir + 1, iln) &
1549 &         + c4 * vpspll(ir + 2, iln)
1550          if(abs(tv)>tol10) then
1551            vpspll_new(irn,iln)=tv
1552            mvpspll=irn
1553          else
1554            vpspll_new(irn,iln)=zero
1555          end if
1556        end if
1557      end do
1558 
1559    else
1560 !    With no mesh multiplication, just copy projectors
1561      ir=irn
1562      iln0=0
1563      do ilmn=1,lmnmax
1564        iln=indlmn(5,ilmn)
1565        if (iln>iln0) then
1566          iln0=iln
1567          tv = vpspll(ir,iln)
1568          if(abs(tv)>tol10) then
1569            vpspll_new(irn,iln)=tv
1570            mvpspll=irn
1571          else
1572            vpspll_new(irn,iln)=zero
1573          end if
1574        end if
1575      end do
1576 
1577    end if
1578  end do !irn
1579 
1580  ABI_ALLOCATE(work,(mvpspll,lnmax))
1581 
1582 !Loop over q values
1583  do iq=1,mqgrid
1584    arg=2.d0*pi*qgrid(iq)
1585 
1586 !  Set up integrands
1587    do  ir=1,mvpspll
1588      call sbf8(lmax+1,arg*rad_new(ir),sb_out)
1589      iln0=0
1590      do ilmn=1,lmnmax
1591        iln=indlmn(5,ilmn)
1592        if (iln>iln0) then
1593          iln0=iln
1594          ll=indlmn(1,ilmn)
1595          work(ir,iln)=sb_out(ll+1)*vpspll_new(ir,iln)*rad_new(ir)
1596        end if
1597      end do
1598    end do !ir
1599 
1600 !  Do integral from zero to rad_new(mvpspll)
1601    iln0=0
1602    do ilmn=1,lmnmax
1603      iln=indlmn(5,ilmn)
1604      if (iln>iln0) then
1605        iln0=iln
1606        call ctrap(mvpspll,work(1,iln),amesh_new,result)
1607        ffspl(iq,1,iln)=result
1608      end if
1609    end do
1610 
1611 !  End loop over q mesh
1612  end do !iq
1613 
1614 !Fit splines for form factors
1615  ABI_ALLOCATE(work2,(mqgrid))
1616  qmesh=qgrid(2)-qgrid(1)
1617 
1618  iln0=0
1619  do ilmn=1,lmnmax
1620    iln=indlmn(5,ilmn)
1621    if (iln>iln0) then
1622      iln0=iln
1623 !    Compute derivatives of form factors at ends of interval
1624      yp1=(-50.d0*ffspl(1,1,iln)+96.d0*ffspl(2,1,iln)-72.d0*ffspl(3,1,iln)&
1625 &     +32.d0*ffspl(4,1,iln)- 6.d0*ffspl(5,1,iln))/(24.d0*qmesh)
1626      ypn=(6.d0*ffspl(mqgrid-4,1,iln)-32.d0*ffspl(mqgrid-3,1,iln)&
1627 &     +72.d0*ffspl(mqgrid-2,1,iln)-96.d0*ffspl(mqgrid-1,1,iln)&
1628 &     +50.d0*ffspl(mqgrid,1,iln))/(24.d0*qmesh)
1629 
1630      call spline(qgrid,ffspl(1,1,iln),mqgrid,yp1,ypn,ffspl(1,2,iln))
1631    end if
1632  end do
1633 
1634  ABI_DEALLOCATE(rad_new)
1635  ABI_DEALLOCATE(vpspll_new)
1636  ABI_DEALLOCATE(work)
1637  ABI_DEALLOCATE(work2)
1638 
1639 end subroutine psp8nl