TABLE OF CONTENTS


ABINIT/m_psp1 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psp1

FUNCTION

  Initialize pspcod=1 or 4 pseudopotential (Teter format)

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, FrD, MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_psp1
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27  use m_splines
28 
29  use m_special_funcs,   only : besjm
30  use m_psptk,           only : psp1cc
31 
32  implicit none
33 
34  private

m_psp1/der_int [ Functions ]

[ Top ] [ m_psp1 ] [ Functions ]

NAME

 der_int

FUNCTION

 Given input function f(i) on Teter radial grid, and grid spacing
 dr(i), compute function derivative df/dr on points from 0 to n.
 Integrate function f(i) on grid r(i) from r(0) to r(nlast).
 Note that array dimensions start at 0.

INPUTS

  f(0 to nlast)=function values on grid
  r(0 to nlast)=radial grid points
  dr(0 to nlast)=INVERSE of spacing on grid
  nlast=radial grid point for upper limit

OUTPUT

  df(0 to n)=derivative $ \frac{df}{dr}$ on grid
  smf= $ \int_{r(0)}^{r(nlast)} f(r) dr $.

SOURCE

853 subroutine der_int(ff,df,rr,dr,nlast,smf)
854 
855 !Arguments ------------------------------------
856 !nmax sets standard number of grid points ! SHOULD BE REMOVED
857 !scalars
858  integer,parameter :: nmax=2000
859  integer,intent(in) :: nlast
860  real(dp),intent(out) :: smf
861 !no_abirules
862 !Note that dimension here starts at 0
863  real(dp), intent(in) :: dr(0:nmax),ff(0:nmax),rr(0:nmax)
864  real(dp), intent(out) :: df(0:nmax)
865 
866 !Local variables-------------------------------
867 !scalars
868  integer :: jj
869  real(dp),parameter :: div12=1.d0/12.d0
870  real(dp) :: hh
871  character(len=500) :: message
872 
873 ! *************************************************************************
874 
875 !Check that nlast lie within 0 to nmax
876  if (nlast<0.or.nlast>nmax) then
877    write(message, '(a,i12,a,i12)' )&
878 &   ' nlast=',nlast,' lies outside range [0,nmax] with dimension nmax=',nmax
879    ABI_BUG(message)
880  end if
881 
882 !Compute derivatives at lower end, near r=0
883  df(0)=-25.d0/12.d0*ff(0)+4.d0*ff(1)-3.d0*ff(2)+4.d0/3.d0*ff(3)&
884 & -1.d0/4.d0*ff(4)
885  df(1)=-1.d0/4.d0*ff(0)-5.d0/6.d0*ff(1)+3.d0/2.d0*ff(2)&
886 & -1.d0/2.d0*ff(3)+1.d0/12.d0*ff(4)
887 
888 !Run over range from just past r=0 to near r(n), using central differences
889  do jj=2,nlast-2
890    df(jj)=(ff(jj-2)-8.d0*(ff(jj-1)-ff(jj+1))-ff(jj+2))*div12
891  end do
892 
893 !Compute derivative at upper end of range
894  if (nlast < 4) then
895    message = ' der_int: ff does not have enough elements. nlast is too low'
896    ABI_ERROR(message)
897  end if
898 
899  df(nlast-1)=-1.d0/12.d0*ff(nlast-4)&
900 & +1.d0/2.d0*ff(nlast-3)&
901 & -3.d0/2.d0*ff(nlast-2)&
902 & +5.d0/6.d0*ff(nlast-1)&
903 & +1.d0/4.d0*ff(nlast)
904  df(nlast)=1.d0/4.d0*ff(nlast-4)&
905 & -4.d0/3.d0*ff(nlast-3)&
906 & +3.d0*ff(nlast-2)&
907 & -4.d0*ff(nlast-1)&
908 & +25.d0/12.d0*ff(nlast)
909 
910 !Apply correct normalization over full range
911  do jj=0,nlast
912    df(jj)=df(jj)*dr(jj)
913  end do
914 
915  smf=0.d0
916  do jj=0,nlast-1
917    hh=rr(jj+1)-rr(jj)
918    smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1)))
919  end do
920  smf=smf/12.d0
921 
922 end subroutine der_int

m_psp1/psp1in [ Functions ]

[ Top ] [ m_psp1 ] [ Functions ]

NAME

 psp1in

FUNCTION

 Initialize pspcod=1 or 4 pseudopotential (Teter format):
 continue to read the corresponding file, then compute the
 local and non-local potentials.

INPUTS

  dq= spacing of the q-grid
  lloc=angular momentum choice of local pseudopotential
  lmax=value of lmax mentioned at the second line of the psp file
  lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
        =if useylm=0, max number of (l,n)   comp. over all type of psps
  lnmax=max. number of (l,n) components over all type of psps
  mmax=maximum number of points in real space grid in the psp file
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=dimension of q (or G) grid for arrays.
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  pspcod=pseudopotential type
  qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
  useylm=governs the way the nonlocal operator is to be applied:
         1=using Ylm, 0=using Legendre polynomials
  zion=nominal valence of atom as specified in psp file
  znucl=atomic number of atom as specified in psp file

