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

NOTES

 Test version by DRH - requires very smooth model core charge

SOURCE

1506 subroutine cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d)
1507 
1508 !Arguments ------------------------------------
1509 ! scalars
1510  integer,intent(in) :: mmax,n1xccc
1511  real(dp),intent(in) :: rchrg
1512 !arrays
1513  real(dp),intent(in) :: rad(mmax),ff(mmax),ff1(mmax),ff2(mmax)
1514  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
1515 
1516 !Local variables-------------------------------
1517 ! scalars
1518  integer :: i1xccc
1519  real(dp) :: der1,dern
1520 !arrays
1521  real(dp),allocatable :: ff3(:),ff4(:),gg(:),gg1(:),gg2(:)
1522  real(dp),allocatable :: gg3(:),gg4(:),work(:),xx(:)
1523 
1524 ! *************************************************************************
1525 
1526  !write(std_out,*) 'cc_derivatives : enter'
1527 
1528  ABI_MALLOC(ff3, (mmax))
1529  ABI_MALLOC(ff4, (mmax))
1530  ABI_MALLOC(gg, (n1xccc))
1531  ABI_MALLOC(gg1, (n1xccc))
1532  ABI_MALLOC(gg2, (n1xccc))
1533  ABI_MALLOC(gg3, (n1xccc))
1534  ABI_MALLOC(gg4, (n1xccc))
1535  ABI_MALLOC(work, (mmax))
1536  ABI_MALLOC(xx, (n1xccc))
1537 
1538  ! calculate third derivative ff3 on logarithmic grid
1539  der1=ff2(1)
1540  dern=ff2(mmax)
1541  call spline(rad,ff1,mmax,der1,dern,ff3)
1542 
1543  ! calculate fourth derivative ff4 on logarithmic grid
1544  der1=0.d0
1545  dern=0.d0
1546  call spline(rad,ff2,mmax,der1,dern,ff4)
1547 
1548  ! generate uniform mesh xx in the box cut by rchrg:
1549  do i1xccc=1,n1xccc
1550    xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1)
1551  end do
1552 
1553  !now interpolate core charge and derivatives on the uniform grid
1554  !core charge, input=ff,  output=gg
1555  call splint(mmax,rad,ff,ff2,n1xccc,xx,gg)
1556 
1557  ! first derivative input=ff1, output=gg1
1558  call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1)
1559 
1560  !normalize gg1
1561  !gg1(:)=gg1(:)*rchrg
1562 
1563  ! second derivative input=ff2, output=gg2
1564  call splint(mmax,rad,ff2,ff4,n1xccc,xx,gg2)
1565 
1566  !normalize gg2
1567  !gg2(:)=gg2(:)*rchrg**2
1568 
1569  ! reallocate work otherwise the calls to spline crash (n1xccc /= mmax)
1570  ABI_FREE(work)
1571  ABI_MALLOC(work, (n1xccc))
1572 
1573 !recalculate 3rd derivative consistent with spline fit to first derivative on linear grid
1574  der1=gg2(1)
1575  dern=gg2(n1xccc)
1576  call spline(xx,gg1,n1xccc,der1,dern,gg3)
1577 
1578 !calculate 4th derivative consistent with spline fit to second derivative on linear grid
1579  der1=0.0d0
1580  dern=0.0d0
1581  call spline(xx,gg2,n1xccc,der1,dern,gg4)
1582 
1583 !now calculate second to fourth derivative by forward differences
1584 !to avoid numerical noise uses a smoothing function
1585 !
1586 !call smooth(gg1,n1xccc,10)
1587 
1588 !gg2(n1xccc)=0.0
1589 !do i1xccc=1,n1xccc-1
1590 !gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1)
1591 !end do
1592 
1593 !call smooth(gg2,n1xccc,10)
1594 
1595 !gg3(n1xccc)=0.0
1596 !do i1xccc=1,n1xccc-1
1597 !gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1)
1598 !end do
1599 
1600 !call smooth(gg3,n1xccc,10)
1601 
1602 !gg4(n1xccc)=0.0
1603 !do i1xccc=1,n1xccc-1
1604 !gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1)
1605 !end do
1606 
1607 !call smooth(gg4,n1xccc,10)
1608 
1609 !write on xcc1d
1610 !normalize to unit range usage later in program
1611  xccc1d(:,1)=gg(:)
1612  xccc1d(:,2)=gg1(:)*rchrg
1613  xccc1d(:,3)=gg2(:)*rchrg**2
1614  xccc1d(:,4)=gg3(:)*rchrg**3
1615  xccc1d(:,5)=gg4(:)*rchrg**4
1616 !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc
1617 
1618 !DEBUG
1619 !note: the normalization condition is the following:
1620 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg
1621 !
1622 !norm=0.d0
1623 !do i1xccc=1,n1xccc
1624 !norm = norm + 4.d0*pi*rchrg/dble(n1xccc-1)*&
1625 !&             xx(i1xccc)**2*xccc1d(i1xccc,1)
1626 !end do
1627 !write(std_out,*) ' norm=',norm
1628 !
1629 !write(std_out,*)' psp6cc_drh : output of core charge density and derivatives '
1630 !write(std_out,*)'   xx          gg           gg1  '
1631 !do i1xccc=1,n1xccc
1632 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
1633 !end do
1634 !write(std_out,*)'   xx          gg2          gg3  '
1635 !do i1xccc=1,n1xccc
1636 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
1637 !end do
1638 !write(std_out,*)'   xx          gg4          gg5  '
1639 !do i1xccc=1,n1xccc
1640 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
1641 !end do
1642 !write(std_out,*)' psp1cc : debug done, stop '
1643 !stop
1644 !ENDDEBUG
1645 
1646  ABI_FREE(ff3)
1647  ABI_FREE(ff4)
1648  ABI_FREE(gg)
1649  ABI_FREE(gg1)
1650  ABI_FREE(gg2)
1651  ABI_FREE(gg3)
1652  ABI_FREE(gg4)
1653  ABI_FREE(work)
1654  ABI_FREE(xx)
1655 
1656 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)

SOURCE

213 subroutine gg1cc(gg1cc_xx,xx)
214 
215 !Arguments ------------------------------------
216 !scalars
217  real(dp),intent(in) :: xx
218  real(dp),intent(out) :: gg1cc_xx
219 
220 !Local variables -------------------------------------------
221 !The c s are coefficients for Taylor expansion of the analytic form near xx=0, 1/2, and 1.
222 !scalars
223  real(dp) :: c21=4.d0/9.d0,c22=-40.d0/27.d0,c23=20.d0/3.d0-16.d0*pi**2/27.d0
224  real(dp) :: c24=-4160.d0/243.d0+160.d0*pi**2/81.d0,c31=1.d0/36.d0
225  real(dp) :: c32=-25.d0/108.d0,c33=485.d0/432.d0-pi**2/27.d0
226  real(dp) :: c34=-4055.d0/972.d0+25.d0*pi**2/81.d0
227 
228 ! *************************************************************************
229 
230 !Cut off beyond 3/gcut=xcccrc
231  if (xx>3.0d0) then
232    gg1cc_xx=0.0d0
233 !  Take care of difficult limits near x=0, 1/2, and 1
234  else if (abs(xx)<=1.d-09) then
235    gg1cc_xx=1.d0
236  else if (abs(xx-0.5d0)<=1.d-04) then
237 !  (this limit and next are more troublesome for numerical cancellation)
238    gg1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*c24))
239  else if (abs(xx-1.d0)<=1.d-04) then
240    gg1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*c34))
241  else
242 !  The following is the square of the Fourier transform of a
243 !  function built out of two spherical bessel functions in G
244 !  space and cut off absolutely beyond gcut
245    gg1cc_xx=(sin(2.0d0*pi*xx)/( (2.0d0*pi*xx) * &
246 &   (1.d0-4.0d0*xx**2)*(1.d0-xx**2) )  )**2
247  end if
248 
249 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)}$

