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

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_psp1
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_splines
33 
34  use m_special_funcs,   only : besjm
35  use m_psptk,           only : psp1cc
36 
37  implicit none
38 
39  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 $.

PARENTS

      psp1lo,psp1nl

CHILDREN

SOURCE

908 subroutine der_int(ff,df,rr,dr,nlast,smf)
909 
910 
911 !This section has been created automatically by the script Abilint (TD).
912 !Do not modify the following lines by hand.
913 #undef ABI_FUNC
914 #define ABI_FUNC 'der_int'
915 !End of the abilint section
916 
917  implicit none
918 
919 !Arguments ------------------------------------
920 !nmax sets standard number of grid points ! SHOULD BE REMOVED
921 !scalars
922  integer,parameter :: nmax=2000
923  integer,intent(in) :: nlast
924  real(dp),intent(out) :: smf
925 !no_abirules
926 !Note that dimension here starts at 0
927  real(dp), intent(in) :: dr(0:nmax),ff(0:nmax),rr(0:nmax)
928  real(dp), intent(out) :: df(0:nmax)
929 
930 !Local variables-------------------------------
931 !scalars
932  integer :: jj
933  real(dp),parameter :: div12=1.d0/12.d0
934  real(dp) :: hh
935  character(len=500) :: message
936 
937 ! *************************************************************************
938 
939 !Check that nlast lie within 0 to nmax
940  if (nlast<0.or.nlast>nmax) then
941    write(message, '(a,i12,a,i12)' )&
942 &   ' nlast=',nlast,' lies outside range [0,nmax] with dimension nmax=',nmax
943    MSG_BUG(message)
944  end if
945 
946 !Compute derivatives at lower end, near r=0
947  df(0)=-25.d0/12.d0*ff(0)+4.d0*ff(1)-3.d0*ff(2)+4.d0/3.d0*ff(3)&
948 & -1.d0/4.d0*ff(4)
949  df(1)=-1.d0/4.d0*ff(0)-5.d0/6.d0*ff(1)+3.d0/2.d0*ff(2)&
950 & -1.d0/2.d0*ff(3)+1.d0/12.d0*ff(4)
951 
952 !Run over range from just past r=0 to near r(n), using central differences
953  do jj=2,nlast-2
954    df(jj)=(ff(jj-2)-8.d0*(ff(jj-1)-ff(jj+1))-ff(jj+2))*div12
955  end do
956 
957 !Compute derivative at upper end of range
958  if (nlast < 4) then
959    message = ' der_int: ff does not have enough elements. nlast is too low'
960    MSG_ERROR(message)
961  end if
962 
963  df(nlast-1)=-1.d0/12.d0*ff(nlast-4)&
964 & +1.d0/2.d0*ff(nlast-3)&
965 & -3.d0/2.d0*ff(nlast-2)&
966 & +5.d0/6.d0*ff(nlast-1)&
967 & +1.d0/4.d0*ff(nlast)
968  df(nlast)=1.d0/4.d0*ff(nlast-4)&
969 & -4.d0/3.d0*ff(nlast-3)&
970 & +3.d0*ff(nlast-2)&
971 & -4.d0*ff(nlast-1)&
972 & +25.d0/12.d0*ff(nlast)
973 
974 !Apply correct normalization over full range
975  do jj=0,nlast
976    df(jj)=df(jj)*dr(jj)
977  end do
978 
979  smf=0.d0
980  do jj=0,nlast-1
981    hh=rr(jj+1)-rr(jj)
982    smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1)))
983  end do
984  smf=smf/12.d0
985 
986 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 .

PARENTS

      pspatm

CHILDREN

      psp1cc,psp1lo,psp1nl,psp4cc,spline,wrtout

SOURCE

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

PARENTS

      psp1in

CHILDREN

      der_int,sincos

SOURCE