OUTPUT

  ekb(lnmax)=Kleinman-Bylander energy,
             {{\ \begin{equation}
               \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
             {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))   dr]}
              \end{equation} }}
             for each (l,n)
  ekb1(mpsang)= Kleinman-Bylander energy from the psp file, for iproj=1
  ekb2(mpsang)= Kleinman-Bylander energy from the psp file, for iproj=2
  epsatm=$(4\pi) \int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
  epspsp(mpsang)=values of epsatm for different angular momenta, from the psp file
  e990(mpsang)=ecut at which 0.99 of the kinetic energy is recovered
  e999(mpsang)=ecut at which 0.999 of the kinetic energy is recovered
  ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum and
   each projector
  indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln  (if useylm=0)
                                           or i=lmn (if useylm=1)
  nproj(mpsang)=number of projection functions for each angular momentum
  qchrg is the total (integrated) core charge
  rcpsp(mpsang)=cut-off radius for each angular momentum
  rms(mpsang)=root mean square of the KB psp
  vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit
  xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
  xcccrc=XC core correction cutoff radius (bohr)

NOTES

 there are only minor differences in the two formats
 1) With pspcod=1, even for the LOCAL angular momentum, there is
    a block for the wfs (can be set to zero, though)
 2) The core charge density differs: for pspcod=1, it is a
    revised expression for core density of 5 Nov 1992, while
    for pspcod=4, it is an older expression, of 7 May 1992 .

SOURCE