SOURCE

273 subroutine gp1cc(gp1cc_xx,xx)
274 
275 !Arguments ------------------------------------
276 !scalars
277  real(dp),intent(in) :: xx
278  real(dp),intent(out) :: gp1cc_xx
279 
280 !Local variables -------------------------------------------
281 !scalars
282  real(dp),parameter :: c11=20.d0-8.d0*pi**2/3.d0
283  real(dp),parameter :: c12=268.d0-160.d0/3.d0*pi**2+128.d0/45.d0*pi**4
284  real(dp),parameter :: c21=-40.d0/27.d0,c22=40.d0/3.d0-32.d0*pi**2/27.d0
285  real(dp),parameter :: c23=-4160.d0/81.d0+160.d0*pi**2/27.d0
286  real(dp),parameter :: c24=157712.d0/729.d0-320.d0*pi**2/9.d0+512.d0*pi**4/405.d0
287  real(dp),parameter :: c25=-452200.d0/729.d0+83200.d0*pi**2/729.d0-1280.d0*pi**4/243.d0
288  real(dp),parameter :: c31=-25.d0/108.d0,c32=485.d0/216.d0-2.d0*pi**2/27.d0
289  real(dp),parameter :: c33=-4055.d0/324.d0+25.d0*pi**2/27.d0
290  real(dp),parameter :: c34=616697.d0/11664.d0-485.d0*pi**2/81.d0+32.d0*pi**4/405.d0
291  real(dp),parameter :: c35=-2933875.d0/15552.d0+20275.d0*pi**2/729.d0-200.d0*pi**4/243.d0
292  real(dp),parameter :: two_pim1=1.0d0/two_pi
293  real(dp) :: denom,phi,phip
294 
295 ! *************************************************************************
296 
297 !Cut off beyond r=3*xcccrc is already done at the calling level
298  if (xx>1.001d0) then
299 !  The part that follows will be repeated later, but written in this way,
300 !  only one "if" condition is tested in most of the cases (1.001 < x < 3.0)
301    denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
302    phi=denom*sin(two_pi*xx)*two_pim1
303    phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
304    gp1cc_xx=2.d0*phi*phip
305 !  Handle limits where denominator vanishes
306  else if (abs(xx)<1.d-03) then
307    gp1cc_xx=xx*(c11+xx**2*c12)
308  else if (abs(xx-0.5d0)<=1.d-03) then
309    gp1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*(c24+(xx-0.5d0)*c25)))
310  else if (abs(xx-1.d0)<=1.d-03) then
311    gp1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*(c34+(xx-1.0d0)*c35)))
312  else
313 !  Here is the repeated part ...
314    denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
315    phi=denom*sin(two_pi*xx)*two_pim1
316    phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
317    gp1cc_xx=2.d0*phi*phip
318  end if
319 
320 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.

SOURCE

339 subroutine gpp1cc(gpp1cc_xx,xx)
340 
341 !Arguments ------------------------------------
342 !scalars
343  real(dp),intent(in) :: xx
344  real(dp),intent(out) :: gpp1cc_xx
345 
346 !Local variables -------------------------------------------
347 !scalars
348  real(dp),parameter :: c1=20.d0-8.d0*pi**2/3.d00
349  real(dp),parameter :: c2=40.d0/3.d0-32.d0*pi**2/27.d0
350  real(dp),parameter :: c3=-8320.d0/81.d0+320.d0*pi**2/27.d0
351  real(dp),parameter :: c4=157712.d0/243.d0-320.d0*pi**2/3.d0+512.d0*pi**4/135.d0
352  real(dp),parameter :: c5=-18088.d2/729.d0+3328.d2*pi**2/729.d0-5120.d0*pi**4/243.d0
353  real(dp),parameter :: c6=485.d0/216.d0-2.d0*pi**2/27.d0
354  real(dp),parameter :: c7=-4055.d0/162.d0+50.d0*pi**2/27.d0
355  real(dp),parameter :: c8=616697.d0/3888.d0-485.d0*pi**2/27.d0+32.d0*pi**4/135.d0
356  real(dp),parameter :: c9=-2933875.d0/3888.d0+81100.d0*pi**2/729.d0-800.d0*pi**4/243.d0
357  real(dp) :: t1,t10,t100,t11,t12,t120,t121,t122,t127,t138,t14,t140,t15,t152
358  real(dp) :: t157,t16,t160,t17,t174,t175,t18,t19,t2,t20,t21,t23,t24,t3,t31,t33
359  real(dp) :: t34,t4,t41,t42,t44,t45,t46,t5,t54,t55,t56,t57,t6,t62,t64,t65,t7
360  real(dp) :: t72,t78,t79,t8,t85,t9,t93
361 
362 ! *************************************************************************
363 
364  if (xx>3.0d0) then
365 !  Cut off beyond 3/gcut=3*xcccrc
366    gpp1cc_xx=0.0d0
367 !  Take care of difficult limits near xx=0, 1/2, and 1
368  else if (abs(xx)<=1.d-09) then
369    gpp1cc_xx=c1
370  else if (abs(xx-0.5d0)<=1.d-04) then
371 !  (this limit and next are more troublesome for numerical cancellation)
372    gpp1cc_xx=c2+(xx-0.5d0)*(c3+(xx-0.5d0)*(c4+(xx-0.5d0)*c5))
373  else if (abs(xx-1.d0)<=1.d-04) then
374    gpp1cc_xx=c6+(xx-1.0d0)*(c7+(xx-1.0d0)*(c8+(xx-1.0d0)*c9))
375  else
376 
377 !  Should fix up this Maple fortran later
378    t1 = xx**2
379    t2 = 1/t1
380    t3 = 1/Pi
381    t4 = 2*xx
382    t5 = t4-1
383    t6 = t5**2
384    t7 = 1/t6
385    t8 = t4+1
386    t9 = t8**2
387    t10 = 1/t9
388    t11 = xx-1
389    t12 = t11**2
390    t14 = 1/t12/t11
391    t15 = xx+1
392    t16 = t15**2
393    t17 = 1/t16
394    t18 = Pi*xx
395    t19 = sin(t18)
396    t20 = cos(t18)
397    t21 = t20**2
398    t23 = t19*t21*t20
399    t24 = t17*t23
400    t31 = t19**2
401    t33 = t31*t19*t20
402    t34 = t17*t33
403    t41 = Pi**2
404    t42 = 1/t41
405    t44 = 1/t16/t15
406    t45 = t31*t21
407    t46 = t44*t45
408    t54 = 1/t1/xx
409    t55 = 1/t12
410    t56 = t55*t46
411    t57 = t10*t56
412    t62 = t9**2
413    t64 = t17*t45
414    t65 = t55*t64
415    t72 = 1/t9/t8
416    t78 = t14*t64
417    t79 = t10*t78
418    t85 = t12**2
419    t93 = t21**2
420    t100 = t31**2
421    t120 = 1/t6/t5
422    t121 = t55*t34
423    t122 = t10*t121
424    t127 = t16**2
425    t138 = t6**2
426    t140 = t10*t65
427    t152 = t72*t65
428    t157 = t7*t140
429    t160 = t1**2
430    t174 = t55*t24
431    t175 = t10*t174
432    gpp1cc_xx = 8*t2*t3*t7*t10*t14*t34+8*t2*t42*t7*t10*t14*t46&
433 &   -8*t2*t3*t7*t10*t14*t24+8*t2*t3*t7*t10*t55*t44*t33+&
434 &   6*t2*t42*t7*t10*t55/t127*t45+24*t2*t42/t138*t140+&
435 &   16*t54*t42*t120*t140+16*t2*t3*t120*t122+16*t2&
436 &   *t42*t7*t72*t78-8*t2*t3*t7*t10*t55*t44*t23-8*t54*t3*t7*t175&
437 &   +2*t2*t7*t10*t55*t17*t100+2*t2*t7*t10*t55*t17*t93+&
438 &   8*t54*t42*t7*t79+16*t2*t42*t7*t72*t56+6*t2*t42*t7*t10/t85&
439 &   *t64+24*t2*t42*t7/t62*t65+8*t54*t42*t7*t57-&
440 &   16*t2*t3*t7*t72*t174+8*t54*t3*t7*t122-16*t2*t3*t120*t175&
441 &   +16*t2*t42*t120*t79+16*t2*t42*t120*t57+16*t54*t42*t7*t152+&
442 &   32*t2*t42*t120*t152+16*t2*t3*t7*t72*t121-12*t2*t157+&
443 &   6/t160*t42*t157
444  end if
445 
446 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-2024 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_psptk
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27  use m_splines
28 
29  use m_numeric_tools,  only : ctrap
30  use m_special_funcs,  only : sbf8
31 
32  implicit none
33 
34  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.

