TABLE OF CONTENTS


ABINIT/m_upf2abinit [ Modules ]

[ Top ] [ Modules ]

NAME

  m_upf2abinit

FUNCTION

COPYRIGHT

  Copyright (C) 2009-2018 ABINIT group (MJV)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

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_upf2abinit
28 
29  use defs_basis
30  use m_splines
31  use m_abicore
32  use m_errors
33  use m_atomdata
34  use pseudo_pwscf ! pwscf module with all data explicit!
35  use m_read_upf_pwscf, only : read_pseudo
36 
37  use defs_datatypes,  only : pseudopotential_type
38  use m_io_tools,      only : open_file
39  use m_numeric_tools, only : smooth, nderiv, ctrap
40  use m_pspheads,      only : upfxc2abi
41  use m_paw_numeric,   only : jbessel => paw_jbessel
42  use m_psptk,         only : cc_derivatives
43 
44  implicit none
45 
46  private

ABINIT/psp11lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp11lo

FUNCTION

 Compute sine transform to transform from V(r) to q^2 V(q).
 Computes integrals on logarithmic grid using related uniform
 grid in exponent and corrected trapezoidal integration.
 Generalized from psp5lo for non log grids using dr / di

INPUTS

  drdi=derivative of radial grid wrt index
  mmax=number of radial r grid points
  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=q grid values (bohr**-1).
  rad(mmax)=r grid values (bohr).
  vloc(mmax)=V(r) on radial grid.
  zion=nominal valence charge of atom.

OUTPUT

  epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
  q2vq(mqgrid)
   =q^2 V(q)
   = -\frac{Zv}{\pi}
     + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
  yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax
   (needed for spline fitter).

PARENTS

      upf2abinit

CHILDREN

      ctrap

SOURCE