106 subroutine psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,&
107 &                  e990,e999,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
108 &                  mmax,mpsang,mqgrid,nproj,n1xccc,pspcod,&
109 &                  qchrg,qgrid,rcpsp,rms,useylm,vlspl,xcccrc,xccc1d,&
110 &                  zion,znucl)
111 
112 !Arguments ------------------------------------
113 !scalars
114  integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc,pspcod
115  integer,intent(in) :: useylm
116  real(dp),intent(in) :: dq,zion,znucl
117  real(dp),intent(out) :: epsatm,qchrg,xcccrc
118 !arrays
119  integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang)
120  real(dp),intent(in) :: qgrid(mqgrid)
121  real(dp),intent(out) :: e990(mpsang),e999(mpsang),ekb(lnmax),ekb1(mpsang)
122  real(dp),intent(out) :: ekb2(mpsang),epspsp(mpsang)
123  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax)
124  real(dp),intent(out) :: rcpsp(mpsang),rms(mpsang),vlspl(mqgrid,2)
125  real(dp),intent(inout) :: xccc1d(n1xccc,6)
126 
127 !Local variables-------------------------------
128 !scalars
129  integer :: ii,iln,index,ipsang,kk,lhigh,ll,mm,nlmax
130  real(dp) :: arg,dq2pi,fchrg,rchrg,xx,yp1,ypn
131  character(len=500) :: message,errmsg
132 !arrays
133  real(dp),allocatable :: drad(:),ekb_tmp(:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:)
134  real(dp),allocatable :: vpspll(:,:),wfll(:,:),wksincos(:,:,:),work_space(:)
135  real(dp),allocatable :: work_spl1(:),work_spl2(:)
136 
137 ! ***************************************************************************
138 
139 !Note: Teter s grid is hard-coded at mmax=2001
140 !mmax was read from the pseudopotential file in the calling routine
141  if (mmax/=2001) then
142    write(message, '(a,i12,a,a,a,a)' )&
143 &   'Using Teter grid (pspcod=1 and 4) but mmax=',mmax,ch10,&
144 &   'mmax must be 2001 for Teter grid.',ch10,&
145 &   'Action: check your pseudopotential input file.'
146    ABI_ERROR(message)
147  end if
148 
149 !File format of formatted Teter psp input (the 3 first lines
150 !have already been read in calling -pspatm- routine) :
151 
152 !(1) title (character) line
153 !(2) znucl,zion,pspdat
154 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well
155 !For each angular momentum :
156 !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll)
157 !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll)
158 !(6) rchrg,fchrg,qchrg
159 !(7) ll
160 !(8) (vpsp(j,ll),j=0,nmax)
161 !Then for iproj=1 to 2
162 !for ll=0,lmax
163 !(10) ll
164 !(11) ((upsp(j,ll,iproj),j=0,nmax)
165 
166  do ipsang=1,lmax+1
167 
168    read (tmp_unit,*,err=10,iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang)
169    write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) &
170 &   ipsang-1,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),&
171 &   'l,e99.0,e99.9,nproj,rcpsp'
172    call wrtout(ab_out,message,'COLL')
173    call wrtout(std_out,  message,'COLL')
174 
175    read (tmp_unit,*,err=10,iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang)
176    write(message, '(4f13.8,t55,a)' ) &
177 &   rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),&
178 &   '   rms, ekb1, ekb2, epsatm'
179    call wrtout(ab_out,message,'COLL')
180    call wrtout(std_out,  message,'COLL')
181 
182  end do
183 
184 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn
185  index=0;iln=0;indlmn(:,:)=0
186  do ipsang=1,lmax+1
187    if(nproj(ipsang)>0)then
188      ll=ipsang-1
189      do kk=1,nproj(ipsang)
190        iln=iln+1
191        do mm=1,2*ll*useylm+1
192          index=index+1
193          indlmn(1,index)=ll
194          indlmn(2,index)=mm-ll*useylm-1
195          indlmn(3,index)=kk
196          indlmn(4,index)=ll*ll+(1-useylm)*ll+mm
197          indlmn(5,index)=iln
198          indlmn(6,index)=1
199        end do
200      end do
201    end if
202  end do
203 
204  read (tmp_unit,*,err=10,iomsg=errmsg) rchrg,fchrg,qchrg
205  write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
206  call wrtout(ab_out,message,'COLL')
207  call wrtout(std_out,  message,'COLL')
208 
209  ! Generate core charge function and derivatives, if needed
210  if(fchrg>1.0d-15)then
211    if(pspcod==1)then
212      call psp1cc(fchrg,n1xccc,xccc1d)
213      ! The core charge function for pspcod=1 becomes zero beyond 3*rchrg only.
214      ! Thus xcccrc must be set equal to 3*rchrg .
215      xcccrc=3*rchrg
216    else if(pspcod==4)then
217      call psp4cc(fchrg,n1xccc,xccc1d)
218      ! For pspcod=4, the core charge cut off exactly beyond rchrg
219      xcccrc=rchrg
220    end if
221  else
222    xcccrc=0.0d0
223    xccc1d(:,:)=0.0d0
224  end if
225 
226 !--------------------------------------------------------------------
227 !Will now proceed at the reading of pots and wfs, as well as their treatment
228 
229 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials
230 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg.
231 !rad(:)=radial grid r(i)
232 !drad(:)= inverse of d(r(i))/d(i) for radial grid
233 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions
234 
235  ABI_MALLOC(vloc,(mmax))
236  ABI_MALLOC(vpspll,(mmax,mpsang))
237  if(lmax==-1) vpspll(:,:)=zero
238 
239 !(1) Read atomic pseudopotential for each l, filling up array vpspll
240 !Note: put each l into vpspll(:,l+1)
241  do ipsang=1,lmax+1
242    read (tmp_unit,*,err=10,iomsg=errmsg) ll
243    read (tmp_unit,*,err=10,iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax)
244  end do
245 
246 !Copy appropriate nonlocal psp for use as local one
247  vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 )
248 
249 !(2) Create radial grid, and associated quantities
250  ABI_MALLOC(rad,(mmax))
251  ABI_MALLOC(drad,(mmax))
252  ABI_MALLOC(wksincos,(mmax,2,2))
253 
254 !Teter grid--need both r and dr in this case
255  do ii=0,mmax-1
256    xx=dble(ii)/dble(mmax-1)
257    rad (ii+1)=100.d0*(xx+.01d0)**5-1.d-8
258    drad(ii+1)=500.d0*(xx+.01d0)**4/dble(mmax-1)
259  end do
260 
261 !here compute sin(r(:)*dq) and cos(r(:)*dq)
262 !NOTE: also invert dr !!
263  dq2pi=2.0d0*pi*dq
264  do ii=1,mmax
265    arg=dq2pi*rad(ii)
266    drad(ii)=1.0d0/drad(ii)
267    wksincos(ii,1,1)=sin(arg)
268    wksincos(ii,2,1)=cos(arg)
269  end do
270 
271 !(3)Carry out calculations for local (lloc) pseudopotential.
272 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q).
273  ABI_MALLOC(work_space,(mqgrid))
274  ABI_MALLOC(work_spl1,(mqgrid))
275  ABI_MALLOC(work_spl2,(mqgrid))
276  call psp1lo(drad,epsatm,mmax,mqgrid,qgrid,&
277 & work_spl1,rad,vloc,wksincos,yp1,ypn,zion)
278 
279 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
280  call spline (qgrid,work_spl1,mqgrid,yp1,ypn,work_spl2)
281  vlspl(:,1)=work_spl1(:)
282  vlspl(:,2)=work_spl2(:)
283 
284  ABI_FREE(work_space)
285  ABI_FREE(work_spl1)
286  ABI_FREE(work_spl2)
287 
288 !(4)Take care of non-local part
289 
290 !Zero out all Kleinman-Bylander energies to initialize
291  ekb(:)=0.0d0
292 !write(std_out,*)' psp1in : before nonlocal corrections '
293 !write(std_out,*)' psp1in : lloc, lmax = ',lloc,lmax
294 
295 !Allow for option of no nonlocal corrections (lloc=lmax=0)
296  if (lloc==0.and.lmax==0) then
297    write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
298    call wrtout(ab_out,message,'COLL')
299    call wrtout(std_out,  message,'COLL')
300 
301  else
302 
303 !  Proceed to make Kleinman-Bylander form factors for each l up to lmax
304 
305 !  Read wavefunctions for each l up to lmax
306    ABI_MALLOC(wfll,(mmax,mpsang))
307    do ipsang=1,lmax+1
308 !    For pspcod==4, wfs for the local angular momentum are not written
309      if (nproj(ipsang)/=0 .or. pspcod==1) then
310        read (tmp_unit,*,err=10,iomsg=errmsg) ll
311        if (ipsang/=ll+1) then
312          write(message, '(a,a,a,a,a,a,2i6,a,a)' )&
313 &         'Pseudopotential input file does not have',ch10,&
314 &         'angular momenta in order expected for first projection',&
315 &         'operator.',ch10,' Values are ',ipsang-1,ll,ch10,&
316 &         'Action: check your pseudopotential input file.'
317          ABI_ERROR(message)
318        end if
319        read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang)
320 
321      else
322        wfll(:,ipsang)=0.0d0
323      end if
324 
325    end do
326 !  ----------------------------------------------------------------------
327 !  Compute KB form factors and fit splines
328 
329 !  nlmax is highest l for which a nonlocal correction is being computed
330    nlmax=lmax
331    if (lloc==lmax) nlmax=lmax-1
332 !  write(std_out,*)' psp1in : lmax,lloc=',lmax,lloc
333    ABI_MALLOC(ekb_tmp,(mpsang,2))
334    ABI_MALLOC(ffspl_tmp,(mqgrid,2,nlmax+1,2))
335 
336    call psp1nl(drad,ekb_tmp(:,1),ffspl_tmp(:,:,:,1),lloc,&
337 &   nlmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos)
338 
339 !  Read second wavefunction for second projection operator
340 !  (only read cases where nproj(ll)=2) --also find highest l for which nproj(l)=2
341    lhigh=-1
342    do ipsang=1,min(lmax+1,mpsang)
343      if (nproj(ipsang)==2) then
344        lhigh=ipsang-1
345        read (tmp_unit,*,err=10,iomsg=errmsg) ll
346        if (ipsang/=ll+1) then
347          write(message, '(a,a,a,a,a,a,2i6,a,a)' )&
348 &         'Pseudopotential input file does not have',ch10,&
349 &         'angular momenta in order expected for second projection',&
350 &         'operator.',ch10,' Values are ',ipsang-1,ll,ch10,&
351 &         'Action: check your pseudopotential input file.'
352          ABI_ERROR(message)
353        end if
354        read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang)
355 
356      else
357        wfll(:,ipsang)=0.0d0
358 
359      end if
360    end do
361 
362 !  Compute KB form factors and fit splines for second wf if any
363 
364    if (lhigh>-1) then
365      call psp1nl(drad,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lloc,&
366 &     lhigh,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos)
367    end if
368 
369 !  Convert ekb and ffspl
370    iln=0
371    do ii=1,lmnmax
372      kk=indlmn(5,ii)
373      if (kk>iln) then
374        iln=kk
375        ekb(kk)=ekb_tmp(1+indlmn(1,ii),indlmn(3,ii))
376 !      write(std_out,*)' psp1in : lmnmax,ii,indlmn(1,ii)=',lmnmax,ii,indlmn(1,ii)
377        ffspl(:,:,kk)=ffspl_tmp(:,:,1+indlmn(1,ii),indlmn(3,ii))
378      end if
379    end do
380 
381    ABI_FREE(ekb_tmp)
382    ABI_FREE(ffspl_tmp)
383    ABI_FREE(wfll)
384  end if
385 
386  ABI_FREE(vpspll)
387  ABI_FREE(rad)
388  ABI_FREE(drad)
389  ABI_FREE(vloc)
390  ABI_FREE(wksincos)
391 
392  return
393 
394  ! Handle IO error
395  10 continue
396  ABI_ERROR(errmsg)
397 
398 end subroutine psp1in