457 subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
458 &  vloc,wksincos,yp1,ypn,zion)
459 
460 
461 !This section has been created automatically by the script Abilint (TD).
462 !Do not modify the following lines by hand.
463 #undef ABI_FUNC
464 #define ABI_FUNC 'psp1lo'
465 !End of the abilint section
466 
467  implicit none
468 
469 !Arguments ------------------------------------
470 !scalars
471  integer,intent(in) :: mmax,mqgrid
472  real(dp),intent(in) :: zion
473  real(dp),intent(out) :: epsatm,yp1,ypn
474 !arrays
475  real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax)
476  real(dp),intent(inout) :: wksincos(mmax,2,2)
477  real(dp),intent(out) :: q2vq(mqgrid)
478 
479 !Local variables-------------------------------
480 !scalars
481  integer,parameter :: mma0=2001
482  integer :: iq,ir,irmax
483  real(dp),parameter :: scale=10.0d0
484  real(dp) :: result,test,tpiq
485 !arrays
486  real(dp) :: wk(mma0),wk1(mma0),wk2(mma0)
487 
488 ! *************************************************************************
489 
490 !Do q=0 separately (compute epsatm)
491 !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr]
492 !Treat r=0 by itself
493  wk(1)=0.0d0
494 
495  do ir=2,mmax
496 !  (at large r do not want prefactor of r^2 and should see
497 !  V(r)+Zv/r go to 0 at large r)
498    test=vloc(ir)+zion/rad(ir)
499 !  write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test
500 !  In this routine, NO cut-off radius is imposed : the input
501 !  vloc MUST be in real(dp) to obtain numerically
502 !  accurate values. The error can be on the order of 0.001 Ha !
503    if (abs(test)<1.0d-20) then
504      wk(ir)=0.0d0
505    else
506      wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
507    end if
508  end do
509 !Do integral from 0 to r(max) (disregard contrib beyond r(max)
510 !(need numerical derivatives to do integral)
511 !Use mmax-1 to convert to Teter s dimensioning starting at 0
512  call der_int(wk,wk2,rad,drad,mmax-1,result)
513 
514  epsatm=4.d0*pi*(result)
515 !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi
516  q2vq(1)=-zion/pi
517 
518 !Prepare loop over q values
519  irmax=mmax+1
520  do ir=mmax,2,-1
521    test=vloc(ir)+zion/rad(ir)
522    wk1(ir)=test*rad(ir)
523 !  Will ignore tail within decade of machine precision
524    if ((scale+abs(test))==scale .and. irmax==ir+1) then
525      irmax=ir
526    end if
527  end do
528 !Increase irmax a bit : this is copied from psp1nl
529  irmax=irmax+4
530  if(irmax>mmax)irmax=mmax
531 
532 !Loop over q values
533  do iq=2,mqgrid
534    tpiq=two_pi*qgrid(iq)
535    call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
536 !  set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral
537 !$\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]$.
538 !  Handle r=0 separately
539    wk(1)=0.0d0
540    do ir=2,irmax
541      wk(ir)=wksincos(ir,1,2)*wk1(ir)
542    end do
543 !  do integral from 0 to r(max)
544    if(irmax>mmax-1)irmax=mmax-1
545 
546    call der_int(wk,wk2,rad,drad,irmax,result)
547 !  store q^2 v(q)
548    q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result
549  end do
550 
551 !Compute derivatives of q^2 v(q) at ends of interval
552  yp1=0.0d0
553 !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]$
554 !integral from r(mmax) to infinity is overkill; ignore
555 !set up integrand
556 !Handle r=0 separately
557  wk(1)=0.0d0
558  tpiq=two_pi*qgrid(mqgrid)
559  do ir=2,mmax
560    test=vloc(ir)+zion/rad(ir)
561 !  Ignore contributions within decade of machine precision
562    if ((scale+abs(test))==scale) then
563      wk(ir)=0.0d0
564    else
565      wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * &
566 &     (rad(ir)*vloc(ir)+zion)
567    end if
568  end do
569  call der_int(wk,wk2,rad,drad,mmax-1,result)
570 
571  ypn=2.0d0*result
572 
573 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.

PARENTS

      psp1in

CHILDREN

      besjm,der_int,sincos,spline

SOURCE

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

PARENTS

      psp1in

CHILDREN

      spline

SOURCE