474 subroutine psp11lo(drdi,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion)
475 
476 
477 !This section has been created automatically by the script Abilint (TD).
478 !Do not modify the following lines by hand.
479 #undef ABI_FUNC
480 #define ABI_FUNC 'psp11lo'
481 !End of the abilint section
482 
483  implicit none
484 
485 !Arguments----------------------------------------------------------
486 !scalars
487  integer,intent(in) :: mmax,mqgrid
488  real(dp),intent(in) :: zion
489  real(dp),intent(out) :: epsatm,yp1,ypn
490 !arrays
491  real(dp),intent(in) :: drdi(mmax)
492  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
493  real(dp),intent(out) :: q2vq(mqgrid)
494 
495 !Local variables-------------------------------
496 !scalars
497  integer :: iq,ir
498  real(dp),parameter :: scale=10.0d0
499  real(dp) :: arg,result_ctrap,test,ztor1
500 !arrays
501  real(dp),allocatable :: work(:)
502 
503 ! *************************************************************************
504 
505  ABI_ALLOCATE(work,(mmax))
506 
507 !Do q=0 separately (compute epsatm)
508 !Do integral from 0 to r1
509  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
510 
511 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
512 !with extra factor of drdi to convert to uniform grid
513  do ir = 1, mmax
514 !  First handle tail region
515    test=vloc(ir)+zion/rad(ir)
516 !  Ignore small contributions, or impose a cut-off in the case
517 !  the pseudopotential data are in single precision.
518 !  (it is indeed expected that vloc is very close to zero beyond 20,
519 !  so a value larger than 2.0d-8 is considered anomalous)
520    if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
521      work(ir)=zero
522    else
523      work(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
524    end if
525    work(ir)=work(ir)*drdi(ir)
526  end do
527 
528 !Do integral from r(1) to r(max)
529  call ctrap(mmax,work,one,result_ctrap)
530  epsatm=4.d0*pi*(result_ctrap+ztor1)
531 
532  q2vq(1)=-zion/pi
533 
534 !Loop over q values
535  do iq=2,mqgrid
536    arg=2.d0*pi*qgrid(iq)
537 !  ztor1=$ -Zv/\pi + 2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
538    ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion)* &
539 &   cos(arg*rad(1)) )/pi
540 
541 !  set up integrand
542    do  ir=1,mmax
543 !    test=vloc(ir)+zion/rad(ir)
544 !    Ignore contributions within decade of machine precision (suppressed ...)
545 !    if ((scale+abs(test)).eq.scale) then
546 !    work(ir)=zero
547 !    else
548      work(ir)=sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
549 !    end if
550      work(ir)=work(ir)*drdi(ir)
551    end do
552 !  do integral from r(1) to r(mmax)
553    call ctrap(mmax,work,one,result_ctrap)
554 
555 !  store q^2 v(q)
556 !  FIXME: I only see one factor q, not q^2, but the same is done in other pspXlo.F90
557    q2vq(iq)=ztor1+2.d0*qgrid(iq)*result_ctrap
558 
559  end do
560 
561 !Compute derivatives of q^2 v(q) at ends of interval
562  yp1=0.0d0
563 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
564 !integral from 0 to r1
565  arg=2.0d0*pi*qgrid(mqgrid)
566  ztor1=zion*rad(1)*sin(arg*rad(1))
567  ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + &
568 & (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1))
569 !integral from r(mmax) to infinity is overkill; ignore
570 !set up integrand
571  do ir=1,mmax
572 !  test=vloc(ir)+zion/rad(ir)
573 !  Ignore contributions within decade of machine precision (supressed ...)
574 !  if ((scale+abs(test)).eq.scale) then
575 !  work(ir)=0.0d0
576 !  else
577    work(ir)=(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * &
578 &   (rad(ir)*vloc(ir)+zion)
579 !  end if
580    work(ir)=work(ir)*drdi(ir)
581  end do
582  call ctrap(mmax,work,one,result_ctrap)
583  ypn=2.0d0 * (ztor1 + result_ctrap)
584 
585  ABI_DEALLOCATE(work)
586 
587 end subroutine psp11lo

ABINIT/upf2abinit [ Functions ]

[ Top ] [ Functions ]

NAME

 upf2abinit

FUNCTION

  This routine wraps a call to a PWSCF module, which reads in
  a UPF (PWSCF / Espresso) format pseudopotential, then transfers
  data to abinit internal variables.
  "UPF PWSCF format" (pspcod=11)

INPUTS

  filpsp = name of file with UPF data
  psps = sturcture with global dimension data for pseudopotentials, header info ...

OUTPUT

  pspxc = index of xc functional for this pseudo
  lmax_ = maximal angular momentum
  lloc = local component chosen for pseudopotential
  mmax = maximum number of points in real space radial grid
  znucl = charge of species nucleus
  zion = valence charge
  epsatm = integral of local potential - coulomb potential of zion
  xcccrc = radius for non linear core correction
  ekb(dimekb)= Kleinman Bylander energies, see pspatm.F90
  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
  indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln  (if useylm=0)
                                           or i=lmn (if useylm=1)
  vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
  ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum and
   each projector; if any, spin-orbit components begin at l=mpsang+1
  nproj= number of projectors for each channel
  xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives,
                              from psp file (used in NC only)

PARENTS

      pspatm

CHILDREN

      atomdata_from_symbol,cc_derivatives,nderiv,psp11lo,psp11nl,read_pseudo
      smooth,spline,upfxc2abi

SOURCE

 99 subroutine upf2abinit (filpsp, znucl, zion, pspxc, lmax_, lloc, mmax, &
100 &  psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj_l, vlspl, xccc1d)
101 
102 
103 !This section has been created automatically by the script Abilint (TD).
104 !Do not modify the following lines by hand.
105 #undef ABI_FUNC
106 #define ABI_FUNC 'upf2abinit'
107 !End of the abilint section
108 
109   implicit none
110 
111 !Arguments -------------------------------
112 
113   character(len=fnlen), intent(in) :: filpsp
114   type(pseudopotential_type),intent(in) :: psps
115   ! used contents:
116   !   psps%lmnmax
117   !   psps%mqgrid_ff
118   !   psps%mqgrid_vl
119   !   psps%dimekb
120   !   psps%n1xccc
121   !   psps%qgrid_ff
122   !   psps%qgrid_vl
123 
124   integer, intent(out) :: pspxc, lmax_, lloc, mmax
125   real(dp), intent(out) :: znucl, zion
126   real(dp), intent(out) :: epsatm, xcccrc
127   !arrays
128   integer, intent(out)  :: indlmn(6,psps%lmnmax)
129   integer, intent(out)  :: nproj_l(psps%mpssoang)
130   real(dp), intent(inout) :: ekb(psps%dimekb) !vz_i
131   real(dp), intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i
132   real(dp), intent(out) :: vlspl(psps%mqgrid_vl,2)
133   real(dp), intent(inout) :: xccc1d(psps%n1xccc,6) !vz_i
134 
135 !Local variables -------------------------
136   integer :: ir, iproj, ll, iunit
137   real(dp) :: yp1, ypn
138   character(len=500) :: msg
139   type(atomdata_t) :: atom
140 
141   logical, allocatable :: found_l(:)
142   real(dp), allocatable :: work_space(:),work_spl(:)
143   real(dp), allocatable :: ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:)
144 
145   ! ######### in module pseudo: ############
146   !
147   !  only npsx = 1 is used here
148   !  grids are allocated for much larger fixed length (ndm=2000)
149   !  number of species (6) and projectors (8) as well...
150   !
151   !  psd(npsx) = specied string
152   !  pseudotype = uspp / nc string
153   !  dft(npsx) = exchange correlation string (20 chars)
154   !  lmax(npsx) = maximum l channel
155   !  mesh(npsx) = number of points for local pot
156   !  nbeta(npsx) = number of projectors (beta functions for uspp)
157   !  nlcc(npsx) = flag for presence of NL core correction
158   !  zp(npsx) = valence ionic charge
159   !  r(ndm,npsx) = radial mesh
160   !  rab(ndm,npsx) = dr / di for radial mesh
161   !  rho_atc(ndm,npsx) = NLCC pseudocharge density
162   !  vloc0(ndm,npsx) = local pseudopotential
163   !  betar(ndm, nbrx, npsx) = projector functions in real space mesh
164   !  lll(nbrx,npsx) = angular momentum channel for each projector
165   !  ikk2(nbrx,npsx) = maximum index for each projector function
166   !  dion(nbrx,nbrx,npsx) = dij or Kleinman Bylander energies
167   !
168   !  ########  end description of pseudo module contents ##########
169 
170 ! *********************************************************************
171 
172 !call pwscf routine for reading in UPF
173  if (open_file (filpsp,msg,newunit=iunit,status='old',form='formatted') /= 0) then
174    MSG_ERROR(msg)
175  end if
176 
177 !read in psp data to static data in pseudo module, for ipsx == 1
178  call read_pseudo(1,iunit)
179  close (iunit)
180 
181 !convert to Ha units
182  vloc0  = half * vloc0
183 !betar = half * betar ! ???
184  dion   = half * dion
185 
186 !if upf file is a USPP one, stop
187  if (pseudotype == 'US') then
188    MSG_ERROR('upf2abinit: USPP UPF files not supported')
189  end if
190 
191 !copy over to abinit internal arrays and vars
192  call upfxc2abi(dft(1), pspxc)
193  lmax_ = lmax(1)
194 
195 !Check if the local component is one of the angular momentum channels
196 !effectively if one of the ll is absent from the NL projectors
197  ABI_ALLOCATE(found_l,(0:lmax_))
198  found_l = .true.
199  do ll = 0, lmax_
200    if (any(lll(1:nbeta(1),1) == ll)) then
201      found_l(ll) = .false.
202    end if
203  end do
204  if (count(found_l) /= 1) then
205    lloc = -1
206  else
207    do ll = 0, lmax_
208      if (found_l(ll)) then
209        lloc = ll
210        exit
211      end if
212    end do
213  end if
214  ABI_DEALLOCATE(found_l)
215 !FIXME: do something about lloc == -1
216 
217  call atomdata_from_symbol(atom,psd(1))
218  znucl = atom%znucl
219  zion = zp(1)
220  mmax = mesh(1)
221 
222  call psp11lo(rab(1:mmax,1),epsatm,mmax,psps%mqgrid_vl,psps%qgrid_vl,&
223 & vlspl(:,1),r(1:mmax,1),vloc0(1:mmax,1),yp1,ypn,zion)
224 
225 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
226  ABI_ALLOCATE(work_space,(psps%mqgrid_vl))
227  ABI_ALLOCATE(work_spl,(psps%mqgrid_vl))
228  call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl)
229  vlspl(:,2)=work_spl(:)
230  ABI_DEALLOCATE(work_space)
231  ABI_DEALLOCATE(work_spl)
232 
233 !this has to do the FT of the projectors to reciprocal space
234 ! allocate proj to avoid temporary copy.
235  ABI_ALLOCATE(proj, (mmax,1:nbeta(1)))
236  proj = betar(1:mmax,1:nbeta(1),1)
237 
238  call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, &
239 & nbeta(1), proj, lll(1:nbeta(1),1), ikk2(1:nbeta(1),1), &
240 & psps%qgrid_ff, r(1:mmax,1), rab(1:mmax,1), psps%useylm)
241 
242  ABI_FREE(proj)
243 
244  nproj_l = 0
245  do iproj = 1, nbeta(1)
246    ll = lll(iproj,1)
247    nproj_l(ll+1) = nproj_l(ll+1) + 1
248  end do
249 
250 !shape = dimekb  vs. shape = n_proj
251  do ll = 1, nbeta(1)
252    ekb(ll) = dion(ll,ll,1)
253  end do
254 
255  xcccrc = zero
256  xccc1d = zero
257 !if we find a core density, do something about it
258 !rho_atc contains the nlcc density
259 !rho_at contains the total density
260  if (nlcc(1)) then
261    ABI_ALLOCATE(ff,(mmax))
262    ABI_ALLOCATE(ff1,(mmax))
263    ABI_ALLOCATE(ff2,(mmax))
264    ff(1:mmax) = rho_atc(1:mmax,1) ! model core charge without derivative factor
265 
266    ff1 = zero
267    call nderiv(one,ff,ff1,mmax,1) ! first derivative
268    ff1(1:mmax) = ff1(1:mmax) / rab(1:mmax,1)
269    call smooth(ff1, mmax, 15) ! run 15 iterations of smoothing
270 
271    ff2 = zero
272    call nderiv(one, ff1, ff2, mmax, 1) ! second derivative
273    ff2(1:mmax) = ff2(1:mmax) / rab(1:mmax,1)
274    call smooth(ff2, mmax, 15) ! run 10 iterations of smoothing?
275 
276 !  determine a good rchrg = xcccrc
277    do ir = mmax, 1, -1
278      if (abs(ff(ir)) > 1.e-6) then
279        xcccrc = r(ir,1)
280        exit
281      end if
282    end do
283    ABI_ALLOCATE(rad_cc,(mmax))
284    rad_cc = r(1:mmax,1)
285    rad_cc(1) = zero ! force this so that the core charge covers whole spline interval.
286    call cc_derivatives(rad_cc,ff,ff1,ff2,mmax,psps%n1xccc,xcccrc,xccc1d)
287    ABI_DEALLOCATE(rad_cc)
288    ABI_DEALLOCATE(ff)
289    ABI_DEALLOCATE(ff1)
290    ABI_DEALLOCATE(ff2)
291 
292  end if !if nlcc present
293 
294 end subroutine upf2abinit