m_psp1/psp1lo [ Functions ]

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

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).

SOURCE

431 subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
432 &  vloc,wksincos,yp1,ypn,zion)
433 
434 !Arguments ------------------------------------
435 !scalars
436  integer,intent(in) :: mmax,mqgrid
437  real(dp),intent(in) :: zion
438  real(dp),intent(out) :: epsatm,yp1,ypn
439 !arrays
440  real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax)
441  real(dp),intent(inout) :: wksincos(mmax,2,2)
442  real(dp),intent(out) :: q2vq(mqgrid)
443 
444 !Local variables-------------------------------
445 !scalars
446  integer,parameter :: mma0=2001
447  integer :: iq,ir,irmax
448  real(dp),parameter :: scale=10.0d0
449  real(dp) :: result,test,tpiq
450 !arrays
451  real(dp) :: wk(mma0),wk1(mma0),wk2(mma0)
452 
453 ! *************************************************************************
454 
455 !Do q=0 separately (compute epsatm)
456 !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr]
457 !Treat r=0 by itself
458  wk(1)=0.0d0
459 
460  do ir=2,mmax
461 !  (at large r do not want prefactor of r^2 and should see
462 !  V(r)+Zv/r go to 0 at large r)
463    test=vloc(ir)+zion/rad(ir)
464 !  write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test
465 !  In this routine, NO cut-off radius is imposed : the input
466 !  vloc MUST be in real(dp) to obtain numerically
467 !  accurate values. The error can be on the order of 0.001 Ha !
468    if (abs(test)<1.0d-20) then
469      wk(ir)=0.0d0
470    else
471      wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
472    end if
473  end do
474 !Do integral from 0 to r(max) (disregard contrib beyond r(max)
475 !(need numerical derivatives to do integral)
476 !Use mmax-1 to convert to Teter s dimensioning starting at 0
477  call der_int(wk,wk2,rad,drad,mmax-1,result)
478 
479  epsatm=4.d0*pi*(result)
480 !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi
481  q2vq(1)=-zion/pi
482 
483 !Prepare loop over q values
484  irmax=mmax+1
485  do ir=mmax,2,-1
486    test=vloc(ir)+zion/rad(ir)
487    wk1(ir)=test*rad(ir)
488 !  Will ignore tail within decade of machine precision
489    if ((scale+abs(test))==scale .and. irmax==ir+1) then
490      irmax=ir
491    end if
492  end do
493 !Increase irmax a bit : this is copied from psp1nl
494  irmax=irmax+4
495  if(irmax>mmax)irmax=mmax
496 
497 !Loop over q values
498  do iq=2,mqgrid
499    tpiq=two_pi*qgrid(iq)
500    call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
501 !  set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral
502 !$\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]$.
503 !  Handle r=0 separately
504    wk(1)=0.0d0
505    do ir=2,irmax
506      wk(ir)=wksincos(ir,1,2)*wk1(ir)
507    end do
508 !  do integral from 0 to r(max)
509    if(irmax>mmax-1)irmax=mmax-1
510 
511    call der_int(wk,wk2,rad,drad,irmax,result)
512 !  store q^2 v(q)
513    q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result
514  end do
515 
516 !Compute derivatives of q^2 v(q) at ends of interval
517  yp1=0.0d0
518 !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]$
519 !integral from r(mmax) to infinity is overkill; ignore
520 !set up integrand
521 !Handle r=0 separately
522  wk(1)=0.0d0
523  tpiq=two_pi*qgrid(mqgrid)
524  do ir=2,mmax
525    test=vloc(ir)+zion/rad(ir)
526 !  Ignore contributions within decade of machine precision
527    if ((scale+abs(test))==scale) then
528      wk(ir)=0.0d0
529    else
530      wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * &
531 &     (rad(ir)*vloc(ir)+zion)
532    end if
533  end do
534  call der_int(wk,wk2,rad,drad,mmax-1,result)
535 
536  ypn=2.0d0*result
537 
538 end subroutine psp1lo

