TABLE OF CONTENTS
ABINIT/psp6cc [ Functions ]
NAME
psp6cc
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for the format 6 of pseudopotentials.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (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 . 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 znucl=nuclear number of atom as specified in psp file
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives Optional output: vh_tnzc(mmax) = Hartree potential induced by density tild_[n_Z+n_core] (pseudized [n_Z+n_core], where n_Z=ions, n_core=core electrons) using a simple pseudization scheme
PARENTS
psp6in
CHILDREN
psden,smooth,spline,splint,vhtnzc
SOURCE
39 #if defined HAVE_CONFIG_H 40 #include "config.h" 41 #endif 42 43 #include "abi_common.h" 44 45 subroutine psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,& 46 & vh_tnzc) ! optional argument 47 48 use defs_basis 49 use m_splines 50 use m_profiling_abi 51 52 use m_numeric_tools, only : smooth 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 'psp6cc' 58 use interfaces_64_psp, except_this_one => psp6cc 59 !End of the abilint section 60 61 implicit none 62 63 !Arguments ------------------------------------ 64 !scalars 65 integer,intent(in) :: mmax,n1xccc 66 real(dp),intent(in) :: rchrg,znucl 67 !arrays 68 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 69 real(dp),intent(out),optional :: vh_tnzc(mmax) 70 71 !Local variables------------------------------- 72 !scalars 73 integer :: i1xccc,irad 74 real(dp) :: der1,dern 75 character(len=500) :: errmsg 76 !arrays 77 real(dp),allocatable :: ff(:),ff1(:),ff2(:),ff3(:),gg(:),gg1(:),gg2(:),gg3(:) 78 real(dp),allocatable :: gg4(:),nc(:),rad(:),work(:),xx(:) 79 80 !********************************************************************** 81 82 ABI_ALLOCATE(ff,(mmax)) 83 ABI_ALLOCATE(ff1,(mmax)) 84 ABI_ALLOCATE(ff2,(mmax)) 85 ABI_ALLOCATE(ff3,(mmax)) 86 ABI_ALLOCATE(rad,(mmax)) 87 ABI_ALLOCATE(gg,(n1xccc)) 88 ABI_ALLOCATE(gg1,(n1xccc)) 89 ABI_ALLOCATE(gg2,(n1xccc)) 90 ABI_ALLOCATE(gg3,(n1xccc)) 91 ABI_ALLOCATE(gg4,(n1xccc)) 92 ABI_ALLOCATE(work,(n1xccc)) 93 ABI_ALLOCATE(xx,(n1xccc)) 94 95 ! 96 !read from pp file the model core charge (ff) and first (ff1) and 97 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid 98 !the input functions contain the 4pi factor, it must be rescaled. 99 100 do irad=1,mmax 101 read(tmp_unit,*, err=10, iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad) 102 ff(irad)=ff(irad)/four_pi 103 ff1(irad)=ff1(irad)/four_pi 104 ff2(irad)=ff2(irad)/four_pi 105 end do 106 107 !Optional output: VHartree(tild_[n_Z+n_core]) 108 if (present(vh_tnzc)) then 109 ABI_ALLOCATE(nc,(mmax)) 110 nc=ff ! n_core 111 call psden(1,ff,mmax,nc,rchrg,rad,ff1=ff1,ff2=ff2) 112 call vhtnzc(ff,rchrg,vh_tnzc,mmax,rad,znucl) 113 ABI_DEALLOCATE(nc) 114 end if 115 116 rad(1)=zero 117 118 !calculate third derivative ff3 on logarithmic grid 119 der1=ff2(1) 120 dern=ff2(mmax) 121 call spline(rad,ff1,mmax,der1,dern,ff3) 122 123 !generate uniform mesh xx in the box cut by rchrg: 124 125 do i1xccc=1,n1xccc 126 xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1) 127 end do 128 129 ! 130 !now interpolate core charge and derivatives on the uniform grid 131 ! 132 !core charge, input=ff, output=gg 133 call splint(mmax,rad,ff,ff2,n1xccc,xx,gg) 134 135 !first derivative input=ff1, output=gg1 136 call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1) 137 138 !normalize gg1 139 gg1(:)=gg1(:)*rchrg 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)=zero 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)=zero 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)=zero 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 xccc1d(:,1)=gg(:) 169 xccc1d(:,2)=gg1(:) 170 xccc1d(:,3)=gg2(:) 171 xccc1d(:,4)=gg3(:) 172 xccc1d(:,5)=gg4(:) 173 174 !WARNING : fifth derivative not yet computed 175 xccc1d(:,6)=zero 176 177 !DEBUG 178 !note: the normalization condition is the following: 179 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg 180 ! 181 !norm=zero 182 !do i1xccc=1,n1xccc 183 !norm = norm + four_pi*rchrg/dble(n1xccc-1)*& 184 !& xx(i1xccc)**2*xccc1d(i1xccc,1) 185 !end do 186 !write(std_out,*) ' norm=',norm 187 ! 188 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 189 !write(std_out,*)' xx gg gg1 ' 190 !do i1xccc=1,n1xccc 191 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 192 !end do 193 !write(std_out,*)' xx gg2 gg3 ' 194 !do i1xccc=1,n1xccc 195 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 196 !end do 197 !write(std_out,*)' xx gg4 gg5 ' 198 !do i1xccc=1,n1xccc 199 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 200 !end do 201 !write(std_out,*)' psp1cc : debug done, stop ' 202 !stop 203 !ENDDEBUG 204 205 ABI_DEALLOCATE(ff) 206 ABI_DEALLOCATE(ff1) 207 ABI_DEALLOCATE(ff2) 208 ABI_DEALLOCATE(ff3) 209 ABI_DEALLOCATE(gg) 210 ABI_DEALLOCATE(gg1) 211 ABI_DEALLOCATE(gg2) 212 ABI_DEALLOCATE(gg3) 213 ABI_DEALLOCATE(gg4) 214 ABI_DEALLOCATE(rad) 215 ABI_DEALLOCATE(work) 216 ABI_DEALLOCATE(xx) 217 218 return 219 220 ! Handle IO error 221 10 continue 222 MSG_ERROR(errmsg) 223 224 end subroutine psp6cc