TABLE OF CONTENTS

ABINIT/psp11lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp11lo

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.
 Generalized from psp5lo for non log grids using dr / di

COPYRIGHT

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

  drdi=derivative of radial grid wrt index
  mmax=number of radial r grid points
  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

      upf2abinit

CHILDREN

      ctrap

SOURCE

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