./64_psp/psp9cc.F90

ABINIT/psp9cc

ABINIT/psp9cc [ Functions ]

NAME

 psp9cc

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for format 9 of the pseudopotentials (PSML).

 Copyright (C) 2017-2018 ABINIT group (YP)
 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

OUTPUT

  rchrg=cut-off radius for the core density
  xccc1d(n1xccc,6)= 1D core charge function and its four first derivatives

NOTES

  This routine will be built only if PSML support is enabled.

PARENTS

      psp9in

CHILDREN

      dgesv

SOURCE

 37 #if defined HAVE_CONFIG_H
 38 #include "config.h"
 39 #endif
 40 
 41 #include "abi_common.h"
 42 
 43 subroutine psp9cc(psxml,mmax,n1xccc,rad,rchrg,xccc1d)
 44 
 45  use defs_basis
 46  use m_errors
 47  use m_profiling_abi
 48  use m_psml
 49 
 50 !This section has been created automatically by the script Abilint (TD).
 51 !Do not modify the following lines by hand.
 52 #undef ABI_FUNC
 53 #define ABI_FUNC 'psp9cc'
 54 !End of the abilint section
 55 
 56  implicit none
 57 
 58 !Arguments ------------------------------------
 59 !scalars
 60  integer,intent(in) :: mmax,n1xccc
 61  real(dp),intent(out) :: rchrg
 62  type(ps_t),intent(in) :: psxml
 63 !arrays
 64  real(dp),intent(in) :: rad(mmax)
 65  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
 66 
 67 !Local variables-------------------------------
 68 !scalars
 69  integer :: i1xccc,idum,irad,jj
 70  real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx,twelvth
 71  character(len=500) :: message
 72 !arrays
 73  integer :: iwork(8)
 74  real(dp) :: rscale(5),dpoly(6,6),vpoly(6)
 75  real(dp),allocatable :: ff(:,:)
 76 
 77 !**********************************************************************
 78 
 79 !Check that rad grid is linear starting at zero
 80  amesh=rad(2)-rad(1)
 81  damesh=zero
 82  do irad=2,mmax-1
 83    damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
 84  end do
 85 
 86  if(damesh>tol8 .or. rad(1)/=zero) then
 87    write(message, '(5a)' )&
 88 &   'Pseudopotential input file requires linear radial mesh',ch10,&
 89 &   'starting at zero.',ch10,&
 90 &   'Action : check your pseudopotential input file.'
 91    MSG_ERROR(message)
 92  end if
 93 
 94  ABI_ALLOCATE(ff,(mmax,5))
 95 
 96  dri = one / amesh
 97  pi4i = quarter / pi
 98  twelvth = one / 12.0_dp
 99 