SOURCE

 88 subroutine psp1cc(fchrg,n1xccc,xccc1d)
 89 
 90 !Arguments ------------------------------------
 91 !scalars
 92  integer,intent(in) :: n1xccc
 93  real(dp),intent(in) :: fchrg
 94 !arrays
 95  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
 96 
 97 !Local variables-------------------------------
 98 !scalars
 99  integer :: i1xccc,ider
100  real(dp) :: der1,dern,factor,gg1cc_xx,gp1cc_xx,gpp1cc_xx,xx
101  character(len=500) :: message
102 !arrays
103  real(dp),allocatable :: ff(:),ff2(:),work(:),yy(:)
104 
105 ! *************************************************************************
106 
107  ABI_MALLOC(ff,(n1xccc))
108  ABI_MALLOC(ff2,(n1xccc))
109  ABI_MALLOC(work,(n1xccc))
110  ABI_MALLOC(yy,(n1xccc))
111 
112  if(n1xccc > 1)then
113    factor=one/dble(n1xccc-1)
114    do i1xccc=1,n1xccc
115      yy(i1xccc)=(i1xccc-1)*factor
116    end do
117  else
118    write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc
119    ABI_BUG(message)
120  end if
121 
122 !Initialization, to avoid some problem with some compilers
123  xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
124 
125 !Take care of each derivative separately
126  do ider=0,2
127 
128    if(ider==0)then
129 !    Generate spline fitting for the function gg
130      do i1xccc=1,n1xccc
131        xx=three*yy(i1xccc)
132        call gg1cc(gg1cc_xx,xx)
133        ff(i1xccc)=fchrg*gg1cc_xx
134      end do
135 !    Complete with derivatives at end points
136      der1=zero
137      call gp1cc(gp1cc_xx,three)
138      dern=three*fchrg*gp1cc_xx
139    else if(ider==1)then
140 !    Generate spline fitting for the function gp
141      do i1xccc=1,n1xccc
142        xx=three*yy(i1xccc)
143        call gp1cc(gp1cc_xx,xx)
144        ff(i1xccc)=three*fchrg*gp1cc_xx
145      end do
146 !    Complete with derivatives at end points, already estimated
147      der1=xccc1d(1,ider+2)
148      dern=xccc1d(n1xccc,ider+2)
149    else if(ider==2)then
150 !    Generate spline fitting for the function gpp
151 !    (note : the function gpp has already been estimated, for the spline
152 !    fitting of the function gg, but it is replaced here by the more
153 !    accurate analytic derivative)
154      do i1xccc=1,n1xccc
155        xx=three*yy(i1xccc)
156        call gpp1cc(gpp1cc_xx,xx)
157        ff(i1xccc)=9.0_dp*fchrg*gpp1cc_xx
158      end do
159 !    Complete with derivatives of end points
160      der1=xccc1d(1,ider+2)
161      dern=xccc1d(n1xccc,ider+2)
162    end if
163 
164 !  Produce second derivative numerically, for use with splines
165    call spline(yy,ff,n1xccc,der1,dern,ff2)
166    xccc1d(:,ider+1)=ff(:)
167    xccc1d(:,ider+3)=ff2(:)
168  end do
169 
170  xccc1d(:,6)=zero
171 
172 !DEBUG
173 !write(std_out,*)' psp1cc : output of core charge density and derivatives '
174 !write(std_out,*)'   yy          gg           gp  '
175 !do i1xccc=1,n1xccc
176 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
177 !end do
178 !write(std_out,*)'   yy          gpp          gg2  '
179 !do i1xccc=1,n1xccc
180 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
181 !end do
182 !write(std_out,*)'   yy          gp2          gpp2  '
183 !do i1xccc=1,n1xccc
184 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
185 !end do
186 !write(std_out,*)' psp1cc : debug done, stop '
187 !stop
188 !ENDDEBUG
189 
190  ABI_FREE(ff)
191  ABI_FREE(ff2)
192  ABI_FREE(work)
193  ABI_FREE(yy)
194 
195 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).

SOURCE