m_psp1/psp1nl [ Functions ]

[ Top ] [ m_psp1 ] [ Functions ]

NAME

 psp1nl

FUNCTION

 Make Kleinman-Bylander form factors f_l(q) for each l from 0 to lmax.
 Vloc is assumed local potential.

INPUTS

  dr(mmax)=inverse of grid spacing for radial grid
  lloc=angular momentum of local channel (avoid doing integrals for this l)
  lmax=maximum ang momentum for which nonlocal form factor is desired.
  mmax=number of radial grid points for atomic grid
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=number of grid points for q grid
  qgrid(mqgrid)=values at which form factors are returned
  rad(mmax)=radial grid values
  vloc(mmax)=local pseudopotential on radial grid
  vpspll(mmax,lmax+1)=nonlocal pseudopotentials for each l on radial grid
  wfll(mmax,lmax+1)=reference state wavefunctions 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

OUTPUT

  ekb(mpsang)=Kleinman-Bylander energy,
              {{\ \begin{equation}
               \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
              {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))   dr]}
               \end{equation} }}
              for each l
  ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum

NOTES

 u_l(r) is reference state wavefunction (input as wfll);
 j_l(q) is a spherical Bessel function;
 dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
 f_l(q) =$ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
 where dvms=$\displaystyle \int_0^{rmax}[(u_l(r) dV_l(r))^2 dr]$ is the mean
 square value of the nonlocal correction for angular momentum l.
 E_KB = $\displaystyle \frac{dvms}{\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]}$.
 This is the eigenvalue of the Kleinman-Bylander operator and sets
 the energy scale of the nonlocal psp corrections.
 Bessel functions replaced by besj, which accomodates args near 0.

SOURCE

