TABLE OF CONTENTS
ABINIT/psp4cc [ Functions ]
NAME
psp4cc
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for the format 4 of pseudopotentials. This is a even polynomial of 24th order for core density, that is cut off exactly beyond rchrg. It has been produced on 7 May 1992 by M. Teter.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, DCA) 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
The argument of xccc1d is assumed to be normalized, and to vary from xx=0 to 1 (from r=0 to r=xcccrc)
WARNINGS
the fifth derivative is not yet delivered.
PARENTS
psp1in
CHILDREN
spline
SOURCE
43 #if defined HAVE_CONFIG_H 44 #include "config.h" 45 #endif 46 47 #include "abi_common.h" 48 49 subroutine psp4cc(fchrg,n1xccc,xccc1d) 50 51 use defs_basis 52 use m_profiling_abi 53 use m_splines 54 use m_errors 55 56 !This section has been created automatically by the script Abilint (TD). 57 !Do not modify the following lines by hand. 58 #undef ABI_FUNC 59 #define ABI_FUNC 'psp4cc' 60 !End of the abilint section 61 62 implicit none 63 64 !Arguments ------------------------------------ 65 !scalars 66 integer,intent(in) :: n1xccc 67 real(dp),intent(in) :: fchrg 68 !arrays 69 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 70 71 !Local variables------------------------------- 72 !scalars 73 integer :: i1xccc,ider 74 real(dp),parameter :: a10=-0.1156854803757563d5,a12=+0.2371534625455588d5 75 real(dp),parameter :: a14=-0.3138755797827918d5,a16=+0.2582842713241039d5 76 real(dp),parameter :: a18=-0.1200356429115204d5,a20=+0.2405099057118771d4 77 real(dp),parameter :: a2=-0.8480751097855989d1,a4=+0.9684600878284791d2 78 real(dp),parameter :: a6=-0.7490894651588015d3,a8=+0.3670890998130434d4 79 real(dp) :: der1,dern,factor 80 character(len=500) :: message 81 !arrays 82 real(dp),allocatable :: ff(:),ff2(:),work(:),xx(:) 83 real(dp) :: x 84 85 ! ************************************************************************* 86 87 ABI_ALLOCATE(ff,(n1xccc)) 88 ABI_ALLOCATE(ff2,(n1xccc)) 89 ABI_ALLOCATE(work,(n1xccc)) 90 ABI_ALLOCATE(xx,(n1xccc)) 91 92 93 if(n1xccc > 1)then 94 factor=1.0d0/dble(n1xccc-1) 95 do i1xccc=1,n1xccc 96 xx(i1xccc)=(i1xccc-1)*factor 97 end do 98 else 99 write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc 100 MSG_BUG(message) 101 end if 102 103 !Initialization, to avoid some problem with some compilers 104 xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero 105 106 !Take care of each derivative separately 107 do ider=0,2 108 109 if(ider==0)then 110 ! Generate spline fitting for the function gg 111 do i1xccc=1,n1xccc 112 ! ff(i1xccc)=fchrg*gg(xx(i1xccc)) 113 ff(i1xccc)=fchrg*gg_psp4(xx(i1xccc)) 114 end do 115 ! Complete with derivatives at end points 116 der1=0.0d0 117 ! dern=fchrg*gp(1.0d0) 118 dern=fchrg*gp_psp4(1.0d0) 119 else if(ider==1)then 120 ! Generate spline fitting for the function gp 121 do i1xccc=1,n1xccc 122 ! ff(i1xccc)=fchrg*gp(xx(i1xccc)) 123 ff(i1xccc)=fchrg*gp_psp4(xx(i1xccc)) 124 end do 125 ! Complete with derivatives at end points, already estimated 126 der1=xccc1d(1,ider+2) 127 dern=xccc1d(n1xccc,ider+2) 128 else if(ider==2)then 129 ! Generate spline fitting for the function gpp 130 ! (note : the function gpp has already been estimated, for the spline 131 ! fitting of the function gg, but it is replaced here by the more 132 ! accurate analytic derivative) 133 do i1xccc=1,n1xccc 134 x=xx(i1xccc) 135 ff(i1xccc)=fchrg*(gpp_1_psp4(x)+gpp_2_psp4(x)+gpp_3_psp4(x)) 136 ! ff(i1xccc)=fchrg*gpp(xx(i1xccc)) 137 end do 138 ! Complete with derivatives of end points 139 der1=xccc1d(1,ider+2) 140 dern=xccc1d(n1xccc,ider+2) 141 end if 142 143 ! Produce second derivative numerically, for use with splines 144 call spline(xx,ff,n1xccc,der1,dern,ff2) 145 xccc1d(:,ider+1)=ff(:) 146 xccc1d(:,ider+3)=ff2(:) 147 148 end do 149 150 xccc1d(:,6)=zero 151 152 ABI_DEALLOCATE(ff) 153 ABI_DEALLOCATE(ff2) 154 ABI_DEALLOCATE(work) 155 ABI_DEALLOCATE(xx) 156 157 !DEBUG 158 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 159 !write(std_out,*)' xx gg gp ' 160 !do i1xccc=1,n1xccc 161 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 162 !end do 163 !write(std_out,*)' xx gpp gg2 ' 164 !do i1xccc=1,n1xccc 165 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 166 !end do 167 !write(std_out,*)' xx gp2 gpp2 ' 168 !do i1xccc=1,n1xccc 169 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 170 !end do 171 !write(std_out,*)' psp1cc : debug done, stop ' 172 !stop 173 !ENDDEBUG 174 175 contains 176 177 function gg_psp4(x) 178 !Expression of 7 May 1992 179 180 !This section has been created automatically by the script Abilint (TD). 181 !Do not modify the following lines by hand. 182 #undef ABI_FUNC 183 #define ABI_FUNC 'gg_psp4' 184 !End of the abilint section 185 186 real(dp) :: gg_psp4 187 real(dp),intent(in) :: x 188 gg_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 189 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 190 & x**2*(a18+x**2*(a20))))))))))) *(1.0d0-x**2)**2 191 end function gg_psp4 192 193 function gp_psp4(x) 194 !gp(x) is the derivative of gg(x) wrt x 195 196 !This section has been created automatically by the script Abilint (TD). 197 !Do not modify the following lines by hand. 198 #undef ABI_FUNC 199 #define ABI_FUNC 'gp_psp4' 200 !End of the abilint section 201 202 real(dp) :: gp_psp4 203 real(dp),intent(in) :: x 204 gp_psp4=2.d0*x*((a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( & 205 & 4.d0*a8+x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( & 206 & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(10.d0*a20))))))))))*& 207 & (1.d0-x**2)**2 & 208 & -2.0d0*(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 209 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 210 & x**2*(a18+x**2*a20)))))))))) *(1.0d0-x**2) ) 211 end function gp_psp4 212 213 function gpp_1_psp4(x) 214 !gpp(x) is the second derivative of gg(x) wrt x 215 216 !This section has been created automatically by the script Abilint (TD). 217 !Do not modify the following lines by hand. 218 #undef ABI_FUNC 219 #define ABI_FUNC 'gpp_1_psp4' 220 !End of the abilint section 221 222 real(dp) :: gpp_1_psp4 223 real(dp),intent(in) :: x 224 gpp_1_psp4= ( 2.d0*a4+ x**2*(3.d0*2.d0*a6 +x**2*( & 225 & 4.d0*3.d0*a8+ x**2*(5.d0*4.d0*a10+x**2*( & 226 & 6.d0*5.d0*a12+x**2*(7.d0*6.d0*a14+x**2*( & 227 & 8.d0*7.d0*a16+x**2*(9.d0*8.d0*a18+x**2*( & 228 & 10.d0*9.d0*a20) & 229 & ))))))))*(2.d0*x*(1.d0-x**2))**2 230 end function gpp_1_psp4 231 232 function gpp_2_psp4(x) 233 234 235 !This section has been created automatically by the script Abilint (TD). 236 !Do not modify the following lines by hand. 237 #undef ABI_FUNC 238 #define ABI_FUNC 'gpp_2_psp4' 239 !End of the abilint section 240 241 real(dp) :: gpp_2_psp4 242 real(dp),intent(in) :: x 243 gpp_2_psp4=(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( & 244 & 4.d0*a8 +x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( & 245 & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*( & 246 & 10.d0*a20) & 247 & )))))))))*(1.d0-x**2)*2*(1.d0-9.d0*x**2) 248 end function gpp_2_psp4 249 250 function gpp_3_psp4(x) 251 252 253 !This section has been created automatically by the script Abilint (TD). 254 !Do not modify the following lines by hand. 255 #undef ABI_FUNC 256 #define ABI_FUNC 'gpp_3_psp4' 257 !End of the abilint section 258 259 real(dp) :: gpp_3_psp4 260 real(dp),intent(in) :: x 261 gpp_3_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 262 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 263 & x**2*(a18+x**2*a20 & 264 & ))))))))))*(1.0d0-3.d0*x**2)*(-4.d0) 265 end function gpp_3_psp4 266 267 end subroutine psp4cc