m_upf2abinit/psp11nl [ Functions ]

[ Top ] [ m_upf2abinit ] [ Functions ]

NAME

 psp11nl

FUNCTION

 Fourier transform the real space UPF projector functions to reciprocal space

INPUTS

  lmax=maximum ang momentum for which nonlocal form factor is desired.
   Usually lmax=1, sometimes = 0 (e.g. for oxygen); lmax <= 2 allowed.
  mmax=number of radial grid points for atomic grid
  lnmax= maximum index for all l channel projectors, dimension of ffspl
  lmnmax= maximum index for all projectors, dimension of indlmn
  mqgrid=number of grid points for q grid
  n_proj = total number of NL projectors read in
  proj = projector data times r, on a real space grid
  proj_l = ang mom channel for each projector
  proj_np = max number of points used for each projector
  qgrid(mqgrid)=values at which form factors are returned
  r(mmax)=radial grid values
  drdi=derivative of grid point wrt index
  useylm = input to use m dependency of NL part, or only Legendre polynomials

OUTPUT

  ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum
  indlmn = indexing of each projector, for n, l, m, s, ln, lmn (see pspatm.F90)

NOTES

PARENTS

      upf2abinit

CHILDREN

      ctrap,jbessel,spline

SOURCE

335 subroutine psp11nl(ffspl,indlmn,mmax,lnmax,lmnmax,mqgrid,n_proj,&
336 &                  proj, proj_l, proj_np, qgrid, r, drdi, useylm)
337 
338 
339 !This section has been created automatically by the script Abilint (TD).
340 !Do not modify the following lines by hand.
341 #undef ABI_FUNC
342 #define ABI_FUNC 'psp11nl'
343 !End of the abilint section
344 
345  implicit none
346 
347 !Arguments ------------------------------------
348 !scalars
349  integer,intent(in) :: mmax, lnmax, lmnmax, mqgrid, useylm, n_proj
350 !arrays
351  integer, intent(in) :: proj_l(n_proj)
352  integer, intent(in) :: proj_np(n_proj)
353  integer, intent(out) :: indlmn(6,lmnmax)
354  real(dp),intent(in) :: r(mmax)
355  real(dp),intent(in) :: drdi(mmax)
356  real(dp),intent(in) :: proj(mmax,n_proj)
357  real(dp),intent(in) :: qgrid(mqgrid)
358  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
359 
360 !Local variables-------------------------------
361 !scalars
362  integer :: iproj, np, ll, llold, ipsang, i_indlmn
363  integer :: iproj_1l, ir, iq, mm
364  integer :: bessorder
365  real(dp) :: res, arg, besfact, dummy, dummy2
366  real(dp), allocatable :: work(:)
367  character(len=500) :: message
368 
369 !*************************************************************************
370  bessorder = 0 ! never calculate derivatives of bessel functions
371 
372  ffspl = zero
373  indlmn = 0
374  i_indlmn = 0
375  llold = -1
376  iproj_1l = 1
377 !big loop over all projectors
378  do iproj = 1, n_proj
379 
380    if (iproj > lmnmax) then
381      write(message,'(a,2i0)') ' Too many projectors found. n_proj, lmnmax =  ',n_proj, lmnmax
382      MSG_ERROR(message)
383    end if
384 
385    np = proj_np(iproj)
386    ABI_ALLOCATE(work,(np))
387    ll = proj_l(iproj)
388    if (ll < llold) then
389      message = 'psp11nl : Error: UPF projectors are not in order of increasing ll'
390      MSG_ERROR(message)
391    else if (ll == llold) then
392      iproj_1l = iproj_1l + 1
393    else
394      iproj_1l = 1
395      llold = ll
396    end if
397 !  determine indlmn for this projector (keep in UPF order and enforce that they are in
398 !  increasing ll)
399    do mm = 1, 2*ll*useylm+1
400      i_indlmn = i_indlmn + 1
401      indlmn(1,i_indlmn) = ll
402      indlmn(2,i_indlmn) = mm-ll*useylm-1
403      indlmn(3,i_indlmn) = iproj_1l
404      indlmn(4,i_indlmn) = ll*ll+(1-useylm)*ll+mm
405      indlmn(5,i_indlmn) = iproj
406      indlmn(6,i_indlmn) = 1 !spin? FIXME: to get j for relativistic cases
407    end do
408 
409 !  FT projectors to reciprocal space q
410    do iq = 1, mqgrid
411      arg = two_pi*qgrid(iq)
412 
413 !    FIXME: add semianalytic form for integral from 0 to first point
414      do ir = 1, np
415        call jbessel(besfact, dummy, dummy2, ll, bessorder, arg*r(ir))
416 !      besfact = sin(arg*r(ir))
417        work(ir) = drdi(ir) * besfact * proj(ir, iproj) * r(ir) !* r(ir)
418      end do
419      call ctrap (np, work, one, res)
420 
421      ffspl(iq, 1, iproj) = res
422    end do
423    ABI_DEALLOCATE(work)
424  end do  ! iproj
425 
426 !add derivative of ffspl(:,1,:) for spline interpolation later
427  ABI_ALLOCATE(work,(mqgrid))
428  do ipsang = 1, lnmax
429    call spline(qgrid,ffspl(:,1,ipsang),mqgrid,zero,zero,ffspl(:,2,ipsang))
430  end do
431  ABI_DEALLOCATE(work)
432 
433 end subroutine psp11nl