590 subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,&
591 &                  qgrid,rad,vloc,vpspll,wfll,wksincos)
592 
593 !Arguments ------------------------------------
594 !scalars
595  integer,intent(in) :: lloc,lmax,mmax,mpsang,mqgrid
596 !arrays
597  real(dp),intent(in) :: dr(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax)
598  real(dp),intent(in) :: vpspll(mmax,mpsang),wfll(mmax,mpsang)
599  real(dp),intent(inout) :: wksincos(mmax,2,2)
600  real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang)
601 
602 !Local variables-------------------------------
603 !scalars
604  integer,parameter :: dpsang=5
605  integer :: iq,ir,irmax,lp1
606  real(dp) :: dvwf,result,test,tpiq,yp1,ypn
607  character(len=500) :: message
608 !arrays
609  real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang)
610  real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:)
611  real(dp),allocatable :: work_spl(:)
612 
613 ! *************************************************************************
614 
615 !Zero out Kleinman-Bylander energies ekb
616  ekb(:)=0.0d0
617 !Zero out eta and other parameters too (so 0 s show up in output later)
618  eta(:)=0.0d0
619  dvms(:)=0.0d0
620  ckb(:)=0.0d0
621 
622 !Allow for no nonlocal correction (lmax=-1)
623  if (lmax/=-1) then
624 
625 !  Check that lmax is within allowed range
626    if (lmax<0.or.lmax>3) then
627      write(message, '(a,i12,a,a,a,a,a,a,a)' )&
628 &     'lmax=',lmax,' is not an allowed value.',ch10,&
629 &     'Allowed values are -1 for no nonlocal correction or else',ch10,&
630 &     '0, 1, 2, or 3 for maximum l nonlocal correction.',ch10,&
631 &     'Action: check the input atomic psp data file for lmax.'
632      ABI_ERROR(message)
633    end if
634 
635 !  Compute normalizing integrals eta=<dV> and mean square
636 !  nonlocal psp correction dvms=<dV^2>
637 !  "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction
638 
639    ABI_MALLOC(work1,(mmax+1))
640    ABI_MALLOC(work2,(mmax+1))
641    ABI_MALLOC(work_spl,(mqgrid))
642    ABI_MALLOC(work5,(mmax))
643    ABI_MALLOC(besjx,(mmax))
644 
645    do lp1=1,lmax+1
646 
647 !    Only do the work if nonlocal correction is nonzero
648      if (lp1 /= lloc+1) then
649 
650 !      integrand for 0 to r(mmax)
651        do ir=1,mmax
652          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
653          work1(ir)=wfll(ir,lp1)*dvwf
654        end do
655 
656 !      do integral
657 !      first need derivative of function; note use of
658 !      shifted indices to accomodate Mike Teter s choice of 0:mmax-1
659        call der_int(work1,work2,rad,dr,mmax-1,result)
660        eta(lp1)=result
661 
662 !      DEBUG
663 !      write(std_out,*)' psp1nl : write eta(lp1)'
664 !      write(std_out,*)result
665 !      do ir=1,mmax,61
666 !      write(std_out,*)vpspll(ir,lp1),vloc(ir),wfll(ir,lp1)
667 !      end do
668 !      write(std_out,*)
669 !      do ir=1,mmax,61
670 !      write(std_out,*)work1(ir),rad(ir),dr(ir)
671 !      end do
672 !      ENDDEBUG
673 
674        do ir=1,mmax
675          dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
676          work1(ir)=dvwf**2
677        end do
678        call der_int(work1,work2,rad,dr,mmax-1,result)
679 
680        dvms(lp1)=result
681 
682 !      If dvms is not 0 for any given angular momentum l,
683 !      compute Xavier Gonze s definition of the Kleinman-Bylander
684 !      energy E_KB = dvms/eta.  In this case also renormalize
685 !      the projection operator to u_KB(r)=$u_l(r) dV(r)/\sqrt{dvms}$.
686 !      This means dvwf gets multiplied by the normalization factor
687 !      "renorm"=$1/\sqrt{dvms}$ as seen below.
688 !      With dvwf=dV(r)*wf(r) for wf(r)=``radial'' wf, the integrand
689 !      for each angular momentum l is
690 !      Bessel_l(2 $\pi$ q r) * wf(r) * dV(r) * r;
691 !      NOTE presence of extra r in integrand.
692 
693        if (dvms(lp1)/=0.0d0) then
694          ekb(lp1)=dvms(lp1)/eta(lp1)
695          renorm(lp1)=1.0d0/sqrt(dvms(lp1))
696 !        ckb is Kleinman-Bylander "cosine" (Xavier Gonze)
697          ckb(lp1)=eta(lp1)/sqrt(dvms(lp1))
698        else
699          ekb(lp1)=0.0d0
700        end if
701      end if
702    end do
703 
704 !  Loop on angular momenta
705    do lp1=1,lmax+1
706 
707 !    Compute form factor if ekb(lp1) not 0
708      if (ekb(lp1)/=0.0d0) then
709 
710 !      do q=0 separately, non-zero if l=0
711        if(lp1==1)then
712          do ir=1,mmax
713            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
714            work1(ir)=rad(ir)*dvwf
715          end do
716          call der_int(work1,work2,rad,dr,mmax-1,result)
717          ffspl(1,1,lp1)=result
718        else
719 !        For l non-zero, f(q=0) vanishes !
720          ffspl(1,1,lp1)=0.0d0
721        end if
722 
723 !      Prepare loop over q values
724        irmax=mmax+1
725        do ir=mmax,2,-1
726          test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir)
727          work5(ir)=test
728          work1(ir)=0.0d0
729 !        Will ignore tail within decade of machine precision
730          if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then
731            irmax=ir
732          end if
733        end do
734 !      Increase irmax a bit
735        irmax=irmax+4
736 !      Ask irmax to be lower than mmax
737        if(irmax>mmax-1)irmax=mmax-1
738 
739        ABI_MALLOC(work3,(irmax-1))
740        ABI_MALLOC(work4,(irmax-1))
741 
742 !      Loop over q values
743        do iq=2,mqgrid
744          tpiq=two_pi*qgrid(iq)
745          call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
746          work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility)
747          work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility)
748 
749 !        Handle r=0 separately
750          work1(1)=0.0d0
751          call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax))
752          do ir=2,irmax
753            work1(ir)=besjx(ir)*work5(ir)
754          end do
755 !        do integral
756          call der_int(work1,work2,rad,dr,irmax,result)
757          ffspl(iq,1,lp1)=result
758        end do
759 
760 !      Compute yp1=derivative of f(q) at q=0
761        if(lp1/=2)then
762 !        For l/=1, yp1=0
763          yp1=0.0d0
764        else
765 !        For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3
766          do ir=1,irmax
767            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
768            work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0
769          end do
770          call der_int(work1,work2,rad,dr,irmax,result)
771          yp1=result
772        end if
773 
774 !      Compute ypn=derivative of f(q) at q=qgrid(mqgrid)
775        tpiq=two_pi*qgrid(mqgrid)
776 !      Treat ir=1, r=0, separately
777        work1(1)=0.0d0
778 !      Here, must distinguish l==0 from others
779        if(lp1==1)then
780 !        l==0 : ypn=$\int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$
781 !        The sine and cosine of the last point were computed in the previous loop
782 !        So, there is no need to call sincos. Note that the rank of besj is 1.
783          call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax))
784          do ir=2,irmax
785            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
786            work1(ir)=-besjx(ir)*two_pi*rad(ir)*rad(ir)*dvwf
787          end do
788        else
789 !        l==1 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$
790 !        l==2 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$
791 !        l==3 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$
792 !        The sine and cosine of the last point were computed in the previous loop
793 !        Store first previously computed value with besj of order l, then use
794 !        besj of order l-1 (=lp1-2)
795          work1(2:irmax)=besjx(2:irmax)
796          call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax))
797          do ir=2,irmax
798            dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
799            work1(ir)=(two_pi*rad(ir)**2)*dvwf*&
800 &           ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) )
801          end do
802        end if
803 !      work1 is ready for integration
804        call der_int(work1,work2,rad,dr,irmax,result)
805        ypn=result
806 
807 !      Fit spline to get second derivatives by spline fit
808        call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,&
809 &       ffspl(:,2,lp1))
810 
811        ABI_FREE(work3)
812        ABI_FREE(work4)
813 
814      else
815 !      KB energy is zero, put nonlocal correction at l=0 to 0
816        ffspl(:,:,lp1)=0.0d0
817      end if
818 
819    end do !    End loop on angular momenta
820 
821    ABI_FREE(work1)
822    ABI_FREE(work2)
823    ABI_FREE(work_spl)
824    ABI_FREE(work5)
825    ABI_FREE(besjx)
826  end if !  End of lmax/=-1 condition
827 
828 end subroutine psp1nl

m_psp1/psp4cc [ Functions ]

[ Top ] [ m_psp1 ] [ Functions ]

NAME

 psp4cc

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for the format 4 of pseudopotentials.
 This is a even polynomial of 24th order for core density,
 that is cut off exactly beyond rchrg.
 It has been produced on 7 May 1992 by M. Teter.

