TABLE OF CONTENTS
- ABINIT/cc_derivatives
- ABINIT/gg1cc
- ABINIT/gp1cc
- ABINIT/gpp1cc
- ABINIT/m_psptk
- ABINIT/psp1cc
- ABINIT/psp5lo
- ABINIT/psp5nl
- ABINIT/psp8lo
- ABINIT/psp8nl
ABINIT/cc_derivatives [ 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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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