480 subroutine psp5lo(al,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
481 &                  vloc,yp1,ypn,zion)
482 
483 !Arguments----------------------------------------------------------
484 !scalars
485  integer,intent(in) :: mmax,mqgrid
486  real(dp),intent(in) :: al,zion
487  real(dp),intent(out) :: epsatm,yp1,ypn
488 !arrays
489  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
490  real(dp),intent(out) :: q2vq(mqgrid)
491 
492 !Local variables-------------------------------
493 !scalars
494  integer :: iq,ir
495  real(dp),parameter :: scale=10.0d0
496  real(dp) :: arg,result,rmtoin,test,ztor1
497 !arrays
498  real(dp),allocatable :: work(:)
499 
500 ! *************************************************************************
501 
502  ABI_MALLOC(work,(mmax))
503 
504 !Do q=0 separately (compute epsatm)
505 !Do integral from 0 to r1
506  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
507 
508 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
509 !with extra factor of r to convert to uniform grid in exponent
510  do ir=1,mmax
511 !  First handle tail region
512    test=vloc(ir)+zion/rad(ir)
513 !  DEBUG
514 !  write(std_out,*)ir,rad(ir),test
515 !  ENDDEBUG
516 !  Ignore small contributions, or impose a cut-off in the case
517 !  the pseudopotential data are in single precision.
518 !  (it is indeed expected that vloc is very close to zero beyond 20,
519 !  so a value larger than 2.0d-8 is considered anomalous)
520    if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
521      work(ir)=zero
522    else
523      work(ir)=(rad(ir)*rad(ir))*(rad(ir)*vloc(ir)+zion)
524    end if
525  end do
526 !DEBUG
527 !write(std_out,*)' psp5lo : stop '
528 !stop
529 !ENDDEBUG
530 
531 !Do integral from r(1) to r(max)
532  call ctrap(mmax,work,al,result)
533 !Do integral from r(mmax) to infinity
534 !compute decay length lambda at r(mmax)
535 !$\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ &
536 !$(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$
537 !rmtoin=$(rad(mmax)*vloc(mmax)+zion)*(rad(mmax)+1.d0/\lambda)/\lambda$
538 !Due to inability to fit exponential decay to r*V(r)+Zv
539 !in tail, NO TAIL CORRECTION IS APPLIED
540 !(numerical trouble might be removed if atomic code is
541 !cleaned up in tail region)
542  rmtoin=0.0d0
543 
544  epsatm=4.d0*pi*(result+ztor1+rmtoin)
545 
546  q2vq(1)=-zion/pi
547 
548 !Loop over q values
549  do iq=2,mqgrid
550    arg=2.d0*pi*qgrid(iq)
551 !  ztor1=$ -Zv/\pi+2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
552    ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion)* &
553 &   cos(arg*rad(1)) )/pi
554 
555 !  set up integrand
556    do  ir=1,mmax
557      test=vloc(ir)+zion/rad(ir)
558 !    Ignore contributions within decade of machine precision
559      if ((scale+abs(test)).eq.scale) then
560        work(ir)=zero
561      else
562        work(ir)=rad(ir)*sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
563      end if
564    end do
565 !  do integral from r(1) to r(mmax)
566    call ctrap(mmax,work,al,result)
567 
568 !  do integral from r(mmax) to infinity
569 !  rmtoin=(r(mmax)*vr(mmax)+zion)*(lambda*sin(arg*r(mmax))+
570 !  arg*cos(arg*r(mmax)))/(arg**2+lambda**2)
571 !  See comment above; no tail correction
572    rmtoin=0.0d0
573 
574 !  store q^2 v(q)
575    q2vq(iq)=ztor1+2.d0*qgrid(iq)*(result+rmtoin)
576 
577  end do
578 
579 !Compute derivatives of q^2 v(q) at ends of interval
580  yp1=0.0d0
581 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
582 !integral from 0 to r1
583  arg=2.0d0*pi*qgrid(mqgrid)
584  ztor1=zion*rad(1)*sin(arg*rad(1))
585  ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + &
586 & (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1))
587 !integral from r(mmax) to infinity is overkill; ignore
588 !set up integrand
589  do ir=1,mmax
590    test=vloc(ir)+zion/rad(ir)
591 !  Ignore contributions within decade of machine precision
592    if ((scale+abs(test)).eq.scale) then
593      work(ir)=0.0d0
594    else
595      work(ir)=rad(ir)*(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * &
596 &     (rad(ir)*vloc(ir)+zion)
597    end if
598  end do
599  call ctrap(mmax,work,al,result)
600  ypn=2.0d0 * (ztor1 + result)
601 
602  ABI_FREE(work)
603 
604 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.

SOURCE

 651 subroutine psp5nl(al,ekb,ffspl,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll)
 652 
 653 !Arguments ------------------------------------
 654 !scalars
 655  real(dp),intent(in) :: al
 656  integer,intent(in) :: lmax,mmax,mpsang,mqgrid
 657 !arrays
 658  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax),vpspll(mmax,mpsang)
 659  real(dp),intent(in) :: wfll(mmax,mpsang)
 660  real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang)
 661 
 662 !Local variables-------------------------------
 663 !scalars
 664  integer,parameter :: dpsang=5
 665  integer :: iq,ir,lp1
 666  real(dp) :: arg,bessel,dvwf,qr,result,yp1,ypn,ztor1
 667  character(len=500) :: message
 668 !arrays
 669  real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang)
 670  real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:)
 671 
 672 !*************************************************************************
 673 
 674 !l=0,1,2 and 3 spherical Bessel functions
 675 !The accuracy of the bes1, bes2, bes3 functions for small arguments
 676 !may be insufficient. In the present version
 677 !of the routines, some care is taken with the value of the argument.
 678 !If smaller than 1.d-3, a two terms
 679 !Taylor series expansion is prefered.
 680 ! bes0(arg)=sin(arg)/arg
 681 ! bes1(arg)=(sin(arg)-arg*cos(arg))/arg**2
 682 ! bes2(arg)=( (3.0d0-arg**2)*sin(arg)-&
 683 !& 3.0d0*arg*cos(arg) )      /arg**3
 684 
 685 ! bes3(arg)=(15.d0*sin(arg)-15.d0*arg*cos(arg) &
 686 !& -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4
 687 
 688 !Zero out Kleinman-Bylander energies ekb
 689  ekb(:)=0.0d0
 690 
 691  ABI_MALLOC(work1,(mmax))
 692  ABI_MALLOC(work2,(mmax))
 693  ABI_MALLOC(work3,(mmax))
 694  ABI_MALLOC(work4,(mmax))
 695 
 696 !Allow for no nonlocal correction (lmax=-1)
 697  if (lmax/=-1) then
 698 
 699 !  Check that lmax is within allowed range
 700    if (lmax<0.or.lmax>3) then
 701      write(message, '(a,i12,a,a,a,a,a,a,a)' )&
 702 &     'lmax=',lmax,' is not an allowed value.',ch10,&
 703 &     'Allowed values are -1 for no nonlocal correction or else',ch10,&
 704 &     '0, 1,2 or 3 for maximum l nonlocal correction.',ch10,&
 705 &     'Action: check the input atomic psp data file for lmax.'
 706      ABI_ERROR(message)
 707    end if
 708 
 709 !  Compute normalizing integrals eta=<dV> and mean square
 710 !  nonlocal psp correction dvms=<dV^2>
 711 !  "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction
 712    do lp1=1,lmax+1
 713 
 714 !    integral from 0 to r1
 715      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)
 716      ztor1=(wfll(1,lp1)*dvwf)*rad(1)/dble(2*(lp1-1)+3)
 717 !    integrand for r1 to r(mmax) (incl extra factor of r)
 718      do ir=1,mmax
 719        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
 720        work1(ir)=rad(ir)*(wfll(ir,lp1)*dvwf)
 721      end do
 722 !    do integral by corrected trapezoidal integration
 723      call ctrap(mmax,work1,al,result)
 724      eta(lp1)=ztor1+result
 725 
 726      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)
 727      ztor1=dvwf**2*rad(1)/dble(2*(lp1-1)+3)
 728      do ir=1,mmax
 729        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
 730        work1(ir)=rad(ir)*(dvwf**2)
 731      end do
 732      call ctrap(mmax,work1,al,result)
 733      dvms(lp1)=ztor1+result
 734 
 735 !    DEBUG
 736 !    Compute the norm of wfll
 737 !    wf=wfll(1,lp1)
 738 !    ztor1=wf**2*rad(1)/dble(2*(lp1-1)+3)
 739 !    do ir=1,mmax
 740 !    wf=wfll(ir,lp1)
 741 !    work1(ir)=rad(ir)*(wf**2)
 742 !    end do
 743 !    call ctrap(mmax,work1,al,result)
 744 !    norm=ztor1+result
 745 !    write(std_out,*)' lp1, norm',lp1,norm
 746 !    ENDDEBUG
 747 
 748 !    If dvms is not 0 for any given angular momentum l,
 749 !    compute Xavier Gonze's definition of the Kleinman-Bylander
 750 !    energy E_KB = dvms/eta.  In this case also renormalize
 751 !    the projection operator to u_KB(r)=$u_l(r)*dV(r)/\sqrt{dvms}$.
 752 !    This means dvwf gets multiplied by the normalization factor
 753 !    "renorm"=$1/\sqrt{dvms}$ as seen below.
 754      if (dvms(lp1)/=0.0d0) then
 755        ekb(lp1)=dvms(lp1)/eta(lp1)
 756        renorm(lp1)=1.0d0/sqrt(dvms(lp1))
 757 !      ckb is Kleinman-Bylander "cosine" (Xavier Gonze)
 758        ckb(lp1)=eta(lp1)/sqrt(dvms(lp1))
 759      else
 760        ekb(lp1)=0.0d0
 761      end if
 762 
 763    end do
 764 
 765 !  l=0 form factor if ekb(1) not 0 (lmax always at least 0)
 766    if (ekb(1)/=0.0d0) then
 767 
 768 !    do q=0 separately
 769      lp1=1
 770 !    0 to r1 integral
 771      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 772      ztor1=(rad(1)*dvwf)*rad(1)/3.0d0
 773 !    integrand
 774      do ir=1,mmax
 775        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 776        work1(ir)=rad(ir)*(rad(ir)*dvwf)
 777      end do
 778      call ctrap(mmax,work1,al,result)
 779      ffspl(1,1,1)=ztor1+result
 780 
 781 !    do rest of q points
 782      do iq=2,mqgrid
 783        arg=two_pi*qgrid(iq)
 784 !      0 to r1 integral
 785        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 786        ztor1=(bes0_psp5(arg*rad(1))*rad(1)*dvwf)*rad(1)/3.0d0
 787 !      integrand
 788        do ir=1,mmax
 789          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 790          work1(ir)=rad(ir)*(rad(ir)*bes0_psp5(arg*rad(ir))*dvwf)
 791        end do
 792        call ctrap(mmax,work1,al,result)
 793        ffspl(iq,1,1)=ztor1+result
 794      end do
 795 
 796 !    Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
 797 !    yp1=0 for l=0
 798      yp1=0.0d0
 799 !    ypn=$ \int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$
 800      arg=two_pi*qgrid(mqgrid)
 801      dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 802      qr=arg*rad(1)
 803      if(qr<1.d-3)then
 804        bessel=(10.d0-qr*qr)*qr/30.0d0
 805      else
 806        bessel=bes1_psp5(qr)
 807      end if
 808 !    ztor1=(-bes1(arg*rad(1))*two_pi*rad(1)*r(1)*dvwf)*rad(1)/5.0d0
 809      ztor1=(-bessel*two_pi*rad(1)*rad(1)*dvwf)*rad(1)/5.0d0
 810      do ir=1,mmax
 811        qr=arg*rad(ir)
 812        if(qr<1.d-3)then
 813          bessel=(10.d0-qr*qr)*qr/30.0d0
 814        else
 815          bessel=bes1_psp5(qr)
 816        end if
 817        dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 818 !      work(ir)=rad(ir)*(-bes1(arg*rad(ir))*two_pi*rad(ir)*rad(ir)*dvwf)
 819        work1(ir)=rad(ir)*(-bessel*two_pi*rad(ir)*rad(ir)*dvwf)
 820      end do
 821      call ctrap(mmax,work1,al,result)
 822      ypn=ztor1+result
 823 
 824 !    Fit spline to get second derivatives by spline fit
 825      call spline(qgrid,ffspl(1,1,1),mqgrid,yp1,ypn,ffspl(1,2,1))
 826 
 827    else
 828 !    or else put nonlocal correction at l=0 to 0
 829      ffspl(:,:,1)=0.0d0
 830    end if
 831 
 832 !  Finished if lmax=0 (highest nonlocal correction)
 833 !  Do l=1 form factor if ekb(2) not 0 and lmax>=1
 834    if (lmax>0)then
 835      if(ekb(2)/=0.0d0) then
 836 
 837        lp1=2
 838 !      do q=0 separately: f_1(q=0) vanishes !
 839        ffspl(1,1,2)=0.0d0
 840 
 841 !      do rest of q points
 842        do iq=2,mqgrid
 843          arg=two_pi*qgrid(iq)
 844          dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 845          qr=arg*rad(1)
 846          if(qr<1.d-3)then
 847            bessel=(10.d0-qr*qr)*qr/30.0d0
 848          else
 849            bessel=bes1_psp5(qr)
 850          end if
 851 !        ztor1=(bes1(arg*rad(1))*rad(1)*dvwf)*rad(1)/5.0d0
 852          ztor1=(bessel*rad(1)*dvwf)*rad(1)/5.0d0
 853 
 854          do ir=1,mmax
 855            qr=arg*rad(ir)
 856            if(qr<1.d-3)then
 857              bessel=(10.d0-qr*qr)*qr/30.0d0
 858            else
 859              bessel=bes1_psp5(qr)
 860            end if
 861            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 862            work2(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
 863          end do
 864 
 865          call ctrap(mmax,work2,al,result)
 866          ffspl(iq,1,2)=ztor1+result
 867        end do
 868 
 869 !      Compute yp1,ypn for l=1
 870 !      yp1=$\displaystyle \int [2\pi r^2 wf(r) dV(r)]/3$
 871        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 872        ztor1=((two_pi*rad(1)**2)*dvwf)*rad(1)/(3.0d0*5.0d0)
 873        do ir=1,mmax
 874          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 875          work2(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf/3.0d0)
 876        end do
 877        call ctrap(mmax,work2,al,result)
 878        yp1=ztor1+result
 879 !      ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$
 880 !      where x=2 Pi qgrid(mqgrid) r
 881        arg=two_pi*qgrid(mqgrid)
 882        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 883        qr=arg*rad(1)
 884        if(qr<1.d-3)then
 885          bessel=(10.d0-3.0d0*qr*qr)/30.0d0
 886        else
 887          bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr
 888        end if
 889 !      ztor1=( (two_pi*rad(1)**2)*dvwf* (bes0(arg*rad(1))-
 890 !      2.0d0*bes1(arg*rad(1))/(arg*rad(1))) ) * rad(1)/5.0d0
 891        ztor1=( (two_pi*rad(1)**2)*dvwf*bessel)*  rad(1)/5.0d0
 892 
 893        do ir=1,mmax
 894          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 895          qr=arg*rad(ir)
 896          if(qr<1.d-3)then
 897            bessel=(10.d0-3.0d0*qr*qr)/30.0d0
 898          else
 899            bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr
 900          end if
 901 !        work(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
 902 !        (bes0(arg*rad(ir))-2.d0*bes1(arg*rad(ir))/(arg*rad(ir))) )
 903          work2(ir)=rad(ir)*(two_pi*rad(ir)**2)*dvwf*bessel
 904        end do
 905        call ctrap(mmax,work2,al,result)
 906        ypn=ztor1+result
 907 
 908 !      Fit spline for l=1 Kleinman-Bylander form factor
 909        call spline(qgrid,ffspl(1,1,2),mqgrid,yp1,ypn,ffspl(1,2,2))
 910 
 911      else
 912 !      or else put form factor to 0 for l=1
 913        ffspl(:,:,2)=0.0d0
 914      end if
 915 !    Endif condition of lmax>0
 916    end if
 917 
 918 !  Finished if lmax=1 (highest nonlocal correction)
 919 !  Do l=2 nonlocal form factor if eta(3) not 0 and lmax>=2
 920    if (lmax>1)then
 921      if(ekb(3)/=0.0d0) then
 922 
 923        lp1=3
 924 !      do q=0 separately; f_2(q=0) vanishes
 925        ffspl(1,1,3)=0.0d0
 926 
 927 !      do rest of q points
 928        do iq=2,mqgrid
 929          arg=two_pi*qgrid(iq)
 930          dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 931          qr=arg*rad(1)
 932          if(qr<1.d-3)then
 933            bessel=qr*qr/15.0d0-qr**4/210.0d0
 934          else
 935            bessel=bes2_psp5(qr)
 936          end if
 937 !        ztor1=(bes2(arg*rad(1))*rad(1)*dvwf)*rad(1)/7.0d0
 938          ztor1=(bessel*rad(1)*dvwf)*rad(1)/7.0d0
 939          do ir=1,mmax
 940            qr=arg*rad(ir)
 941            if(qr<1.d-3)then
 942              bessel=qr*qr/15.0d0-qr**4/210.0d0
 943            else
 944              bessel=bes2_psp5(qr)
 945            end if
 946            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 947 !          work(ir)=rad(ir)*(r(ir)*bes2(arg*rad(ir))*dvwf)
 948            work3(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
 949          end do
 950          call ctrap(mmax,work3,al,result)
 951          ffspl(iq,1,3)=ztor1+result
 952        end do
 953 
 954 !      Compute yp1,ypn for l=2
 955 !      yp1=0 for l=2
 956        yp1=0.0d0
 957 !      ypn=$\int [2 \pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$
 958 !      where x=2 Pi qgrid(mqgrid) r
 959        arg=two_pi*qgrid(mqgrid)
 960        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
 961        qr=arg*rad(1)
 962        if(qr<1.d-3)then
 963          bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0
 964        else
 965          bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr
 966        end if
 967 !      ztor1=( (two_pi*rad(1)**2)*dvwf* (bes1(arg*rad(1))-
 968 !      3.0d0*bes2(arg*rad(1))/(arg*rad(1))) ) * rad(1)/7.0d0
 969        ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/7.0d0
 970        do ir=1,mmax
 971          qr=arg*rad(ir)
 972          if(qr<1.d-3)then
 973            bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0
 974          else
 975            bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr
 976          end if
 977          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
 978 !        work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
 979 !        (bes1(arg*rad(ir))-3.d0*bes2(arg*rad(ir))/(arg*rad(ir))) )
 980          work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel)
 981        end do
 982        call ctrap(mmax,work3,al,result)
 983        ypn=ztor1+result
 984 
 985 !      Fit spline for l=2 Kleinman-Bylander form factor
 986        call spline(qgrid,ffspl(1,1,3),mqgrid,yp1,ypn,ffspl(1,2,3))
 987 
 988      else
 989 !      or else put form factor to 0 for l=1
 990        ffspl(:,:,3)=0.0d0
 991      end if
 992 !    Endif condition of lmax>1
 993    end if
 994 
 995 !  Finished if lmax=2 (highest nonlocal correction)
 996 !  Do l=3 nonlocal form factor if eta(4) not 0 and lmax>=3
 997    if (lmax>2)then
 998      if(ekb(4)/=0.0d0) then
 999 
1000        lp1=4
1001 !      do q=0 separately; f_3(q=0) vanishes
1002        ffspl(1,1,4)=0.0d0
1003 
1004 !      do rest of q points
1005        do iq=2,mqgrid
1006          arg=two_pi*qgrid(iq)
1007          dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1008          qr=arg*rad(1)
1009          if(qr<1.d-3)then
1010            bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0
1011          else
1012            bessel=bes3_psp5(qr)
1013          end if
1014 !        ztor1=(bes3(arg*rad(1))*rad(1)*dvwf)*rad(1)/9.0d0
1015          ztor1=(bessel*rad(1)*dvwf)*rad(1)/9.0d0
1016          do ir=1,mmax
1017            qr=arg*rad(ir)
1018            if(qr<1.d-3)then
1019              bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0
1020            else
1021              bessel=bes3_psp5(qr)
1022            end if
1023            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1024 !          work(ir)=rad(ir)*(rad(ir)*bes3(arg*rad(ir))*dvwf)
1025            work4(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
1026          end do
1027          call ctrap(mmax,work4,al,result)
1028          ffspl(iq,1,4)=ztor1+result
1029        end do
1030 
1031 !      Compute yp1,ypn for l=3
1032 !      yp1=0 for l=3
1033        yp1=0.0d0
1034 !      ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$
1035 !      where x=2 Pi qgrid(mqgrid) r
1036        arg=two_pi*qgrid(mqgrid)
1037        dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1038        qr=arg*rad(1)
1039        if(qr<1.d-3)then
1040          bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0
1041        else
1042          bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr
1043        end if
1044 !      ztor1=( (two_pi*rad(1)**2)*dvwf* (bes2(arg*rad(1))-
1045 !      3.0d0*bes3(arg*rad(1))/(arg*rad(1))) ) * rad(1)/9.0d0
1046        ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/9.0d0
1047        do ir=1,mmax
1048          qr=arg*rad(ir)
1049          if(qr<1.d-3)then
1050            bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0
1051          else
1052            bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr
1053          end if
1054          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1055 !        work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
1056 !        (bes2(arg*rad(ir))-4.d0*bes3(arg*rad(ir))/(arg*rad(ir))) )
1057          work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel)
1058        end do
1059        call ctrap(mmax,work4,al,result)
1060        ypn=ztor1+result
1061 
1062 !      Fit spline for l=3 Kleinman-Bylander form factor
1063        call spline(qgrid,ffspl(1,1,4),mqgrid,yp1,ypn,ffspl(1,2,4))
1064 
1065      else
1066 !      or else put form factor to 0 for l=3
1067        ffspl(:,:,4)=0.0d0
1068      end if
1069 !    Endif condition of lmax>2
1070    end if
1071 
1072 !  Endif condition lmax/=-1
1073  end if
1074 
1075 !DEBUG
1076 !write(std_out,*) 'EKB=',(ekb(iq),iq=1,3)
1077 !write(std_out,*) 'COSKB=',(ckb(iq),iq=1,3)
1078 !ENDDEBUG
1079 
1080  ABI_FREE(work1)
1081  ABI_FREE(work2)
1082  ABI_FREE(work3)
1083  ABI_FREE(work4)
1084 
1085 contains
1086 
1087 function  bes0_psp5(arg)
1088  real(dp) :: bes0_psp5,arg
1089  bes0_psp5=sin(arg)/arg
1090 end function bes0_psp5
1091 
1092 function bes1_psp5(arg)
1093   real(dp) :: bes1_psp5,arg
1094   bes1_psp5=(sin(arg)-arg*cos(arg))/arg**2
1095 end function bes1_psp5
1096 
1097 function bes2_psp5(arg)
1098   real(dp) :: bes2_psp5,arg
1099   bes2_psp5=( (3.0d0-arg**2)*sin(arg)- 3.0d0*arg*cos(arg))/arg**3
1100 end function bes2_psp5
1101 
1102 function bes3_psp5(arg)
1103   real(dp) :: bes3_psp5, arg
1104   bes3_psp5=(15.d0*sin(arg)-15.d0*arg*cos(arg) &
1105               -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4
1106 end function bes3_psp5
1107 
1108 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).

