TABLE OF CONTENTS
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)
PARENTS
psp1cc
CHILDREN
SOURCE
206 subroutine gg1cc(gg1cc_xx,xx) 207 208 use defs_basis 209 use m_profiling_abi 210 211 !This section has been created automatically by the script Abilint (TD). 212 !Do not modify the following lines by hand. 213 #undef ABI_FUNC 214 #define ABI_FUNC 'gg1cc' 215 !End of the abilint section 216 217 implicit none 218 219 !Arguments ------------------------------------ 220 !scalars 221 real(dp),intent(in) :: xx 222 real(dp),intent(out) :: gg1cc_xx 223 224 !Local variables ------------------------------------------- 225 !The c s are coefficients for Taylor expansion of the analytic 226 !form near xx=0, 1/2, and 1. 227 !scalars 228 real(dp) :: c21=4.d0/9.d0,c22=-40.d0/27.d0,c23=20.d0/3.d0-16.d0*pi**2/27.d0 229 real(dp) :: c24=-4160.d0/243.d0+160.d0*pi**2/81.d0,c31=1.d0/36.d0 230 real(dp) :: c32=-25.d0/108.d0,c33=485.d0/432.d0-pi**2/27.d0 231 real(dp) :: c34=-4055.d0/972.d0+25.d0*pi**2/81.d0 232 233 ! ************************************************************************* 234 235 !Cut off beyond 3/gcut=xcccrc 236 if (xx>3.0d0) then 237 gg1cc_xx=0.0d0 238 ! Take care of difficult limits near x=0, 1/2, and 1 239 else if (abs(xx)<=1.d-09) then 240 gg1cc_xx=1.d0 241 else if (abs(xx-0.5d0)<=1.d-04) then 242 ! (this limit and next are more troublesome for numerical cancellation) 243 gg1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*c24)) 244 else if (abs(xx-1.d0)<=1.d-04) then 245 gg1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*c34)) 246 else 247 ! The following is the square of the Fourier transform of a 248 ! function built out of two spherical bessel functions in G 249 ! space and cut off absolutely beyond gcut 250 gg1cc_xx=(sin(2.0d0*pi*xx)/( (2.0d0*pi*xx) * & 251 & (1.d0-4.0d0*xx**2)*(1.d0-xx**2) ) )**2 252 end if 253 254 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)}$
PARENTS
psp1cc
CHILDREN
SOURCE
283 subroutine gp1cc(gp1cc_xx,xx) 284 285 use defs_basis 286 use m_profiling_abi 287 288 !This section has been created automatically by the script Abilint (TD). 289 !Do not modify the following lines by hand. 290 #undef ABI_FUNC 291 #define ABI_FUNC 'gp1cc' 292 !End of the abilint section 293 294 implicit none 295 296 !Arguments ------------------------------------ 297 !scalars 298 real(dp),intent(in) :: xx 299 real(dp),intent(out) :: gp1cc_xx 300 301 !Local variables ------------------------------------------- 302 !scalars 303 real(dp),parameter :: c11=20.d0-8.d0*pi**2/3.d0 304 real(dp),parameter :: c12=268.d0-160.d0/3.d0*pi**2+128.d0/45.d0*pi**4 305 real(dp),parameter :: c21=-40.d0/27.d0,c22=40.d0/3.d0-32.d0*pi**2/27.d0 306 real(dp),parameter :: c23=-4160.d0/81.d0+160.d0*pi**2/27.d0 307 real(dp),parameter :: c24=157712.d0/729.d0-320.d0*pi**2/9.d0+512.d0*pi**4/405.d0 308 real(dp),parameter :: c25=-452200.d0/729.d0+83200.d0*pi**2/729.d0-1280.d0*pi**4/243.d0 309 real(dp),parameter :: c31=-25.d0/108.d0,c32=485.d0/216.d0-2.d0*pi**2/27.d0 310 real(dp),parameter :: c33=-4055.d0/324.d0+25.d0*pi**2/27.d0 311 real(dp),parameter :: c34=616697.d0/11664.d0-485.d0*pi**2/81.d0+32.d0*pi**4/405.d0 312 real(dp),parameter :: c35=-2933875.d0/15552.d0+20275.d0*pi**2/729.d0-200.d0*pi**4/243.d0 313 real(dp),parameter :: two_pim1=1.0d0/two_pi 314 real(dp) :: denom,phi,phip 315 316 ! ************************************************************************* 317 318 !Cut off beyond r=3*xcccrc is already done at the calling level 319 if (xx>1.001d0) then 320 ! The part that follows will be repeated later, but written in this way, 321 ! only one "if" condition is tested in most of the cases (1.001 < x < 3.0) 322 denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2)) 323 phi=denom*sin(two_pi*xx)*two_pim1 324 phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi) 325 gp1cc_xx=2.d0*phi*phip 326 ! Handle limits where denominator vanishes 327 else if (abs(xx)<1.d-03) then 328 gp1cc_xx=xx*(c11+xx**2*c12) 329 else if (abs(xx-0.5d0)<=1.d-03) then 330 gp1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*(c24+(xx-0.5d0)*c25))) 331 else if (abs(xx-1.d0)<=1.d-03) then 332 gp1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*(c34+(xx-1.0d0)*c35))) 333 else 334 ! Here is the repeated part ... 335 denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2)) 336 phi=denom*sin(two_pi*xx)*two_pim1 337 phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi) 338 gp1cc_xx=2.d0*phi*phip 339 end if 340 341 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.
PARENTS
psp1cc
CHILDREN
SOURCE
365 subroutine gpp1cc(gpp1cc_xx,xx) 366 367 use defs_basis 368 use m_profiling_abi 369 370 !This section has been created automatically by the script Abilint (TD). 371 !Do not modify the following lines by hand. 372 #undef ABI_FUNC 373 #define ABI_FUNC 'gpp1cc' 374 !End of the abilint section 375 376 implicit none 377 378 !Arguments ------------------------------------ 379 !scalars 380 real(dp),intent(in) :: xx 381 real(dp),intent(out) :: gpp1cc_xx 382 383 !Local variables ------------------------------------------- 384 !scalars 385 real(dp),parameter :: c1=20.d0-8.d0*pi**2/3.d00 386 real(dp),parameter :: c2=40.d0/3.d0-32.d0*pi**2/27.d0 387 real(dp),parameter :: c3=-8320.d0/81.d0+320.d0*pi**2/27.d0 388 real(dp),parameter :: c4=157712.d0/243.d0-320.d0*pi**2/3.d0+512.d0*pi**4/135.d0 389 real(dp),parameter :: c5=-18088.d2/729.d0+3328.d2*pi**2/729.d0-5120.d0*pi**4/243.d0 390 real(dp),parameter :: c6=485.d0/216.d0-2.d0*pi**2/27.d0 391 real(dp),parameter :: c7=-4055.d0/162.d0+50.d0*pi**2/27.d0 392 real(dp),parameter :: c8=616697.d0/3888.d0-485.d0*pi**2/27.d0+32.d0*pi**4/135.d0 393 real(dp),parameter :: c9=-2933875.d0/3888.d0+81100.d0*pi**2/729.d0-800.d0*pi**4/243.d0 394 real(dp) :: t1,t10,t100,t11,t12,t120,t121,t122,t127,t138,t14,t140,t15,t152 395 real(dp) :: t157,t16,t160,t17,t174,t175,t18,t19,t2,t20,t21,t23,t24,t3,t31,t33 396 real(dp) :: t34,t4,t41,t42,t44,t45,t46,t5,t54,t55,t56,t57,t6,t62,t64,t65,t7 397 real(dp) :: t72,t78,t79,t8,t85,t9,t93 398 399 ! ************************************************************************* 400 401 if (xx>3.0d0) then 402 ! Cut off beyond 3/gcut=3*xcccrc 403 gpp1cc_xx=0.0d0 404 ! Take care of difficult limits near xx=0, 1/2, and 1 405 else if (abs(xx)<=1.d-09) then 406 gpp1cc_xx=c1 407 else if (abs(xx-0.5d0)<=1.d-04) then 408 ! (this limit and next are more troublesome for numerical cancellation) 409 gpp1cc_xx=c2+(xx-0.5d0)*(c3+(xx-0.5d0)*(c4+(xx-0.5d0)*c5)) 410 else if (abs(xx-1.d0)<=1.d-04) then 411 gpp1cc_xx=c6+(xx-1.0d0)*(c7+(xx-1.0d0)*(c8+(xx-1.0d0)*c9)) 412 else 413 414 ! Should fix up this Maple fortran later 415 t1 = xx**2 416 t2 = 1/t1 417 t3 = 1/Pi 418 t4 = 2*xx 419 t5 = t4-1 420 t6 = t5**2 421 t7 = 1/t6 422 t8 = t4+1 423 t9 = t8**2 424 t10 = 1/t9 425 t11 = xx-1 426 t12 = t11**2 427 t14 = 1/t12/t11 428 t15 = xx+1 429 t16 = t15**2 430 t17 = 1/t16 431 t18 = Pi*xx 432 t19 = sin(t18) 433 t20 = cos(t18) 434 t21 = t20**2 435 t23 = t19*t21*t20 436 t24 = t17*t23 437 t31 = t19**2 438 t33 = t31*t19*t20 439 t34 = t17*t33 440 t41 = Pi**2 441 t42 = 1/t41 442 t44 = 1/t16/t15 443 t45 = t31*t21 444 t46 = t44*t45 445 t54 = 1/t1/xx 446 t55 = 1/t12 447 t56 = t55*t46 448 t57 = t10*t56 449 t62 = t9**2 450 t64 = t17*t45 451 t65 = t55*t64 452 t72 = 1/t9/t8 453 t78 = t14*t64 454 t79 = t10*t78 455 t85 = t12**2 456 t93 = t21**2 457 t100 = t31**2 458 t120 = 1/t6/t5 459 t121 = t55*t34 460 t122 = t10*t121 461 t127 = t16**2 462 t138 = t6**2 463 t140 = t10*t65 464 t152 = t72*t65 465 t157 = t7*t140 466 t160 = t1**2 467 t174 = t55*t24 468 t175 = t10*t174 469 gpp1cc_xx = 8*t2*t3*t7*t10*t14*t34+8*t2*t42*t7*t10*t14*t46& 470 & -8*t2*t3*t7*t10*t14*t24+8*t2*t3*t7*t10*t55*t44*t33+& 471 & 6*t2*t42*t7*t10*t55/t127*t45+24*t2*t42/t138*t140+& 472 & 16*t54*t42*t120*t140+16*t2*t3*t120*t122+16*t2& 473 & *t42*t7*t72*t78-8*t2*t3*t7*t10*t55*t44*t23-8*t54*t3*t7*t175& 474 & +2*t2*t7*t10*t55*t17*t100+2*t2*t7*t10*t55*t17*t93+& 475 & 8*t54*t42*t7*t79+16*t2*t42*t7*t72*t56+6*t2*t42*t7*t10/t85& 476 & *t64+24*t2*t42*t7/t62*t65+8*t54*t42*t7*t57-& 477 & 16*t2*t3*t7*t72*t174+8*t54*t3*t7*t122-16*t2*t3*t120*t175& 478 & +16*t2*t42*t120*t79+16*t2*t42*t120*t57+16*t54*t42*t7*t152+& 479 & 32*t2*t42*t120*t152+16*t2*t3*t7*t72*t121-12*t2*t157+& 480 & 6/t160*t42*t157 481 end if 482 483 end subroutine gpp1cc
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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, DCA, MM) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
fchrg=magnitude of the core charge correction n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
NOTES
This is a revised expression for core density (5 Nov 1992) :
density(r)=fchrg*gg(xx)
with
$ gg(xx)=(\frac{\sin(2\pi xx)}{(2\pi xx)(1-4 xx^2)(1-xx^2)})^2 $
and
$ xx=\frac{r}{rchrg}=\frac{r}{xcccrc/3.0d0}=3*\frac{r}{xcccrc}=3*yy $
Code for gg(xx), gp(xx), and gpp(xx) has been tested by numerical
derivatives--looks ok. gpp(x) should still be rewritten.
The argument of xccc1d is assumed to be normalized, and to vary
from yy=0 to 1 (from r=0 to r=xcccrc, or from xx=0 to 3)
Thus :
xccc1d(yy)=fchrg*[\frac{\sin(2*\pi*(3yy))}
{(6*\pi*(3yy))(1-4*(3yy)^2)(1-(3yy)^2)}]^2
WARNINGS
Warning : the fifth derivative is not yet delivered.
PARENTS
psp1in,psp5in
CHILDREN
gg1cc,gp1cc,gpp1cc,spline
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 subroutine psp1cc(fchrg,n1xccc,xccc1d) 62 63 use defs_basis 64 use m_splines 65 use m_errors 66 use m_profiling_abi 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'psp1cc' 72 use interfaces_64_psp, except_this_one => psp1cc 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: n1xccc 80 real(dp),intent(in) :: fchrg 81 !arrays 82 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 83 84 !Local variables------------------------------- 85 !scalars 86 integer :: i1xccc,ider 87 real(dp) :: der1,dern,factor,gg1cc_xx,gp1cc_xx,gpp1cc_xx,xx 88 character(len=500) :: message 89 !arrays 90 real(dp),allocatable :: ff(:),ff2(:),work(:),yy(:) 91 92 ! ************************************************************************* 93 94 ABI_ALLOCATE(ff,(n1xccc)) 95 ABI_ALLOCATE(ff2,(n1xccc)) 96 ABI_ALLOCATE(work,(n1xccc)) 97 ABI_ALLOCATE(yy,(n1xccc)) 98 99 if(n1xccc > 1)then 100 factor=one/dble(n1xccc-1) 101 do i1xccc=1,n1xccc 102 yy(i1xccc)=(i1xccc-1)*factor 103 end do 104 else 105 write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc 106 MSG_BUG(message) 107 end if 108 109 !Initialization, to avoid some problem with some compilers 110 xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero 111 112 !Take care of each derivative separately 113 do ider=0,2 114 115 if(ider==0)then 116 ! Generate spline fitting for the function gg 117 do i1xccc=1,n1xccc 118 xx=three*yy(i1xccc) 119 call gg1cc(gg1cc_xx,xx) 120 ff(i1xccc)=fchrg*gg1cc_xx 121 end do 122 ! Complete with derivatives at end points 123 der1=zero 124 call gp1cc(gp1cc_xx,three) 125 dern=three*fchrg*gp1cc_xx 126 else if(ider==1)then 127 ! Generate spline fitting for the function gp 128 do i1xccc=1,n1xccc 129 xx=three*yy(i1xccc) 130 call gp1cc(gp1cc_xx,xx) 131 ff(i1xccc)=three*fchrg*gp1cc_xx 132 end do 133 ! Complete with derivatives at end points, already estimated 134 der1=xccc1d(1,ider+2) 135 dern=xccc1d(n1xccc,ider+2) 136 else if(ider==2)then 137 ! Generate spline fitting for the function gpp 138 ! (note : the function gpp has already been estimated, for the spline 139 ! fitting of the function gg, but it is replaced here by the more 140 ! accurate analytic derivative) 141 do i1xccc=1,n1xccc 142 xx=three*yy(i1xccc) 143 call gpp1cc(gpp1cc_xx,xx) 144 ff(i1xccc)=9.0_dp*fchrg*gpp1cc_xx 145 end do 146 ! Complete with derivatives of end points 147 der1=xccc1d(1,ider+2) 148 dern=xccc1d(n1xccc,ider+2) 149 end if 150 151 ! Produce second derivative numerically, for use with splines 152 call spline(yy,ff,n1xccc,der1,dern,ff2) 153 xccc1d(:,ider+1)=ff(:) 154 xccc1d(:,ider+3)=ff2(:) 155 156 end do 157 158 xccc1d(:,6)=zero 159 160 !DEBUG 161 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 162 !write(std_out,*)' yy gg gp ' 163 !do i1xccc=1,n1xccc 164 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 165 !end do 166 !write(std_out,*)' yy gpp gg2 ' 167 !do i1xccc=1,n1xccc 168 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 169 !end do 170 !write(std_out,*)' yy gp2 gpp2 ' 171 !do i1xccc=1,n1xccc 172 !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 173 !end do 174 !write(std_out,*)' psp1cc : debug done, stop ' 175 !stop 176 !ENDDEBUG 177 178 ABI_DEALLOCATE(ff) 179 ABI_DEALLOCATE(ff2) 180 ABI_DEALLOCATE(work) 181 ABI_DEALLOCATE(yy) 182 183 end subroutine psp1cc