TABLE OF CONTENTS

ABINIT/psp8lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp8lo

FUNCTION

 Compute sine transform to transform from V(r) to q^2 V(q).
 Computes integrals on linear grid interpolated from the linear input
 grid with a spacing adjusted to ensure convergence at the maximum
 wavevector using corrected trapezoidal integration.

COPYRIGHT

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

  amesh=spacing for linear radial atomic grid.
  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

      psp8in,psp9in

CHILDREN

      ctrap,splfit,spline

SOURCE

 47 #if defined HAVE_CONFIG_H
 48 #include "config.h"
 49 #endif
 50 
 51 #include "abi_common.h"
 52 
 53 subroutine psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion)
 54 
 55  use defs_basis
 56  use m_splines
 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 'psp8lo'
 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) :: amesh,zion
 73  real(dp),intent(out) :: epsatm,yp1,ypn
 74 !arrays
 75  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
 76  real(dp),intent(out) :: q2vq(mqgrid)
 77 
 78 !Local variables-------------------------------
 79 !Following parameter controls accuracy of Fourier transform based on qmax
 80 !and represents the minimun number of integration points in one period.
 81 !scalars
 82  integer,parameter :: NPT_IN_2PI=200
 83  integer :: ider,iq,ir,irmu,irn,mesh_mult,mmax_new
 84  real(dp) :: amesh_new,arg,fp1,fpn,qmesh,result,ztor1
 85 !arrays
 86  real(dp),allocatable :: rad_new(:),rvlpz(:),rvlpz_new(:),sprvlpz(:,:),work(:)
 87 
 88 ! *************************************************************************
 89 
 90  ABI_ALLOCATE(work,(mmax))
 91  ABI_ALLOCATE(rvlpz,(mmax))
 92 
 93 !Do q=0 separately (compute epsatm)
 94  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
 95 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
 96  do ir=1,mmax
 97    rvlpz(ir)=rad(ir)*vloc(ir)+zion
 98    work(ir)=rad(ir)*rvlpz(ir)
 99  end do
100 
101 !Do integral from zero to r(max)
102  call ctrap(mmax,work,amesh,result)
103 
104  epsatm=4.d0*pi*result
105  q2vq(1)=-zion/pi
106 
107 !Find r mesh spacing necessary for accurate integration at qmax
108  amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
109 
110 !Choose submultiple of input mesh
111  mesh_mult=int(amesh/amesh_new) + 1
112  mmax_new=mesh_mult*(mmax-1)+1
113  amesh_new=amesh/dble(mesh_mult)
114 
115  ABI_ALLOCATE(rad_new,(mmax_new))
116  ABI_ALLOCATE(rvlpz_new,(mmax_new))
117 
118  if(mesh_mult==1) then
119    rad_new(:)=rad(:)
120    rvlpz_new(:)=rvlpz(:)
121  else
122 !  Set up spline and interpolate to finer mesh.
123 !  First, compute derivatives at end points
124    fp1=(-50.d0*rvlpz(1)+96.d0*rvlpz(2)-72.d0*rvlpz(3)+32.d0*rvlpz(4)&
125 &   -6.d0*rvlpz(5))/(24.d0*amesh)
126    fpn=(6.d0*rvlpz(mmax-4)-32.d0*rvlpz(mmax-3)+72.d0*rvlpz(mmax-2)&
127 &   -96.d0*rvlpz(mmax-1)+50.d0*rvlpz(mmax))/(24.d0*amesh)
128    ABI_ALLOCATE(sprvlpz,(mmax,2))
129    work(:)=zero
130 
131 !  Spline fit
132    call spline(rad, rvlpz,mmax,fp1,fpn,sprvlpz(:,2))
133    sprvlpz(:,1)=rvlpz(:)
134 
135 !  Set up new radial mesh
136    irn=1
137    do ir=1,mmax-1
138      do irmu=0,mesh_mult-1
139        rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
140        irn=irn+1
141      end do
142    end do
143    rad_new(mmax_new)=rad(mmax)
144 
145    ider=0
146    call splfit(rad,work,sprvlpz,ider,rad_new,rvlpz_new,mmax,mmax_new)
147 
148    ABI_DEALLOCATE(sprvlpz)
149    ABI_DEALLOCATE(work)
150    ABI_ALLOCATE(work,(mmax_new))
151  end if
152 
153 !Loop over q values
154  do iq=2,mqgrid
155    arg=2.d0*pi*qgrid(iq)
156 
157 !  Set up integrand
158    do  ir=1,mmax_new
159      work(ir)=sin(arg*rad_new(ir))*rvlpz_new(ir)
160    end do
161 
162 !  Do integral from zero to rad(mmax)
163    call ctrap(mmax_new,work,amesh_new,result)
164 
165 !  Store q^2 v(q)
166    q2vq(iq)=q2vq(1)+2.d0*qgrid(iq)*result
167 
168  end do
169 
170 !Compute derivatives of q^2 v(q) at ends of interval
171  qmesh=qgrid(2)-qgrid(1)
172  yp1=(-50.d0*q2vq(1)+96.d0*q2vq(2)-72.d0*q2vq(3)+32.d0*q2vq(4)&
173 & -6.d0*q2vq(5))/(24.d0*qmesh)
174  ypn=(6.d0*q2vq(mqgrid-4)-32.d0*q2vq(mqgrid-3)+72.d0*q2vq(mqgrid-2)&
175 & -96.d0*q2vq(mqgrid-1)+50.d0*q2vq(mqgrid))/(24.d0*qmesh)
176 
177  ABI_DEALLOCATE(work)
178  ABI_DEALLOCATE(rad_new)
179  ABI_DEALLOCATE(rvlpz_new)
180  ABI_DEALLOCATE(rvlpz)
181 
182 end subroutine psp8lo