TABLE OF CONTENTS


ABINIT/psp6cc [ Functions ]

[ Top ] [ 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