INPUTS

  fchrg=magnitude of the core charge correction
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used

OUTPUT

  xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives

NOTES

 The argument of xccc1d is assumed to be normalized, and to vary
 from xx=0 to 1 (from r=0 to r=xcccrc)

WARNINGS

 the fifth derivative is not yet delivered.

SOURCE

1047 subroutine psp4cc(fchrg,n1xccc,xccc1d)
1048 
1049 !Arguments ------------------------------------
1050 !scalars
1051  integer,intent(in) :: n1xccc
1052  real(dp),intent(in) :: fchrg
1053 !arrays
1054  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
1055 
1056 !Local variables-------------------------------
1057 !scalars
1058  integer :: i1xccc,ider
1059  real(dp),parameter :: a10=-0.1156854803757563d5,a12=+0.2371534625455588d5
1060  real(dp),parameter :: a14=-0.3138755797827918d5,a16=+0.2582842713241039d5
1061  real(dp),parameter :: a18=-0.1200356429115204d5,a20=+0.2405099057118771d4
1062  real(dp),parameter :: a2=-0.8480751097855989d1,a4=+0.9684600878284791d2
1063  real(dp),parameter :: a6=-0.7490894651588015d3,a8=+0.3670890998130434d4
1064  real(dp) :: der1,dern,factor
1065  character(len=500) :: message
1066 !arrays
1067  real(dp),allocatable :: ff(:),ff2(:),work(:),xx(:)
1068  real(dp) :: x
1069 
1070 ! *************************************************************************
1071 
1072  ABI_MALLOC(ff,(n1xccc))
1073  ABI_MALLOC(ff2,(n1xccc))
1074  ABI_MALLOC(work,(n1xccc))
1075  ABI_MALLOC(xx,(n1xccc))
1076 
1077 
1078  if(n1xccc > 1)then
1079    factor=1.0d0/dble(n1xccc-1)
1080    do i1xccc=1,n1xccc
1081      xx(i1xccc)=(i1xccc-1)*factor
1082    end do
1083  else
1084    write(message, '(a,i0)' )'  n1xccc should larger than 1, while it is n1xccc=',n1xccc
1085    ABI_BUG(message)
1086  end if
1087 
1088 !Initialization, to avoid some problem with some compilers
1089  xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
1090 
1091 !Take care of each derivative separately
1092  do ider=0,2
1093 
1094    if(ider==0)then
1095 !    Generate spline fitting for the function gg
1096      do i1xccc=1,n1xccc
1097 !      ff(i1xccc)=fchrg*gg(xx(i1xccc))
1098        ff(i1xccc)=fchrg*gg_psp4(xx(i1xccc))
1099      end do
1100 !    Complete with derivatives at end points
1101      der1=0.0d0
1102 !    dern=fchrg*gp(1.0d0)
1103      dern=fchrg*gp_psp4(1.0d0)
1104    else if(ider==1)then
1105 !    Generate spline fitting for the function gp
1106      do i1xccc=1,n1xccc
1107 !      ff(i1xccc)=fchrg*gp(xx(i1xccc))
1108        ff(i1xccc)=fchrg*gp_psp4(xx(i1xccc))
1109      end do
1110 !    Complete with derivatives at end points, already estimated
1111      der1=xccc1d(1,ider+2)
1112      dern=xccc1d(n1xccc,ider+2)
1113    else if(ider==2)then
1114 !    Generate spline fitting for the function gpp
1115 !    (note : the function gpp has already been estimated, for the spline
1116 !    fitting of the function gg, but it is replaced here by the more
1117 !    accurate analytic derivative)
1118      do i1xccc=1,n1xccc
1119        x=xx(i1xccc)
1120        ff(i1xccc)=fchrg*(gpp_1_psp4(x)+gpp_2_psp4(x)+gpp_3_psp4(x))
1121 !      ff(i1xccc)=fchrg*gpp(xx(i1xccc))
1122      end do
1123 !    Complete with derivatives of end points
1124      der1=xccc1d(1,ider+2)
1125      dern=xccc1d(n1xccc,ider+2)
1126    end if
1127 
1128 !  Produce second derivative numerically, for use with splines
1129    call spline(xx,ff,n1xccc,der1,dern,ff2)
1130    xccc1d(:,ider+1)=ff(:)
1131    xccc1d(:,ider+3)=ff2(:)
1132  end do
1133 
1134  xccc1d(:,6)=zero
1135 
1136  ABI_FREE(ff)
1137  ABI_FREE(ff2)
1138  ABI_FREE(work)
1139  ABI_FREE(xx)
1140 
1141 !DEBUG
1142 !write(std_out,*)' psp1cc : output of core charge density and derivatives '
1143 !write(std_out,*)'   xx          gg           gp  '
1144 !do i1xccc=1,n1xccc
1145 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
1146 !end do
1147 !write(std_out,*)'   xx          gpp          gg2  '
1148 !do i1xccc=1,n1xccc
1149 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
1150 !end do
1151 !write(std_out,*)'   xx          gp2          gpp2  '
1152 !do i1xccc=1,n1xccc
1153 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
1154 !end do
1155 !write(std_out,*)' psp1cc : debug done, stop '
1156 !stop
1157 !ENDDEBUG
1158 
1159  contains
1160 
1161    function gg_psp4(x)
1162 !Expression of 7 May 1992
1163    real(dp) :: gg_psp4
1164    real(dp),intent(in) :: x
1165    gg_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + &
1166 &   x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ &
1167 &   x**2*(a18+x**2*(a20)))))))))))          *(1.0d0-x**2)**2
1168  end function gg_psp4
1169 
1170    function gp_psp4(x)
1171 !gp(x) is the derivative of gg(x) wrt x
1172    real(dp) :: gp_psp4
1173    real(dp),intent(in) :: x
1174    gp_psp4=2.d0*x*((a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*(              &
1175 &   4.d0*a8+x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*(                     &
1176 &   7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(10.d0*a20))))))))))*&
1177 &   (1.d0-x**2)**2                                                &
1178 &   -2.0d0*(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 +            &
1179 &   x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+            &
1180 &   x**2*(a18+x**2*a20))))))))))        *(1.0d0-x**2) )
1181  end function gp_psp4
1182 
1183    function gpp_1_psp4(x)
1184 !gpp(x) is the second derivative of gg(x) wrt x
1185    real(dp) :: gpp_1_psp4
1186    real(dp),intent(in) :: x
1187    gpp_1_psp4= ( 2.d0*a4+ x**2*(3.d0*2.d0*a6 +x**2*(               &
1188 &   4.d0*3.d0*a8+ x**2*(5.d0*4.d0*a10+x**2*(               &
1189 &   6.d0*5.d0*a12+x**2*(7.d0*6.d0*a14+x**2*(               &
1190 &   8.d0*7.d0*a16+x**2*(9.d0*8.d0*a18+x**2*(               &
1191 &   10.d0*9.d0*a20)                                        &
1192 &   ))))))))*(2.d0*x*(1.d0-x**2))**2
1193  end function gpp_1_psp4
1194 
1195    function gpp_2_psp4(x)
1196 
1197    real(dp) :: gpp_2_psp4
1198    real(dp),intent(in) :: x
1199    gpp_2_psp4=(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*(                 &
1200 &   4.d0*a8 +x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*(          &
1201 &   7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(          &
1202 &   10.d0*a20)                                             &
1203 &   )))))))))*(1.d0-x**2)*2*(1.d0-9.d0*x**2)
1204  end function gpp_2_psp4
1205 
1206    function gpp_3_psp4(x)
1207 
1208    real(dp) :: gpp_3_psp4
1209    real(dp),intent(in) :: x
1210    gpp_3_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 +         &
1211 &   x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+         &
1212 &   x**2*(a18+x**2*a20                               &
1213 &   ))))))))))*(1.0d0-3.d0*x**2)*(-4.d0)
1214  end function gpp_3_psp4
1215 
1216 end subroutine psp4cc