100 !Read from pp file the model core charge and calculate its first 4 derivatives
101 !assumed to be on a linear grid starting at zero.
102 !The input functions contain the 4pi factor, and must be rescaled.
103 
104 !Store the value of the pseudo-core charge.
105  ff(:,:) = zero
106  do jj=1,mmax
107    ff(jj,1) = ps_CoreCharge_Value(psxml,rad(jj))
108  end do
109 
110 !Calculate 4 first derivatives with 5-point stencil, except borders
111  do irad=3,mmax-2
112    ff(irad,2) = (-ff(irad+2,1) + 8.0d0*ff(irad+1,1) - &
113 &   8.0d0*ff(irad-1,1) + ff(irad-2,1)) * twelvth * dri
114    ff(irad,3) = (-ff(irad+2,1) + 16.0d0*ff(irad+1,1) - 30.0d0*ff(irad,1) + &
115 &   16.0d0*ff(irad-1,1) - ff(irad-2,1)) * twelvth * dri * dri
116    ff(irad,4) = (ff(irad+2,1) - 2.0d0*ff(irad+1,1) + &
117 &   2.0d0*ff(irad-1,1) - ff(irad-2,1)) * half * dri * dri * dri
118    ff(irad,5) = (ff(irad+2,1) - 4.0d0*ff(irad+1,1) + 6.0d0*ff(irad,1) - &
119 &   4.0d0*ff(irad-1,1) + ff(irad-2,1)) * dri * dri * dri * dri
120  end do
121 
122 !Add border near zero using polynomial fit
123  dpoly(:,:) = zero
124  dpoly(:,1) = one
125  vpoly(:) = zero
126  vpoly(1) = ff(1,1)
127  do irad=2,6
128    do jj=1,6
129      dpoly(irad,jj) = rad(irad)**(jj-1)
130    end do
131    vpoly(irad) = ff(irad,1)
132  end do
133  call dgesv(6,1,dpoly,6,iwork,vpoly,6,idum)
134 
135  do irad=1,2
136    ff(irad,2) = &
137 &   vpoly(2) + 2.0d0*vpoly(3)*rad(irad) + &
138 &   3.0d0*vpoly(4)*rad(irad)*rad(irad) + &
139 &   4.0d0*vpoly(5)*rad(irad)*rad(irad)*rad(irad) + &
140 &   5.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)*rad(irad)
141    ff(irad,3) = &
142 &   2.0d0*vpoly(3)*rad(irad) + &
143 &   6.0d0*vpoly(4)*rad(irad) + &
144 &   12.0d0*vpoly(5)*rad(irad)*rad(irad) + &
145 &   20.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)
146    ff(irad,4) = &
147 &   6.0d0*vpoly(4) + &
148 &   24.0d0*vpoly(5)*rad(irad) + &
149 &   60.0d0*vpoly(6)*rad(irad)*rad(irad)
150    ff(irad,5) = 24.0d0*vpoly(5) + &
151 &   120.0d0*vpoly(6)*rad(irad)
152  end do
153 
154 !Make linear approximation for the tail near mmax
155  do irad=1,2
156    ff(mmax-2+irad,2) = ff(mmax-2,2) + irad * (ff(mmax-2,2) - ff(mmax-3,2))
157    ff(mmax-2+irad,3) = ff(mmax-2,3) + irad * (ff(mmax-2,3) - ff(mmax-3,3))
158    ff(mmax-2+irad,4) = ff(mmax-2,4) + irad * (ff(mmax-2,4) - ff(mmax-3,4))
159    ff(mmax-2+irad,5) = ff(mmax-2,5) + irad * (ff(mmax-2,5) - ff(mmax-3,5))
160  end do
161 
162 !Renormalize core charge
163 ! ff(:,:) = ff(:,:) * pi4i
164 
165 !determine xcccrc where the pseudocore becomes 0
166 !This is a difference with respect the Hamann's treatment of the core
167 !charge when reading PSP8.
168 !In Hamann's case (PSP8), xcccrc = rchrg, and this value is
169 !introduced in the pseudopotential input file.
170 !rchrg is not included in the PSML format
171  rchrg = zero
172  do jj=mmax,1,-1
173    if (ff(jj,1) > tol13) then
174      rchrg=rad(jj)
175      exit
176    end if
177  end do
178 
179 !Check that input rchrg is consistent with last grid point
180  if(rchrg>rad(mmax)) then
181    write(message, '(5a)' )&
182 &   'Pseudopotential input file core charge mesh',ch10,&
183 &   'is inconsistent with rchrg in header.',ch10,&
184 &   'Action : check your pseudopotential input file.'
185    MSG_ERROR(message)
186  end if
187 
188 !Factors for unit range scaling
189  do jj = 1, 5
190    rscale(jj)=rchrg**(jj-1)
191  end do
192 
193 !Generate uniform mesh xx in the box cut by rchrg
194 !and interpolate the core charge and derivatives
195 !Cubic polynomial interpolation is used which is consistent
196 !with the original interpolation of these functions from
197 !a log grid to the input linear grid.
198 
199  dri=1.d0/amesh
200  do i1xccc=1,n1xccc
201    xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
202 
203 !  index to find bracketing input mesh points
204    irad = int(dri * xx) + 1
205    irad = max(irad,2)
206    irad = min(irad,mmax-2)
207 !  interpolation coefficients
208    xp = dri * (xx - rad(irad))
209    xpp1 = xp + one
210    xpm1 = xp - one
211    xpm2 = xp - two
212    c1 = -xp * xpm1 * xpm2 * sixth
213    c2 = xpp1 * xpm1 * xpm2 * half
214    c3 = - xp * xpp1 * xpm2 * half
215    c4 = xp * xpp1 * xpm1 * sixth
216 !  Now do the interpolation on all derivatives for this grid point
217 !  Include 1/4pi normalization and unit range scaling
218    do jj=1,5
219      tff =  c1 * ff(irad - 1, jj) &
220 &     + c2 * ff(irad    , jj) &
221 &     + c3 * ff(irad + 1, jj) &
222 &     + c4 * ff(irad + 2, jj)
223      xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff
224    end do
225  end do
226 
227 !5th derivative is apparently not in use, so set to zero
228  xccc1d(:,6)=zero
229 
230  ABI_DEALLOCATE(ff)
231 
232 end subroutine psp9cc

TABLE OF CONTENTS

ABINIT/psp9cc [ Functions ]