TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (AF,DRH) 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
mmax=maximum number of points in real space grid in the psp file n1xccc=dimension of xccc1d ; 0 if no XC core correction is used rchrg=cut-off radius for the core density rad=radial grid points ff=core charge at points in rad ff1=first derivative of ff on log grid ff2=second derivative of ff on log grid
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
PARENTS
psp6cc_drh,upf2abinit
CHILDREN
spline,splint
NOTES
Test version by DRH - requires very smooth model core charge
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 subroutine cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d) 49 50 use defs_basis 51 use m_splines 52 use m_profiling_abi 53 54 !This section has been created automatically by the script Abilint (TD). 55 !Do not modify the following lines by hand. 56 #undef ABI_FUNC 57 #define ABI_FUNC 'cc_derivatives' 58 !End of the abilint section 59 60 implicit none 61 62 !Arguments ------------------------------------ 63 ! scalars 64 integer,intent(in) :: mmax,n1xccc 65 real(dp),intent(in) :: rchrg 66 !arrays 67 real(dp),intent(in) :: rad(mmax),ff(mmax),ff1(mmax),ff2(mmax) 68 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 69 70 !Local variables------------------------------- 71 ! scalars 72 integer :: i1xccc 73 real(dp) :: der1,dern 74 !arrays 75 real(dp),allocatable :: ff3(:),ff4(:),gg(:),gg1(:),gg2(:) 76 real(dp),allocatable :: gg3(:),gg4(:),work(:),xx(:) 77 78 ! ************************************************************************* 79 ABI_ALLOCATE(ff3,(mmax)) 80 ABI_ALLOCATE(ff4,(mmax)) 81 ABI_ALLOCATE(gg,(n1xccc)) 82 ABI_ALLOCATE(gg1,(n1xccc)) 83 ABI_ALLOCATE(gg2,(n1xccc)) 84 ABI_ALLOCATE(gg3,(n1xccc)) 85 ABI_ALLOCATE(gg4,(n1xccc)) 86 ABI_ALLOCATE(work,(mmax)) 87 ABI_ALLOCATE(xx,(n1xccc)) 88 89 write(std_out,*) 'cc_derivatives : enter' 90 91 !calculate third derivative ff3 on logarithmic grid 92 der1=ff2(1) 93 dern=ff2(mmax) 94 call spline(rad,ff1,mmax,der1,dern,ff3) 95 96 !calculate fourth derivative ff4 on logarithmic grid 97 der1=0.d0 98 dern=0.d0 99 call spline(rad,ff2,mmax,der1,dern,ff4) 100 101 !generate uniform mesh xx in the box cut by rchrg: 102 103 do i1xccc=1,n1xccc 104 xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1) 105 end do 106 ! 107 !now interpolate core charge and derivatives on the uniform grid 108 ! 109 !core charge, input=ff, output=gg 110 call splint(mmax,rad,ff,ff2,n1xccc,xx,gg) 111 112 !first derivative input=ff1, output=gg1 113 call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1) 114 115 !normalize gg1 116 !gg1(:)=gg1(:)*rchrg 117 118 !second derivative input=ff2, output=gg2 119 call splint(mmax,rad,ff2,ff4,n1xccc,xx,gg2) 120 121 !normalize gg2 122 !gg2(:)=gg2(:)*rchrg**2 123 124 !reallocate work otherwise the calls to spline crash (n1xccc /= mmax) 125 ABI_DEALLOCATE(work) 126 ABI_ALLOCATE(work,(n1xccc)) 127 128 !recalculate 3rd derivative consistent with spline fit to first derivative 129 !on linear grid 130 der1=gg2(1) 131 dern=gg2(n1xccc) 132 call spline(xx,gg1,n1xccc,der1,dern,gg3) 133 134 !calculate 4th derivative consistent with spline fit to second derivative 135 !on linear grid 136 der1=0.0d0 137 dern=0.0d0 138 call spline(xx,gg2,n1xccc,der1,dern,gg4) 139 140 141 !now calculate second to fourth derivative by forward differences 142 !to avoid numerical noise uses a smoothing function 143 ! 144 !call smooth(gg1,n1xccc,10) 145 146 !gg2(n1xccc)=0.0 147 !do i1xccc=1,n1xccc-1 148 !gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1) 149 !end do 150 151 !call smooth(gg2,n1xccc,10) 152 153 !gg3(n1xccc)=0.0 154 !do i1xccc=1,n1xccc-1 155 !gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1) 156 !end do 157 158 !call smooth(gg3,n1xccc,10) 159 160 !gg4(n1xccc)=0.0 161 !do i1xccc=1,n1xccc-1 162 !gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1) 163 !end do 164 165 !call smooth(gg4,n1xccc,10) 166 167 !write on xcc1d 168 !normalize to unit range usage later in program 169 xccc1d(:,1)=gg(:) 170 xccc1d(:,2)=gg1(:)*rchrg 171 xccc1d(:,3)=gg2(:)*rchrg**2 172 xccc1d(:,4)=gg3(:)*rchrg**3 173 xccc1d(:,5)=gg4(:)*rchrg**4 174 !***drh test 175 !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc 176 !***end drh test 177 178 179 !DEBUG 180 !note: the normalization condition is the following: 181 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg 182 ! 183 !norm=0.d0 184 !do i1xccc=1,n1xccc 185 !norm = norm + 4.d0*pi*rchrg/dble(n1xccc-1)*& 186 !& xx(i1xccc)**2*xccc1d(i1xccc,1) 187 !end do 188 !write(std_out,*) ' norm=',norm 189 ! 190 !write(std_out,*)' psp6cc_drh : output of core charge density and derivatives ' 191 !write(std_out,*)' xx gg gg1 ' 192 !do i1xccc=1,n1xccc 193 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 194 !end do 195 !write(std_out,*)' xx gg2 gg3 ' 196 !do i1xccc=1,n1xccc 197 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 198 !end do 199 !write(std_out,*)' xx gg4 gg5 ' 200 !do i1xccc=1,n1xccc 201 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 202 !end do 203 !write(std_out,*)' psp1cc : debug done, stop ' 204 !stop 205 !ENDDEBUG 206 207 ABI_DEALLOCATE(ff3) 208 ABI_DEALLOCATE(ff4) 209 ABI_DEALLOCATE(gg) 210 ABI_DEALLOCATE(gg1) 211 ABI_DEALLOCATE(gg2) 212 ABI_DEALLOCATE(gg3) 213 ABI_DEALLOCATE(gg4) 214 ABI_DEALLOCATE(work) 215 ABI_DEALLOCATE(xx) 216 217 end subroutine cc_derivatives