./64_psp/psp8cc.F90

ABINIT/psp8cc

ABINIT/psp8cc [ Functions ]

NAME

 psp8cc

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for format 8 of the pseudopotentials.

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

OUTPUT

  xccc1d(n1xccc,6)= 1D core charge function and its four first derivatives

PARENTS

      psp8in

CHILDREN

SOURCE

 33 #if defined HAVE_CONFIG_H
 34 #include "config.h"
 35 #endif
 36 
 37 #include "abi_common.h"
 38 
 39 subroutine psp8cc(mmax,n1xccc,rchrg,xccc1d)
 40 
 41  use defs_basis
 42  use m_errors
 43  use m_profiling_abi
 44 
 45 !This section has been created automatically by the script Abilint (TD).
 46 !Do not modify the following lines by hand.
 47 #undef ABI_FUNC
 48 #define ABI_FUNC 'psp8cc'
 49 !End of the abilint section
 50 
 51  implicit none
 52 
 53 !Arguments ------------------------------------
 54 !scalars
 55  integer,intent(in) :: mmax,n1xccc
 56  real(dp),intent(in) :: rchrg
 57 !arrays
 58  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
 59 
 60 !Local variables-------------------------------
 61 !scalars
 62  integer :: i1xccc,idum,irad,jj
 63  real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx
 64  character(len=500) :: message,errmsg
 65 !arrays
 66  real(dp) :: rscale(5)
 67  real(dp),allocatable :: ff(:,:),rad(:)
 68 
 69 !**********************************************************************
 70 
 71  ABI_ALLOCATE(ff,(mmax,5))
 72  ABI_ALLOCATE(rad,(mmax))
 73 
 74  pi4i=quarter/pi
 75 !
 76 !Read from pp file the model core charge and its first 4 derivatives
 77 !assumed to be on a linear grid starting at zero.
 78 !The input functions contain the 4pi factor, and must be rescaled.
 79 
 80  do irad=1,mmax
 81    read(tmp_unit,*, err=10, iomsg=errmsg) idum,rad(irad),(ff(irad,jj),jj=1,5)
 82  end do
 83 
 84 
 85 !Check that rad grid is linear starting at zero
 86  amesh=rad(2)-rad(1)
 87  damesh=zero
 88  do irad=2,mmax-1
 89    damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
 90  end do
 91 
 92  if(damesh>tol8 .or. rad(1)/=zero) then
 93    write(message, '(5a)' )&
 94 &   'Pseudopotential input file requires linear radial mesh',ch10,&
 95 &   'starting at zero.',ch10,&
 96 &   'Action : check your pseudopotential input file.'
 97    MSG_ERROR(message)
 98  end if
 99 
100 !Check that input rchrg is consistent with last grid point
101  if(rchrg>rad(mmax)) then
102    write(message, '(5a)' )&
103 &   'Pseudopotential input file core charge mesh',ch10,&
104 &   'is inconsistent with rchrg in header.',ch10,&
105 &   'Action : check your pseudopotential input file.'
106    MSG_ERROR(message)
107  end if
108 
109 !Factors for unit range scaling
110  do jj = 1, 5
111    rscale(jj)=rchrg**(jj-1)
112  end do
113 
114 !Generate uniform mesh xx in the box cut by rchrg
115 !and interpolate the core charge and derivatives
116 !Cubic polynomial interpolation is used which is consistent
117 !with the original interpolation of these functions from
118 !a log grid to the input linear grid.
119 
120  dri=1.d0/amesh
121  do i1xccc=1,n1xccc
122    xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
123 
124 !  index to find bracketing input mesh points
125    irad = int(dri * xx) + 1
126    irad = max(irad,2)
127    irad = min(irad,mmax-2)
128 !  interpolation coefficients
129    xp = dri * (xx - rad(irad))
130    xpp1 = xp + one
131    xpm1 = xp - one
132    xpm2 = xp - two
133    c1 = -xp * xpm1 * xpm2 * sixth
134    c2 = xpp1 * xpm1 * xpm2 * half
135    c3 = - xp * xpp1 * xpm2 * half
136    c4 = xp * xpp1 * xpm1 * sixth
137 !  Now do the interpolation on all derivatives for this grid point
138 !  Include 1/4pi normalization and unit range scaling
139    do jj=1,5
140      tff =  c1 * ff(irad - 1, jj) &
141 &     + c2 * ff(irad    , jj) &
142 &     + c3 * ff(irad + 1, jj) &
143 &     + c4 * ff(irad + 2, jj)
144      xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff
145    end do
146  end do
147 
148 !5th derivative is apparently not in use, so set to zero
149  xccc1d(:,6)=zero
150 
151  ABI_DEALLOCATE(ff)
152  ABI_DEALLOCATE(rad)
153 
154  return 
155 
156  ! Handle IO error
157  10 continue
158  MSG_ERROR(errmsg)
159 
160 end subroutine psp8cc

TABLE OF CONTENTS

ABINIT/psp8cc [ Functions ]