1146 subroutine psp4cc(fchrg,n1xccc,xccc1d)
1147 
1148 
1149 !This section has been created automatically by the script Abilint (TD).
1150 !Do not modify the following lines by hand.
1151 #undef ABI_FUNC
1152 #define ABI_FUNC 'psp4cc'
1153 !End of the abilint section
1154 
1155  implicit none
1156 
1157 !Arguments ------------------------------------
1158 !scalars
1159  integer,intent(in) :: n1xccc
1160  real(dp),intent(in) :: fchrg
1161 !arrays
1162  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
1163 
1164 !Local variables-------------------------------
1165 !scalars
1166  integer :: i1xccc,ider
1167  real(dp),parameter :: a10=-0.1156854803757563d5,a12=+0.2371534625455588d5
1168  real(dp),parameter :: a14=-0.3138755797827918d5,a16=+0.2582842713241039d5
1169  real(dp),parameter :: a18=-0.1200356429115204d5,a20=+0.2405099057118771d4
1170  real(dp),parameter :: a2=-0.8480751097855989d1,a4=+0.9684600878284791d2
1171  real(dp),parameter :: a6=-0.7490894651588015d3,a8=+0.3670890998130434d4
1172  real(dp) :: der1,dern,factor
1173  character(len=500) :: message
1174 !arrays
1175  real(dp),allocatable :: ff(:),ff2(:),work(:),xx(:)
1176  real(dp) :: x
1177 
1178 ! *************************************************************************
1179 
1180  ABI_ALLOCATE(ff,(n1xccc))
1181  ABI_ALLOCATE(ff2,(n1xccc))
1182  ABI_ALLOCATE(work,(n1xccc))
1183  ABI_ALLOCATE(xx,(n1xccc))
1184 
1185 
1186  if(n1xccc > 1)then
1187    factor=1.0d0/dble(n1xccc-1)
1188    do i1xccc=1,n1xccc
1189      xx(i1xccc)=(i1xccc-1)*factor
1190    end do
1191  else
1192    write(message, '(a,i0)' )'  n1xccc should larger than 1, while it is n1xccc=',n1xccc
1193    MSG_BUG(message)
1194  end if
1195 
1196 !Initialization, to avoid some problem with some compilers
1197  xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
1198 
1199 !Take care of each derivative separately
1200  do ider=0,2
1201 
1202    if(ider==0)then
1203 !    Generate spline fitting for the function gg
1204      do i1xccc=1,n1xccc
1205 !      ff(i1xccc)=fchrg*gg(xx(i1xccc))
1206        ff(i1xccc)=fchrg*gg_psp4(xx(i1xccc))
1207      end do
1208 !    Complete with derivatives at end points
1209      der1=0.0d0
1210 !    dern=fchrg*gp(1.0d0)
1211      dern=fchrg*gp_psp4(1.0d0)
1212    else if(ider==1)then
1213 !    Generate spline fitting for the function gp
1214      do i1xccc=1,n1xccc
1215 !      ff(i1xccc)=fchrg*gp(xx(i1xccc))
1216        ff(i1xccc)=fchrg*gp_psp4(xx(i1xccc))
1217      end do
1218 !    Complete with derivatives at end points, already estimated
1219      der1=xccc1d(1,ider+2)
1220      dern=xccc1d(n1xccc,ider+2)
1221    else if(ider==2)then
1222 !    Generate spline fitting for the function gpp
1223 !    (note : the function gpp has already been estimated, for the spline
1224 !    fitting of the function gg, but it is replaced here by the more
1225 !    accurate analytic derivative)
1226      do i1xccc=1,n1xccc
1227        x=xx(i1xccc)
1228        ff(i1xccc)=fchrg*(gpp_1_psp4(x)+gpp_2_psp4(x)+gpp_3_psp4(x))
1229 !      ff(i1xccc)=fchrg*gpp(xx(i1xccc))
1230      end do
1231 !    Complete with derivatives of end points
1232      der1=xccc1d(1,ider+2)
1233      dern=xccc1d(n1xccc,ider+2)
1234    end if
1235 
1236 !  Produce second derivative numerically, for use with splines
1237    call spline(xx,ff,n1xccc,der1,dern,ff2)
1238    xccc1d(:,ider+1)=ff(:)
1239    xccc1d(:,ider+3)=ff2(:)
1240  end do
1241 
1242  xccc1d(:,6)=zero
1243 
1244  ABI_DEALLOCATE(ff)
1245  ABI_DEALLOCATE(ff2)
1246  ABI_DEALLOCATE(work)
1247  ABI_DEALLOCATE(xx)
1248 
1249 !DEBUG
1250 !write(std_out,*)' psp1cc : output of core charge density and derivatives '
1251 !write(std_out,*)'   xx          gg           gp  '
1252 !do i1xccc=1,n1xccc
1253 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
1254 !end do
1255 !write(std_out,*)'   xx          gpp          gg2  '
1256 !do i1xccc=1,n1xccc
1257 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
1258 !end do
1259 !write(std_out,*)'   xx          gp2          gpp2  '
1260 !do i1xccc=1,n1xccc
1261 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
1262 !end do
1263 !write(std_out,*)' psp1cc : debug done, stop '
1264 !stop
1265 !ENDDEBUG
1266 
1267  contains
1268 
1269    function gg_psp4(x)
1270 !Expression of 7 May 1992
1271 
1272 !This section has been created automatically by the script Abilint (TD).
1273 !Do not modify the following lines by hand.
1274 #undef ABI_FUNC
1275 #define ABI_FUNC 'gg_psp4'
1276 !End of the abilint section
1277 
1278    real(dp) :: gg_psp4
1279    real(dp),intent(in) :: x
1280    gg_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + &
1281 &   x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ &
1282 &   x**2*(a18+x**2*(a20)))))))))))          *(1.0d0-x**2)**2
1283  end function gg_psp4
1284 
1285    function gp_psp4(x)
1286 !gp(x) is the derivative of gg(x) wrt x
1287 
1288 !This section has been created automatically by the script Abilint (TD).
1289 !Do not modify the following lines by hand.
1290 #undef ABI_FUNC
1291 #define ABI_FUNC 'gp_psp4'
1292 !End of the abilint section
1293 
1294    real(dp) :: gp_psp4
1295    real(dp),intent(in) :: x
1296    gp_psp4=2.d0*x*((a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*(              &
1297 &   4.d0*a8+x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*(                     &
1298 &   7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(10.d0*a20))))))))))*&
1299 &   (1.d0-x**2)**2                                                &
1300 &   -2.0d0*(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 +            &
1301 &   x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+            &
1302 &   x**2*(a18+x**2*a20))))))))))        *(1.0d0-x**2) )
1303  end function gp_psp4
1304 
1305    function gpp_1_psp4(x)
1306 !gpp(x) is the second derivative of gg(x) wrt x
1307 
1308 !This section has been created automatically by the script Abilint (TD).
1309 !Do not modify the following lines by hand.
1310 #undef ABI_FUNC
1311 #define ABI_FUNC 'gpp_1_psp4'
1312 !End of the abilint section
1313 
1314    real(dp) :: gpp_1_psp4
1315    real(dp),intent(in) :: x
1316    gpp_1_psp4= ( 2.d0*a4+ x**2*(3.d0*2.d0*a6 +x**2*(               &
1317 &   4.d0*3.d0*a8+ x**2*(5.d0*4.d0*a10+x**2*(               &
1318 &   6.d0*5.d0*a12+x**2*(7.d0*6.d0*a14+x**2*(               &
1319 &   8.d0*7.d0*a16+x**2*(9.d0*8.d0*a18+x**2*(               &
1320 &   10.d0*9.d0*a20)                                        &
1321 &   ))))))))*(2.d0*x*(1.d0-x**2))**2
1322  end function gpp_1_psp4
1323 
1324    function gpp_2_psp4(x)
1325 
1326 
1327 !This section has been created automatically by the script Abilint (TD).
1328 !Do not modify the following lines by hand.
1329 #undef ABI_FUNC
1330 #define ABI_FUNC 'gpp_2_psp4'
1331 !End of the abilint section
1332 
1333    real(dp) :: gpp_2_psp4
1334    real(dp),intent(in) :: x
1335    gpp_2_psp4=(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*(                 &
1336 &   4.d0*a8 +x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*(          &
1337 &   7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(          &
1338 &   10.d0*a20)                                             &
1339 &   )))))))))*(1.d0-x**2)*2*(1.d0-9.d0*x**2)
1340  end function gpp_2_psp4
1341 
1342    function gpp_3_psp4(x)
1343 
1344 
1345 !This section has been created automatically by the script Abilint (TD).
1346 !Do not modify the following lines by hand.
1347 #undef ABI_FUNC
1348 #define ABI_FUNC 'gpp_3_psp4'
1349 !End of the abilint section
1350 
1351    real(dp) :: gpp_3_psp4
1352    real(dp),intent(in) :: x
1353    gpp_3_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 +         &
1354 &   x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+         &
1355 &   x**2*(a18+x**2*a20                               &
1356 &   ))))))))))*(1.0d0-3.d0*x**2)*(-4.d0)
1357  end function gpp_3_psp4
1358 
1359 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.