SOURCE

1143 subroutine psp8lo(amesh, epsatm, mmax, mqgrid, qgrid, q2vq, rad, vloc, yp1, ypn, zion)
1144 
1145 !Arguments----------------------------------------------------------
1146 !scalars
1147  integer,intent(in) :: mmax,mqgrid
1148  real(dp),intent(in) :: amesh,zion
1149  real(dp),intent(out) :: epsatm,yp1,ypn
1150 !arrays
1151  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
1152  real(dp),intent(out) :: q2vq(mqgrid)
1153 
1154 !Local variables-------------------------------
1155 !Following parameter controls accuracy of Fourier transform based on qmax
1156 !and represents the minimun number of integration points in one period scalars
1157  integer,parameter :: NPT_IN_2PI=200
1158  integer :: ider,iq,ir,irmu,irn,mesh_mult,mmax_new
1159  real(dp) :: amesh_new,arg,fp1,fpn,qmesh,result,ztor1
1160 !arrays
1161  real(dp),allocatable :: rad_new(:),rvlpz(:),rvlpz_new(:),sprvlpz(:,:),work(:)
1162 
1163 ! *************************************************************************
1164 
1165  ABI_MALLOC(work,(mmax))
1166  ABI_MALLOC(rvlpz,(mmax))
1167 
1168 !Do q=0 separately (compute epsatm)
1169  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
1170 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
1171  do ir=1,mmax
1172    rvlpz(ir)=rad(ir)*vloc(ir)+zion
1173    work(ir)=rad(ir)*rvlpz(ir)
1174  end do
1175 
1176 !Do integral from zero to r(max)
1177  call ctrap(mmax,work,amesh,result)
1178 
1179  epsatm=4.d0*pi*result
1180  q2vq(1)=-zion/pi
1181 
1182 !Find r mesh spacing necessary for accurate integration at qmax
1183  amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1184 
1185 !Choose submultiple of input mesh
1186  mesh_mult=int(amesh/amesh_new) + 1
1187  !mesh_mult = 1  ! DEBUG
1188  mmax_new=mesh_mult*(mmax-1)+1
1189  amesh_new=amesh/dble(mesh_mult)
1190 
1191  ABI_MALLOC(rad_new,(mmax_new))
1192  ABI_MALLOC(rvlpz_new,(mmax_new))
1193 
1194  !print *, "in psp8lo with mesh_mult:", mesh_mult
1195  !print *, "in psp8lo with mmax_new:", mmax_new
1196  !print *, "in psp8lo with amesh_new:", amesh_new
1197 
1198  if(mesh_mult==1) then
1199    rad_new(:)=rad(:)
1200    rvlpz_new(:)=rvlpz(:)
1201  else
1202 !  Set up spline and interpolate to finer mesh.
1203 !  First, compute derivatives at end points
1204    fp1=(-50.d0*rvlpz(1)+96.d0*rvlpz(2)-72.d0*rvlpz(3)+32.d0*rvlpz(4)&
1205 &   -6.d0*rvlpz(5))/(24.d0*amesh)
1206    fpn=(6.d0*rvlpz(mmax-4)-32.d0*rvlpz(mmax-3)+72.d0*rvlpz(mmax-2)&
1207 &   -96.d0*rvlpz(mmax-1)+50.d0*rvlpz(mmax))/(24.d0*amesh)
1208    ABI_MALLOC(sprvlpz,(mmax,2))
1209    work(:)=zero
1210 
1211 !  Spline fit
1212    call spline(rad, rvlpz,mmax,fp1,fpn,sprvlpz(:,2))
1213    sprvlpz(:,1)=rvlpz(:)
1214 
1215 !  Set up new radial mesh
1216    irn=1
1217    do ir=1,mmax-1
1218      do irmu=0,mesh_mult-1
1219        rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
1220        irn=irn+1
1221      end do
1222    end do
1223    rad_new(mmax_new)=rad(mmax)
1224 
1225    ider=0
1226    call splfit(rad,work,sprvlpz,ider,rad_new,rvlpz_new,mmax,mmax_new)
1227 
1228    ABI_FREE(sprvlpz)
1229    ABI_FREE(work)
1230    ABI_MALLOC(work,(mmax_new))
1231  end if
1232 
1233 !Loop over q values
1234  do iq=2,mqgrid
1235    arg=2.d0*pi*qgrid(iq)
1236 
1237 !  Set up integrand
1238    do  ir=1,mmax_new
1239      work(ir)=sin(arg*rad_new(ir))*rvlpz_new(ir)
1240    end do
1241 
1242 !  Do integral from zero to rad(mmax)
1243    call ctrap(mmax_new,work,amesh_new,result)
1244 
1245 !  Store q^2 v(q)
1246    q2vq(iq)=q2vq(1)+2.d0*qgrid(iq)*result
1247 
1248  end do
1249 
1250 !Compute derivatives of q^2 v(q) at ends of interval
1251  qmesh=qgrid(2)-qgrid(1)
1252  yp1=(-50.d0*q2vq(1)+96.d0*q2vq(2)-72.d0*q2vq(3)+32.d0*q2vq(4)&
1253 & -6.d0*q2vq(5))/(24.d0*qmesh)
1254  ypn=(6.d0*q2vq(mqgrid-4)-32.d0*q2vq(mqgrid-3)+72.d0*q2vq(mqgrid-2)&
1255 & -96.d0*q2vq(mqgrid-1)+50.d0*q2vq(mqgrid))/(24.d0*qmesh)
1256 
1257  ABI_FREE(work)
1258  ABI_FREE(rad_new)
1259  ABI_FREE(rvlpz_new)
1260  ABI_FREE(rvlpz)
1261 
1262 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,lmnmax)= 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
  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.

