TABLE OF CONTENTS

ABINIT/psp5lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp5lo

FUNCTION

 Compute sine transform to transform from V(r) to q^2 V(q).
 Computes integrals on logarithmic grid using related uniform
 grid in exponent and corrected trapezoidal integration.

COPYRIGHT

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

  al=spacing in exponent for radial atomic grid.
  mmax=number of radial r grid points (logarithmic atomic grid).
  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=q grid values (bohr**-1).
  rad(mmax)=r grid values (bohr).
  vloc(mmax)=V(r) on radial grid.
  zion=nominal valence charge of atom.

OUTPUT

  epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
  q2vq(mqgrid)
   =q^2 V(q)
   = -\frac{Zv}{\pi}
     + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
  yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax
   (needed for spline fitter).

PARENTS

      psp5in,psp6in

CHILDREN

      ctrap

SOURCE

 46 #if defined HAVE_CONFIG_H
 47 #include "config.h"
 48 #endif
 49 
 50 #include "abi_common.h"
 51 
 52 subroutine psp5lo(al,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
 53 &                  vloc,yp1,ypn,zion)
 54 
 55  use defs_basis
 56  use m_profiling_abi
 57 
 58  use m_numeric_tools, only : ctrap
 59 
 60 !This section has been created automatically by the script Abilint (TD).
 61 !Do not modify the following lines by hand.
 62 #undef ABI_FUNC
 63 #define ABI_FUNC 'psp5lo'
 64 !End of the abilint section
 65 
 66  implicit none
 67 
 68 !Arguments----------------------------------------------------------
 69 !scalars
 70  integer,intent(in) :: mmax,mqgrid
 71  real(dp),intent(in) :: al,zion
 72  real(dp),intent(out) :: epsatm,yp1,ypn
 73 !arrays
 74  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
 75  real(dp),intent(out) :: q2vq(mqgrid)
 76 
 77 !Local variables-------------------------------
 78 !scalars
 79  integer :: iq,ir
 80  real(dp),parameter :: scale=10.0d0
 81  real(dp) :: arg,result,rmtoin,test,ztor1
 82 !arrays
 83  real(dp),allocatable :: work(:)
 84 
 85 ! *************************************************************************
 86 
 87  ABI_ALLOCATE(work,(mmax))
 88 
 89 !Do q=0 separately (compute epsatm)
 90 !Do integral from 0 to r1
 91  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
 92 
 93 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
 94 !with extra factor of r to convert to uniform grid in exponent
 95  do ir=1,mmax
 96 !  First handle tail region
 97    test=vloc(ir)+zion/rad(ir)
 98 !  DEBUG
 99 !  write(std_out,*)ir,rad(ir),test
100 !  ENDDEBUG
101 !  Ignore small contributions, or impose a cut-off in the case
102 !  the pseudopotential data are in single precision.
103 !  (it is indeed expected that vloc is very close to zero beyond 20,
104 !  so a value larger than 2.0d-8 is considered anomalous)
105    if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
106      work(ir)=zero
107    else
108      work(ir)=(rad(ir)*rad(ir))*(rad(ir)*vloc(ir)+zion)
109    end if
110  end do
111 !DEBUG
112 !write(std_out,*)' psp5lo : stop '
113 !stop
114 !ENDDEBUG
115 
116 !Do integral from r(1) to r(max)
117  call ctrap(mmax,work,al,result)
118 !Do integral from r(mmax) to infinity
119 !compute decay length lambda at r(mmax)
120 !$\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ &
121 !$(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$
122 !rmtoin=$(rad(mmax)*vloc(mmax)+zion)*(rad(mmax)+1.d0/\lambda)/\lambda$
123 !Due to inability to fit exponential decay to r*V(r)+Zv
124 !in tail, NO TAIL CORRECTION IS APPLIED
125 !(numerical trouble might be removed if atomic code is
126 !cleaned up in tail region)
127  rmtoin=0.0d0
128 
129  epsatm=4.d0*pi*(result+ztor1+rmtoin)
130 
131  q2vq(1)=-zion/pi
132 
133 !Loop over q values
134  do iq=2,mqgrid
135    arg=2.d0*pi*qgrid(iq)
136 !  ztor1=$ -Zv/\pi+2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
137    ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion)* &
138 &   cos(arg*rad(1)) )/pi
139 
140 !  set up integrand
141    do  ir=1,mmax
142      test=vloc(ir)+zion/rad(ir)
143 !    Ignore contributions within decade of machine precision
144      if ((scale+abs(test)).eq.scale) then
145        work(ir)=zero
146      else
147        work(ir)=rad(ir)*sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
148      end if
149    end do
150 !  do integral from r(1) to r(mmax)
151    call ctrap(mmax,work,al,result)
152 
153 !  do integral from r(mmax) to infinity
154 !  rmtoin=(r(mmax)*vr(mmax)+zion)*(lambda*sin(arg*r(mmax))+
155 !  arg*cos(arg*r(mmax)))/(arg**2+lambda**2)
156 !  See comment above; no tail correction
157    rmtoin=0.0d0
158 
159 !  store q^2 v(q)
160    q2vq(iq)=ztor1+2.d0*qgrid(iq)*(result+rmtoin)
161 
162  end do
163 
164 !Compute derivatives of q^2 v(q) at ends of interval
165  yp1=0.0d0
166 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
167 !integral from 0 to r1
168  arg=2.0d0*pi*qgrid(mqgrid)
169  ztor1=zion*rad(1)*sin(arg*rad(1))
170  ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + &
171 & (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1))
172 !integral from r(mmax) to infinity is overkill; ignore
173 !set up integrand
174  do ir=1,mmax
175    test=vloc(ir)+zion/rad(ir)
176 !  Ignore contributions within decade of machine precision
177    if ((scale+abs(test)).eq.scale) then
178      work(ir)=0.0d0
179    else
180      work(ir)=rad(ir)*(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * &
181 &     (rad(ir)*vloc(ir)+zion)
182    end if
183  end do
184  call ctrap(mmax,work,al,result)
185  ypn=2.0d0 * (ztor1 + result)
186 
187  ABI_DEALLOCATE(work)
188 
189 end subroutine psp5lo