PARENTS

      psp1lo,psp1nl

CHILDREN

SOURCE

1022 subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq)
1023 
1024 
1025 !This section has been created automatically by the script Abilint (TD).
1026 !Do not modify the following lines by hand.
1027 #undef ABI_FUNC
1028 #define ABI_FUNC 'sincos'
1029 !End of the abilint section
1030 
1031  implicit none
1032 
1033 !Arguments ------------------------------------
1034 !scalars
1035  integer,intent(in) :: iq,irmax,mmax
1036  real(dp),intent(in) :: tpiq
1037 !arrays
1038  real(dp),intent(in) :: rad(mmax)
1039  real(dp),intent(inout) :: pspwk(mmax,2,2)
1040 
1041 !Local variables-------------------------------
1042 !scalars
1043  integer :: ir,nstep
1044  real(dp) :: prevcos,prevsin
1045  logical :: testmipspro
1046 #if defined HAVE_LINALG_MLIB
1047  real(dp) :: halfpi
1048 #endif
1049 
1050 
1051 ! *************************************************************************
1052 
1053 #if defined HAVE_LINALG_MLIB
1054  halfpi=asin(1.0d0)
1055 #endif
1056 
1057  if(iq==2)then
1058 
1059 !  Here set up the sin and cos at iq=2
1060    do ir=2,irmax
1061      pspwk(ir,1,2)=pspwk(ir,1,1)
1062      pspwk(ir,2,2)=pspwk(ir,2,1)
1063    end do
1064 
1065  else
1066 !
1067 !  The sensitivity of the algorithm to changes of nstep
1068 !  has been tested : for all the machines except SGI - R10000 ,
1069 !  either using only the hard way, or
1070 !  using up to nstep=40 causes changes at the level
1071 !  of 1.0d-16 in the total energy. Larger values of
1072 !  nstep might be possible, but the associated residual
1073 !  is already very small ! The accelerated computation of
1074 !  sine and cosine is essential for a good speed on IBM, but,
1075 !  fortunately, on the SGI - R10000 the normal computation is fast enough.
1076 
1077    testmipspro=.false.
1078 #ifdef FC_MIPSPRO
1079    testmipspro=.true.
1080 #endif
1081    nstep=40
1082    if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then
1083 
1084 !    Every nstep steps, uses the hard way
1085      do ir=2,irmax
1086 #if defined HAVE_LINALG_MLIB
1087 !      There is a bug in the hp library !! Sine is slightly inaccurate !
1088        pspwk(ir,1,2)=cos(tpiq*rad(ir)-halfpi)
1089 #else
1090        pspwk(ir,1,2)=sin(tpiq*rad(ir))
1091 #endif
1092        pspwk(ir,2,2)=cos(tpiq*rad(ir))
1093      end do
1094 
1095    else
1096 
1097 !    Here the fastest way, iteratively
1098      do ir=2,irmax
1099        prevsin=pspwk(ir,1,2)
1100        prevcos=pspwk(ir,2,2)
1101        pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1)
1102        pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1)
1103      end do
1104 
1105    end if
1106 
1107  end if ! iq==2
1108 
1109 end subroutine sincos