TABLE OF CONTENTS

ABINIT/psp1lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp1lo

FUNCTION

 Compute sine transform to transform from v(r) to q^2 v(q)
 using subroutines related to Teter atomic structure grid.

COPYRIGHT

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

  drad(mmax)=inverse of r grid spacing at each point
  mmax=number of radial r grid points (Teter 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.
  wksincos(mmax,2,2)=contains sine and cosine of 2*pi*r(:)*dq and 2*pi*r(:)*q
    at input :  wksincos(:,1,1)=sine of 2*pi*r(:)*dq
                wksincos(:,2,1)=cosine of 2*pi*r(:)*dq
    wksincos(:,:,2) is not initialized, will be used inside the routine
  zion=nominal valence charge of atom.

OUTPUT

  epsatm= $4\pi \int[r^2 (v(r)+Zv/r) dr]$
  q2vq(mqgrid)=$q^2 v(q)$
  =$\displaystyle -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

      psp1in

CHILDREN

      der_int,sincos

SOURCE

 45 #if defined HAVE_CONFIG_H
 46 #include "config.h"
 47 #endif
 48 
 49 #include "abi_common.h"
 50 
 51 
 52 subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
 53 &  vloc,wksincos,yp1,ypn,zion)
 54 
 55  use defs_basis
 56  use m_profiling_abi
 57 
 58 !This section has been created automatically by the script Abilint (TD).
 59 !Do not modify the following lines by hand.
 60 #undef ABI_FUNC
 61 #define ABI_FUNC 'psp1lo'
 62  use interfaces_64_psp, except_this_one => psp1lo
 63 !End of the abilint section
 64 
 65  implicit none
 66 
 67 !Arguments ------------------------------------
 68 !scalars
 69  integer,intent(in) :: mmax,mqgrid
 70  real(dp),intent(in) :: zion
 71  real(dp),intent(out) :: epsatm,yp1,ypn
 72 !arrays
 73  real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax)
 74  real(dp),intent(inout) :: wksincos(mmax,2,2)
 75  real(dp),intent(out) :: q2vq(mqgrid)
 76 
 77 !Local variables-------------------------------
 78 !scalars
 79  integer,parameter :: mma0=2001
 80  integer :: iq,ir,irmax
 81  real(dp),parameter :: scale=10.0d0
 82  real(dp) :: result,test,tpiq
 83 !arrays
 84  real(dp) :: wk(mma0),wk1(mma0),wk2(mma0)
 85 
 86 ! *************************************************************************
 87 
 88 !DEBUG
 89 !write(std_out,*)' psp1lo : enter '
 90 !if(.true.)stop
 91 !ENDDEBUG
 92 
 93 !Do q=0 separately (compute epsatm)
 94 !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr]
 95 !Treat r=0 by itself
 96  wk(1)=0.0d0
 97 
 98  do ir=2,mmax
 99 !  (at large r do not want prefactor of r^2 and should see
100 !  V(r)+Zv/r go to 0 at large r)
101    test=vloc(ir)+zion/rad(ir)
102 !  DEBUG
103 !  write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test
104 !  ENDDEBUG
105 !  In this routine, NO cut-off radius is imposed : the input
106 !  vloc MUST be in real(dp) to obtain numerically
107 !  accurate values. The error can be on the order of 0.001 Ha !
108    if (abs(test)<1.0d-20) then
109      wk(ir)=0.0d0
110    else
111      wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
112    end if
113  end do
114 !Do integral from 0 to r(max) (disregard contrib beyond r(max)
115 !(need numerical derivatives to do integral)
116 !Use mmax-1 to convert to Teter s dimensioning starting at 0
117  call der_int(wk,wk2,rad,drad,mmax-1,result)
118 
119 !DEBUG
120 !write(std_out,*)' psp1lo : result ',result
121 !stop
122 !ENDDEBUG
123 
124  epsatm=4.d0*pi*(result)
125 !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi
126  q2vq(1)=-zion/pi
127 
128 !Prepare loop over q values
129  irmax=mmax+1
130  do ir=mmax,2,-1
131    test=vloc(ir)+zion/rad(ir)
132    wk1(ir)=test*rad(ir)
133 !  Will ignore tail within decade of machine precision
134    if ((scale+abs(test))==scale .and. irmax==ir+1) then
135      irmax=ir
136    end if
137  end do
138 !Increase irmax a bit : this is copied from psp1nl
139  irmax=irmax+4
140  if(irmax>mmax)irmax=mmax
141 
142 !Loop over q values
143  do iq=2,mqgrid
144    tpiq=two_pi*qgrid(iq)
145    call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
146 !  set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral
147 !$\displaystyle -Zv/\pi + q^2 4\pi \int[\frac{\sin(2\pi q r)}{2\pi q r}(r^2 v(r)+r Zv)dr]$.
148 !  Handle r=0 separately
149    wk(1)=0.0d0
150    do ir=2,irmax
151      wk(ir)=wksincos(ir,1,2)*wk1(ir)
152    end do
153 !  do integral from 0 to r(max)
154    if(irmax>mmax-1)irmax=mmax-1
155 
156    call der_int(wk,wk2,rad,drad,irmax,result)
157 !  store q^2 v(q)
158    q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result
159  end do
160 
161 !Compute derivatives of q^2 v(q) at ends of interval
162  yp1=0.0d0
163 !ypn=$\displaystyle 2\int_0^\infty (\sin (2\pi qmax r)+(2\pi qmax r)\cos (2\pi qmax r)(r V(r)+Z)dr]$
164 !integral from r(mmax) to infinity is overkill; ignore
165 !set up integrand
166 !Handle r=0 separately
167  wk(1)=0.0d0
168  tpiq=two_pi*qgrid(mqgrid)
169  do ir=2,mmax
170    test=vloc(ir)+zion/rad(ir)
171 !  Ignore contributions within decade of machine precision
172    if ((scale+abs(test))==scale) then
173      wk(ir)=0.0d0
174    else
175      wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * &
176 &     (rad(ir)*vloc(ir)+zion)
177    end if
178  end do
179  call der_int(wk,wk2,rad,drad,mmax-1,result)
180 
181  ypn=2.0d0*result
182 
183 end subroutine psp1lo