SOURCE

1310 subroutine psp8nl(amesh, ffspl, indlmn, lmax, lmnmax, lnmax, mmax, mqgrid, qgrid, rad, vpspll)
1311 
1312 !Arguments----------------------------------------------------------
1313 !scalars
1314  integer,intent(in) :: lmax,lmnmax,lnmax,mmax,mqgrid
1315  real(dp),intent(in) :: amesh
1316 !arrays
1317  integer,intent(in) :: indlmn(6,lmnmax)
1318  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vpspll(mmax,lnmax)
1319  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
1320 
1321 !Local variables-------------------------------
1322 !Following parameter controls accuracy of Fourier transform based on qmax
1323 !and represents the minimun number of integration points in one period.
1324 !scalars
1325  integer,parameter :: NPT_IN_2PI=200
1326  integer :: iln,iln0,ilmn,iq,ir,irmu,irn,ll,mesh_mult,mmax_new,mvpspll
1327  real(dp) :: amesh_new,arg,c1,c2,c3,c4,dri,qmesh,result,tv,xp,xpm1,xpm2,xpp1,yp1,ypn
1328 !arrays
1329  real(dp) :: sb_out(4)
1330  real(dp),allocatable :: rad_new(:),vpspll_new(:,:),work(:,:),work2(:)
1331 
1332 ! *************************************************************************
1333 
1334  ! Find r mesh spacing necessary for accurate integration at qmax
1335  amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1336 
1337  ! Choose submultiple of input mesh
1338  mesh_mult=int(amesh/amesh_new) + 1
1339  mmax_new=mesh_mult*(mmax-1)+1
1340  amesh_new=amesh/dble(mesh_mult)
1341 
1342  ABI_MALLOC(rad_new,(mmax_new))
1343  ABI_MALLOC(vpspll_new,(mmax_new,lnmax))
1344 
1345  if (mesh_mult == 1) then
1346    rad_new(:)=rad(:)
1347  else
1348    ! Set up new radial mesh
1349    irn=1
1350    do ir=1,mmax-1
1351      do irmu=0,mesh_mult-1
1352        rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
1353        irn=irn+1
1354      end do
1355    end do
1356    rad_new(mmax_new)=rad(mmax)
1357  end if
1358 
1359  ! Interpolate projectors onto new grid if called for
1360  ! Cubic polynomial interpolation is used which is consistent
1361  ! with the original interpolation of these functions from
1362  ! a log grid to the input linear grid.
1363  dri = one/amesh
1364  do irn=1,mmax_new
1365    ! index to find bracketing input mesh points
1366    if(mesh_mult>1) then
1367      ir = irn/mesh_mult + 1
1368      ir = max(ir,2)
1369      ir = min(ir,mmax-2)
1370      ! interpolation coefficients
1371      xp = dri * (rad_new(irn) - rad(ir))
1372      xpp1 = xp + one
1373      xpm1 = xp - one
1374      xpm2 = xp - two
1375      c1 = -xp * xpm1 * xpm2 * sixth
1376      c2 = xpp1 * xpm1 * xpm2 * half
1377      c3 = - xp * xpp1 * xpm2 * half
1378      c4 = xp * xpp1 * xpm1 * sixth
1379 
1380      ! Now do the interpolation on all projectors for this grid point
1381      iln0=0
1382      do ilmn=1,lmnmax
1383        iln=indlmn(5,ilmn)
1384        if (iln>iln0) then
1385          iln0=iln
1386          tv =  c1 * vpspll(ir - 1, iln) &
1387 &         + c2 * vpspll(ir    , iln) &
1388 &         + c3 * vpspll(ir + 1, iln) &
1389 &         + c4 * vpspll(ir + 2, iln)
1390          if(abs(tv)>tol10) then
1391            vpspll_new(irn,iln)=tv
1392            mvpspll=irn
1393          else
1394            vpspll_new(irn,iln)=zero
1395          end if
1396        end if
1397      end do
1398 
1399    else
1400      ! With no mesh multiplication, just copy projectors
1401      ir=irn
1402      iln0=0
1403      do ilmn=1,lmnmax
1404        iln=indlmn(5,ilmn)
1405        if (iln>iln0) then
1406          iln0=iln
1407          tv = vpspll(ir,iln)
1408          if(abs(tv)>tol10) then
1409            vpspll_new(irn,iln)=tv
1410            mvpspll=irn
1411          else
1412            vpspll_new(irn,iln)=zero
1413          end if
1414        end if
1415      end do
1416 
1417    end if
1418  end do ! irn
1419 
1420  ABI_MALLOC(work, (mvpspll,lnmax))
1421 
1422  ! Loop over q values
1423  do iq=1,mqgrid
1424    arg=2.d0*pi*qgrid(iq)
1425 
1426    ! Set up integrands
1427    do ir=1,mvpspll
1428      call sbf8(lmax+1,arg*rad_new(ir),sb_out)
1429      iln0=0
1430      do ilmn=1,lmnmax
1431        iln=indlmn(5,ilmn)
1432        if (iln>iln0) then
1433          iln0=iln
1434          ll=indlmn(1,ilmn)
1435          work(ir,iln)=sb_out(ll+1)*vpspll_new(ir,iln)*rad_new(ir)
1436        end if
1437      end do
1438    end do !ir
1439 
1440    ! Do integral from zero to rad_new(mvpspll)
1441    iln0=0
1442    do ilmn=1,lmnmax
1443      iln=indlmn(5,ilmn)
1444      if (iln>iln0) then
1445        iln0=iln
1446        call ctrap(mvpspll,work(1,iln),amesh_new,result)
1447        ffspl(iq,1,iln)=result
1448      end if
1449    end do
1450  end do ! iq mesh
1451 
1452  ! Fit splines for form factors
1453  ABI_MALLOC(work2,(mqgrid))
1454  qmesh=qgrid(2)-qgrid(1)
1455 
1456  iln0=0
1457  do ilmn=1,lmnmax
1458    iln=indlmn(5,ilmn)
1459    if (iln>iln0) then
1460      iln0=iln
1461      ! Compute derivatives of form factors at ends of interval
1462      yp1=(-50.d0*ffspl(1,1,iln)+96.d0*ffspl(2,1,iln)-72.d0*ffspl(3,1,iln)&
1463 &     +32.d0*ffspl(4,1,iln)- 6.d0*ffspl(5,1,iln))/(24.d0*qmesh)
1464      ypn=(6.d0*ffspl(mqgrid-4,1,iln)-32.d0*ffspl(mqgrid-3,1,iln)&
1465 &     +72.d0*ffspl(mqgrid-2,1,iln)-96.d0*ffspl(mqgrid-1,1,iln)&
1466 &     +50.d0*ffspl(mqgrid,1,iln))/(24.d0*qmesh)
1467 
1468      call spline(qgrid,ffspl(1,1,iln),mqgrid,yp1,ypn,ffspl(1,2,iln))
1469    end if
1470  end do
1471 
1472  ABI_FREE(rad_new)
1473  ABI_FREE(vpspll_new)
1474  ABI_FREE(work)
1475  ABI_FREE(work2)
1476 
1477 end subroutine psp8nl