m_psp1/sincos [ Functions ]

[ Top ] [ m_psp1 ] [ Functions ]

NAME

 sincos

FUNCTION

 Update the sine and cosine values, needed inside the
 pseudopotential routines psp1lo and psp1nl.

INPUTS

  iq  = number of current wavevector q
  irmax = number of values  of r on the radial grid to be computed
  mmax = dimension of pspwk and rad
  pspwk(:,1,1) and pspwk(:,2,1) : sine and cosine of 2$\pi$ dq * rad
  pspwk(:,1,2) and pspwk(:,2,2) : sine and cosine of 2$\pi$ previous q * rad
  rad(mmax) radial grid
  tpiq = 2 $\pi$ * current wavevector q

OUTPUT

  pspwk(*,1,2) and pspwk(*,2,2) : sine and cosine of 2$\pi$ current q * rad

NOTES

 The speed was a special concern, so iterative computation
 based on addition formula is possible. Interestingly,
 this algorithm places strong constraints on accuracy,
 so this routine is machine-dependent.

SOURCE

 953 subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq)
 954 
 955 !Arguments ------------------------------------
 956 !scalars
 957  integer,intent(in) :: iq,irmax,mmax
 958  real(dp),intent(in) :: tpiq
 959 !arrays
 960  real(dp),intent(in) :: rad(mmax)
 961  real(dp),intent(inout) :: pspwk(mmax,2,2)
 962 
 963 !Local variables-------------------------------
 964 !scalars
 965  integer :: ir,nstep
 966  real(dp) :: prevcos,prevsin
 967  logical :: testmipspro
 968 
 969 
 970 ! *************************************************************************
 971 
 972  if(iq==2)then
 973 
 974 !  Here set up the sin and cos at iq=2
 975    do ir=2,irmax
 976      pspwk(ir,1,2)=pspwk(ir,1,1)
 977      pspwk(ir,2,2)=pspwk(ir,2,1)
 978    end do
 979 
 980  else
 981 !
 982 !  The sensitivity of the algorithm to changes of nstep
 983 !  has been tested : for all the machines except SGI - R10000 ,
 984 !  either using only the hard way, or
 985 !  using up to nstep=40 causes changes at the level
 986 !  of 1.0d-16 in the total energy. Larger values of
 987 !  nstep might be possible, but the associated residual
 988 !  is already very small ! The accelerated computation of
 989 !  sine and cosine is essential for a good speed on IBM, but,
 990 !  fortunately, on the SGI - R10000 the normal computation is fast enough.
 991 
 992    testmipspro=.false.
 993    nstep=40
 994    if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then
 995 
 996 !    Every nstep steps, uses the hard way
 997      do ir=2,irmax
 998        pspwk(ir,1,2)=sin(tpiq*rad(ir))
 999        pspwk(ir,2,2)=cos(tpiq*rad(ir))
1000      end do
1001 
1002    else
1003 
1004 !    Here the fastest way, iteratively
1005      do ir=2,irmax
1006        prevsin=pspwk(ir,1,2)
1007        prevcos=pspwk(ir,2,2)
1008        pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1)
1009        pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1)
1010      end do
1011 
1012    end if
1013 
1014  end if ! iq==2
1015 
1016 end subroutine sincos