TABLE OF CONTENTS


ABINIT/m_upf2abinit [ Modules ]

[ Top ] [ Modules ]

NAME

  m_upf2abinit

FUNCTION

  Procedures to read NC pseudos in UPF1/UPF2 format and convert data into
  the internal ABINIT representation.

COPYRIGHT

  Copyright (C) 2009-2022 ABINIT group (MJV, MG, DRH)
  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

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_upf2abinit
24 
25  use defs_basis
26  use m_abicore
27  use m_errors
28  use m_atomdata
29  use m_splines
30 
31  use defs_datatypes,  only : pseudopotential_type
32  use m_io_tools,      only : open_file
33  use m_numeric_tools, only : smooth, nderiv, ctrap
34  use m_copy,          only : alloc_copy
35  use m_paw_numeric,   only : jbessel => paw_jbessel
36  use m_pawpsp,        only : pawpsp_nl
37  use m_pawrad,        only : pawrad_type, pawrad_init, pawrad_free, simp_gen
38  use m_psptk,         only : cc_derivatives, psp8lo, psp8nl
39 
40  implicit none
41 
42  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 weights.

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

SOURCE

861 subroutine psp11lo(drdi,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion)
862 
863 !Arguments----------------------------------------------------------
864 !scalars
865  integer,intent(in) :: mmax,mqgrid
866  real(dp),intent(in) :: zion
867  real(dp),intent(out) :: epsatm,yp1,ypn
868 !arrays
869  real(dp),intent(in) :: drdi(mmax)
870  real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
871  real(dp),intent(out) :: q2vq(mqgrid)
872 
873 !Local variables-------------------------------
874 !scalars
875  integer :: iq,ir
876  real(dp),parameter :: scale=10.0d0
877  real(dp) :: arg,result_ctrap,test,ztor1
878 !arrays
879  real(dp),allocatable :: work(:)
880 
881 ! *************************************************************************
882 
883  ABI_MALLOC(work,(mmax))
884 
885  ! Do q=0 separately (compute epsatm)
886  ! Do integral from 0 to r1
887  ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
888 
889  ! Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
890  ! with extra factor of drdi to convert to uniform grid
891  do ir = 1, mmax
892    ! First handle tail region
893    test=vloc(ir)+zion/rad(ir)
894    ! Ignore small contributions, or impose a cut-off in the case
895    ! the pseudopotential data are in single precision.
896    ! (it is indeed expected that vloc is very close to zero beyond 20,
897    ! so a value larger than 2.0d-8 is considered anomalous)
898    if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
899      work(ir)=zero
900    else
901      work(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
902    end if
903    work(ir)=work(ir)*drdi(ir)
904  end do
905 
906  ! Do integral from r(1) to r(max)
907  call ctrap(mmax,work,one,result_ctrap)
908  epsatm=4.d0*pi*(result_ctrap+ztor1)
909 
910  q2vq(1)=-zion/pi
911 
912  ! Loop over q values
913  do iq=2,mqgrid
914    arg=2.d0*pi*qgrid(iq)
915    ! ztor1=$ -Zv/\pi + 2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
916    ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion) * cos(arg*rad(1)) )/pi
917 
918    ! set up integrand
919    do  ir=1,mmax
920      !test=vloc(ir)+zion/rad(ir)
921      !Ignore contributions within decade of machine precision (suppressed ...)
922      !if ((scale+abs(test)).eq.scale) then
923      !work(ir)=zero
924      !else
925      work(ir)=sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
926      !end if
927      work(ir)=work(ir)*drdi(ir)
928    end do
929    ! do integral from r(1) to r(mmax)
930    call ctrap(mmax,work,one,result_ctrap)
931 
932    ! store q^2 v(q)
933    ! FIXME: I only see one factor q, not q^2, but the same is done in other pspXlo.F90
934    q2vq(iq)=ztor1+2.d0*qgrid(iq)*result_ctrap
935 
936  end do
937 
938  ! Compute derivatives of q^2 v(q) at ends of interval
939  yp1=0.0d0
940  !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
941  !integral from 0 to r1
942  arg=2.0d0*pi*qgrid(mqgrid)
943  ztor1=zion*rad(1)*sin(arg*rad(1))
944  ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1))
945  !integral from r(mmax) to infinity is overkill; ignore
946  !set up integrand
947  do ir=1,mmax
948    !test=vloc(ir)+zion/rad(ir)
949    !Ignore contributions within decade of machine precision (supressed ...)
950    !if ((scale+abs(test)).eq.scale) then
951    !work(ir)=0.0d0
952    !else
953    work(ir)=(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * (rad(ir)*vloc(ir)+zion)
954    !end if
955    work(ir)=work(ir)*drdi(ir)
956  end do
957 
958  call ctrap(mmax,work,one,result_ctrap)
959  ypn=2.0d0 * (ztor1 + result_ctrap)
960  ABI_FREE(work)
961 
962 end subroutine psp11lo

ABINIT/upf1_to_abinit [ Functions ]

[ Top ] [ Functions ]

NAME

 upf1_to_abinit

FUNCTION

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

INPUTS

  filpsp = name of file with UPF data
  psps = sturcture with global dimension data for pseudopotentials, header info ...
    used contents:
       psps%lmnmax
       psps%mqgrid_ff
       psps%mqgrid_vl
       psps%dimekb
       psps%n1xccc
       psps%qgrid_ff
       psps%qgrid_vl

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)

SOURCE

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

ABINIT/upf2_to_abinit [ Functions ]

[ Top ] [ Functions ]

NAME

 upf2_to_abinit

FUNCTION

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

INPUTS

  filpsp = name of file with UPF2 data
  psps = sturcture with global dimension data for pseudopotentials, header info ...
    used contents:
       psps%lmnmax
       psps%mqgrid_ff
       psps%mqgrid_vl
       psps%dimekb
       psps%n1xccc
       psps%qgrid_ff
       psps%qgrid_vl

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)
  nctab<nctab_t>=NC tables
    %has_tvale=True if the pseudo contains the pseudo valence charge
    %tvalespl(mqgrid_vl,2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid

SOURCE

331 subroutine upf2_to_abinit(ipsp, filpsp, znucl, zion, pspxc, lmax, lloc, mmax, &
332                           psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj_l, vlspl, xccc1d, nctab, maxrad)
333 
334  use pseudo_types,        only : pseudo_upf, deallocate_pseudo_upf !, pseudo_config
335  use read_upf_new_module, only : read_upf_new
336  use defs_datatypes,      only : nctab_t
337  use m_psps,              only : nctab_eval_tvalespl
338  use m_pspheads,          only : upfdft_to_ixc, upf2_jl2srso
339 
340 !Arguments -------------------------------
341  integer,intent(in) :: ipsp
342  character(len=fnlen), intent(in) :: filpsp
343  type(pseudopotential_type),intent(in) :: psps
344  type(nctab_t),intent(inout) :: nctab
345  integer, intent(out) :: pspxc, lmax, lloc, mmax
346  real(dp), intent(out) :: znucl, zion
347  real(dp), intent(out) :: epsatm, xcccrc, maxrad
348  !arrays
349  integer, intent(out)  :: indlmn(6,psps%lmnmax), nproj_l(psps%mpssoang)
350  real(dp), intent(inout) :: ekb(psps%dimekb)
351  real(dp), intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax)
352  real(dp), intent(out) :: vlspl(psps%mqgrid_vl,2)
353  real(dp), intent(inout) :: xccc1d(psps%n1xccc,6)
354 
355 !Local variables -------------------------
356  integer :: ierr, ir, irad, iprj, il, ll, smooth_niter, nso, nn, iln, kk, mm, pspindex, iwfc !, iq
357  integer :: atmwfc_lmax
358  real(dp) :: yp1, ypn, amesh, damesh, intg
359  character(len=500) :: msg
360  logical :: linear_mesh
361  type(pseudo_upf) :: upf
362  type(atomdata_t) :: atom
363  type(pawrad_type) :: mesh
364  integer :: my_nproj_l(0:3), my_nprojso_l(1:3)
365  !integer :: nproj_tmp(psps%mpssoang)
366  integer,allocatable :: awfc_indlmn(:,:)
367  logical,allocatable :: found_l(:)
368  real(dp),allocatable :: work_space(:), work_spl(:)
369  real(dp),allocatable :: ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:)
370  real(dp),allocatable :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:)
371 
372 ! *********************************************************************
373 
374  ! See also https://github.com/QEF/qeschemas/blob/master/UPF/qe_pp-0.99.xsd
375  ! and https://github.com/QEF/qeschemas/files/9497267/pp.md
376  call read_upf_new(filpsp, upf, ierr)
377  ABI_CHECK(ierr == 0, sjoin("read_upf_new returned ierr:", itoa(ierr)))
378 
379  call atomdata_from_symbol(atom, upf%psd)
380  znucl = atom%znucl
381  zion = upf%zp
382  mmax = upf%mesh
383  ! FIXME Temporary hack
384  !mmax = 300  ! H
385  !mmax = 600  ! Si
386  maxrad = upf%rmax
387  !maxrad = upf%r(mmax)
388 
389  ABI_CHECK(upfdft_to_ixc(upf%dft, pspxc, msg) == 0, msg)
390  lmax = upf%lmax
391 
392  ! Write some description of file
393  write(msg, '(3(a,1x))' ) '-',trim(upf%psd), trim(upf%generated)
394  call wrtout([std_out, ab_out], msg)
395  write(msg,'(a,f9.5,f10.5,2x,a,t47,a)')'-',znucl,zion,trim(upf%date),'znucl, zion, pspdat'
396  call wrtout([std_out, ab_out], msg)
397  write(msg, '(5(i0,1x),t47,a)' ) 12, pspxc, lmax, upf%lloc, mmax,'pspcod,pspxc,lmax,lloc,mmax'
398  call wrtout([std_out, ab_out], msg)
399 
400  ! Check that rad grid is linear starting at zero
401  linear_mesh = .True.
402  amesh = upf%r(2) - upf%r(1); damesh = zero
403  do irad=2,mmax-1
404    damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1)))
405  end do
406  linear_mesh = damesh < tol8
407 
408  if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then
409    write(msg,'(3a)')&
410    'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
411    'Action: check your pseudopotential file.'
412    ABI_ERROR(msg)
413  end if
414 
415  ! Check if the local component is one of the angular momentum channels
416  ! effectively if one of the ll is absent from the NL projectors
417  ABI_MALLOC(found_l, (0:lmax))
418  found_l = .true.
419  do ll=0,lmax
420    if (any(upf%lll(1:upf%nbeta) == ll)) found_l(ll) = .false.
421  end do
422 
423  if (count(found_l) /= 1) then
424    lloc = -1
425  else
426    do ll=0,lmax
427      if (found_l(ll)) then
428        lloc = ll
429        exit
430      end if
431    end do
432  end if
433  ABI_FREE(found_l)
434  !FIXME: do something about lloc == -1
435 
436  ! convert vloc from Rydberg to Ha
437  upf%vloc = half * upf%vloc
438 
439  if (linear_mesh) then
440    call psp8lo(amesh, epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl, &
441                vlspl(:,1), upf%r(1:mmax), upf%vloc(1:mmax), yp1, ypn, zion)
442  else
443    call psp11lo(upf%rab(1:mmax), epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl,&
444                 vlspl(:,1), upf%r(1:mmax), upf%vloc(1:mmax), yp1, ypn, zion)
445  end if
446 
447  ! Fit spline to q^2 V(q) (Numerical Recipes subroutine)
448  ABI_MALLOC(work_space, (psps%mqgrid_vl))
449  ABI_MALLOC(work_spl, (psps%mqgrid_vl))
450 
451  call spline(psps%qgrid_vl, vlspl(:,1), psps%mqgrid_vl, yp1, ypn, work_spl)
452 
453  vlspl(:,2) = work_spl(:)
454  ABI_FREE(work_space)
455  ABI_FREE(work_spl)
456 
457  nproj_l = 0
458 
459  if (.not. upf%has_so) then
460    do iprj=1,upf%nbeta
461      ll = upf%lll(iprj)
462      nproj_l(ll+1) = nproj_l(ll+1) + 1
463    end do
464    write(msg, '(a,*(i6))' ) '     nproj',nproj_l
465    call wrtout([std_out, ab_out], msg)
466 
467    ! shape = dimekb  vs. shape = n_proj
468    ! convert from Rydberg to Ha
469    do ll=1,upf%nbeta
470      ekb(ll) = upf%dion(ll,ll) * half
471    end do
472 
473    ! this has to do the FT of the projectors to reciprocal space
474    ! allocate proj to avoid temporary copy.
475    ABI_MALLOC(proj, (mmax,1:upf%nbeta))
476    proj = upf%beta(1:mmax,1:upf%nbeta)
477 
478    call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, &
479                 upf%nbeta, proj, upf%lll(1:upf%nbeta), upf%kbeta(1:upf%nbeta), &
480                 psps%qgrid_ff, upf%r(1:mmax), upf%rab(1:mmax), psps%useylm)
481 
482    ! This to reproduce psp8in version with linear meshes.
483    ! Compute Vanderbilt-KB form factors and fit splines
484    call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, &
485                psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj)
486 
487    ABI_FREE(proj)
488 
489  else
490    call upf2_jl2srso(upf, my_nproj_l, my_nprojso_l, vsr, esr, vso, eso)
491 
492    write(msg, '(a,*(i6))' ) '     nproj',my_nproj_l
493    call wrtout([std_out, ab_out], msg)
494    write(msg, '(5x,a)' ) "spin-orbit psp"
495    call wrtout([std_out, ab_out], msg)
496    write(msg, '(5x,a,*(i6))' ) '   nprojso',my_nprojso_l
497    call wrtout([std_out, ab_out], msg)
498 
499    ABI_CALLOC(proj, (mmax, psps%lnmax))
500    pspindex = 0; iln=0; indlmn(:,:)=0
501    nso = 2
502 
503    if (psps%pspso(ipsp) == 0) then
504      write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,&
505      ' but internal variable pspso is 0'
506      ABI_COMMENT(msg)
507      nso = 1
508    end if
509 
510    do nn=1,nso
511      !do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
512      !  ll = ipsang-(nn-1)*lmax-1
513      do il=1,lmax+1
514        ll = il - 1
515        if (nn == 1) then
516          ! SR part
517          do iprj=1,my_nproj_l(ll)
518            iln = iln + 1
519            ekb(iln) = esr(iprj, il)
520            proj(:,iln) = vsr(:, iprj, il)
521            nproj_l(il) = my_nproj_l(ll)
522            kk = iprj
523            do mm=1,2*ll*psps%useylm+1
524              pspindex = pspindex + 1
525              indlmn(1,pspindex) = ll
526              indlmn(2,pspindex) = mm-ll*psps%useylm-1
527              indlmn(3,pspindex) = kk
528              indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm
529              indlmn(5,pspindex) = iln
530              indlmn(6,pspindex) = nn
531              !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex
532            end do
533          end do
534 
535        else
536          ! SOC part
537          if (ll == 0) cycle
538          do iprj=1,my_nprojso_l(ll)
539            iln = iln + 1
540            ekb(iln) = eso(iprj, il)
541            proj(:,iln) = vso(:,iprj,il)
542            ! Note ll in nproj_l i.e. the s channel in the SOC part is not included in nproj.
543            nproj_l(ll + psps%mpsang) = my_nprojso_l(ll)
544            kk = iprj !+ my_nproj_l(ll)
545            do mm=1,2*ll*psps%useylm+1
546              pspindex = pspindex + 1
547              indlmn(1,pspindex) = ll
548              indlmn(2,pspindex) = mm-ll*psps%useylm-1
549              indlmn(3,pspindex) = kk
550              indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm
551              indlmn(5,pspindex) = iln
552              indlmn(6,pspindex) = nn
553              !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex
554            end do
555          end do
556        end if
557 
558      end do ! il
559    end do ! nn
560 
561    ! This to reproduce psp8in version with linear meshes.
562    ! Compute Vanderbilt-KB form factors and fit splines
563    call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, &
564                psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj)
565 
566    ABI_FREE(proj)
567    ABI_FREE(vsr)
568    ABI_FREE(esr)
569    ABI_FREE(vso)
570    ABI_FREE(eso)
571    !ABI_WARNING("upf2_to_abinit: UPF2 with SOC")
572  end if
573 
574  ! if we find a core density, do something about it
575  ! rho_atc contains the nlcc density
576  ! rho_at contains the total density
577  xcccrc = zero; xccc1d = zero
578 
579  if (upf%nlcc) then
580    ABI_MALLOC(ff, (mmax))
581    ABI_MALLOC(ff1, (mmax))
582    ABI_MALLOC(ff2, (mmax))
583    ff(1:mmax) = upf%rho_atc(1:mmax) ! model core charge without derivative factor
584    !smooth_niter = 15 ! run 15 iterations of smoothing?
585    smooth_niter = 0   ! Don't smooth core charges to be consistent with the treatment done in psp8in
586 
587    ff1 = zero
588    call nderiv(one, ff, ff1, mmax, 1) ! first derivative
589    ff1(1:mmax) = ff1(1:mmax) / upf%rab(1:mmax)
590    call smooth(ff1, mmax, smooth_niter)
591 
592    ff2 = zero
593    call nderiv(one, ff1, ff2, mmax, 1) ! second derivative
594    ff2(1:mmax) = ff2(1:mmax) / upf%rab(1:mmax)
595    call smooth(ff2, mmax, smooth_niter)
596 
597    ! determine a good rchrg = xcccrc
598    do ir = mmax, 1, -1
599      !if (abs(ff(ir)) > 1.e-6) then
600      if (abs(ff(ir)) > tol20) then
601        xcccrc = upf%r(ir)
602        exit
603      end if
604    end do
605    !xcccrc = upf%r(mmax)
606 
607    ABI_MALLOC(rad_cc, (mmax))
608    rad_cc = upf%r(1:mmax)
609    rad_cc(1) = zero ! force this so that the core charge covers whole spline interval.
610 
611    call cc_derivatives(rad_cc, ff, ff1, ff2, mmax, psps%n1xccc, xcccrc, xccc1d)
612    !call psp8cc(mmax, psps%n1xccc, xcccrc, xccc1d)
613 
614    ABI_FREE(rad_cc)
615    ABI_FREE(ff)
616    ABI_FREE(ff1)
617    ABI_FREE(ff2)
618  end if ! nlcc present
619 
620  ! Read pseudo valence charge in real space on the linear mesh
621  ! and transform it to reciprocal space on a regular grid
622  ! TODO: Spline input data on linear mesh if not linear
623  ABI_MALLOC(ff, (mmax))
624  ff = upf%rho_at / four_pi
625  where (abs(upf%r) > tol16)
626    ff = ff / upf%r ** 2
627  else where
628    ff = zero
629  end where
630 
631  ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
632  call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh)
633  call nctab_eval_tvalespl(nctab, zion, mesh, ff, psps%mqgrid_vl, psps%qgrid_vl)
634 
635  nctab%num_tphi = upf%nwfc
636  if (upf%nwfc > 0) then
637    ! Store atomic wavefunctions and metadata in nctab, then compute form factors for spline.
638    ! NB: lchi is the radial part of the KS equation, multiplied by r.
639    ABI_MALLOC(awfc_indlmn, (6, upf%nwfc))
640    awfc_indlmn = huge(1)
641    atmwfc_lmax = -1
642 
643    do iwfc=1, upf%nwfc
644      !print *, "label: ", upf%els(iwfc)
645      !print *, "n", upf%nchi(iwfc)
646      !print *, "nn", upf%nn(iwfc)
647      !print *, "j", upf%jchi(iwfc)
648      !print *, "l", upf%lchi(iwfc)
649      !print *, "occ:", upf%oc(iwfc)
650      atmwfc_lmax = max(atmwfc_lmax, upf%lchi(iwfc))
651      ff = upf%chi(:, iwfc) ** 2; call simp_gen(intg, ff, mesh)
652      write(std_out, *)" wavefunction (before rescaling) integrates to: ",intg
653      upf%chi(:, iwfc) = upf%chi(:, iwfc) / sqrt(intg)
654      ff = upf%chi(:, iwfc) ** 2; call simp_gen(intg, ff, mesh)
655      write(std_out, *)" wavefunction (after rescaling) integrates to: ",intg
656 
657      ! NB: we only need ll (1), and iln (5) in psp8nl
658      awfc_indlmn(1, iwfc) = upf%lchi(iwfc)
659      awfc_indlmn(5, iwfc) = iwfc
660    end do
661 
662    ! All this sfree/remalloc stuff is for handling memory in multi dataset mode!
663    ABI_REMALLOC(nctab%tphi_qspl, (psps%mqgrid_ff, 2, upf%nwfc))
664    ABI_SFREE(nctab%tphi_n)
665    ABI_SFREE(nctab%tphi_l)
666    ABI_SFREE(nctab%tphi_occ)
667    ABI_SFREE(nctab%tphi_jtot)
668    call alloc_copy(upf%nchi, nctab%tphi_n)
669    call alloc_copy(upf%lchi, nctab%tphi_l)
670    call alloc_copy(upf%oc, nctab%tphi_occ)
671    nctab%has_jtot = upf%has_so
672    if (upf%has_so) call alloc_copy(upf%jchi, nctab%tphi_jtot)
673 
674    !call psp8nl(amesh, nctab%tphi_qspl, awfc_indlmn, atmwfc_lmax, upf%nwfc, upf%nwfc, mmax, &
675    !            psps%mqgrid_ff, psps%qgrid_ff, upf%r, upf%chi)
676 
677    !do iwfc=1, upf%nwfc
678    !  do iq=1,psps%mqgrid_ff
679    !    write(555, *) nctab%tphi_qspl(iq,:,iwfc)
680    !  end do
681    !end do
682 
683    !do iwfc=1, upf%nwfc
684    !  where (abs(upf%r) > tol16)
685    !    upf%chi(:,iwfc) = upf%chi(:,iwfc) / upf%r
686    !  else where
687    !    upf%chi(:,iwfc) = zero
688    !  end where
689    !end do
690 
691    call pawpsp_nl(nctab%tphi_qspl, awfc_indlmn, upf%nwfc, upf%nwfc, psps%mqgrid_ff, psps%qgrid_ff, mesh, upf%chi)
692 
693    !do iwfc=1, upf%nwfc
694    !  do iq=1,psps%mqgrid_ff
695    !    write(666, *) nctab%tphi_qspl(iq,:,iwfc)
696    !  end do
697    !end do
698 
699    ABI_FREE(awfc_indlmn)
700    !stop "nwfc"
701  end if ! upf%nwfc > 0
702 
703  ABI_FREE(ff)
704  call pawrad_free(mesh)
705  call deallocate_pseudo_upf(upf)
706 
707 end subroutine upf2_to_abinit

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 functions times r, on a real space grid
  proj_l = angular momentum channel for each projector
  proj_nr = max number of r-points used for each projector
  qgrid(mqgrid)=q-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)

SOURCE

740 subroutine psp11nl(ffspl,indlmn, mmax, lnmax, lmnmax, mqgrid, n_proj, &
741                    proj, proj_l, proj_nr, qgrid, r, drdi, useylm)
742 
743 !Arguments ------------------------------------
744 !scalars
745  integer,intent(in) :: mmax, lnmax, lmnmax, mqgrid, useylm, n_proj
746 !arrays
747  integer, intent(in) :: proj_l(n_proj)
748  integer, intent(in) :: proj_nr(n_proj)
749  integer, intent(out) :: indlmn(6,lmnmax)
750  real(dp),intent(in) :: r(mmax)
751  real(dp),intent(in) :: drdi(mmax)
752  real(dp),intent(in) :: proj(mmax,n_proj)
753  real(dp),intent(in) :: qgrid(mqgrid)
754  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
755 
756 !Local variables-------------------------------
757 !scalars
758  integer,parameter :: bessorder = 0 ! never calculate derivatives of bessel functions
759  integer :: iproj, nr, ll, llold, ipsang, i_indlmn
760  integer :: iproj_1l, ir, iq, mm
761  real(dp) :: res, arg, besfact, dummy, dummy2
762  real(dp), allocatable :: work(:)
763  character(len=500) :: msg
764 
765 !*************************************************************************
766 
767  ffspl = zero; indlmn = 0; i_indlmn = 0
768  llold = -1; iproj_1l = 1
769 
770  ! loop over all projectors
771  do iproj=1,n_proj
772    if (iproj > lmnmax) then
773      write(msg,'(a,2i0)') ' Too many projectors found. n_proj, lmnmax = ',n_proj, lmnmax
774      ABI_ERROR(msg)
775    end if
776 
777    nr = proj_nr(iproj)
778    ABI_MALLOC(work, (nr))
779    ll = proj_l(iproj)
780    if (ll < llold) then
781      ABI_ERROR('UPF projectors are not in order of increasing ll')
782    else if (ll == llold) then
783      iproj_1l = iproj_1l + 1
784    else
785      iproj_1l = 1
786      llold = ll
787    end if
788 
789    ! determine indlmn for this projector (keep in UPF order and enforce that they are in increasing ll)
790    ! indlmn(6,lmnmax,ntypat)
791    ! For each type of psp,
792    ! array giving l,m,n,lm,ln,spin for i=ln  (if useylm=0)
793    !                                or i=lmn (if useylm=1)
794    ! NB: spin is used for NC pseudos with SOC term: 1 if scalar term (spin diagonal), 2 if SOC term.
795    do mm=1,2*ll*useylm+1
796      i_indlmn = i_indlmn + 1
797      indlmn(1, i_indlmn) = ll
798      indlmn(2, i_indlmn) = mm-ll*useylm-1
799      indlmn(3, i_indlmn) = iproj_1l
800      indlmn(4, i_indlmn) = ll*ll+(1-useylm)*ll+mm
801      indlmn(5, i_indlmn) = iproj
802      indlmn(6, i_indlmn) = 1 !spin? FIXME: to get j for relativistic cases
803    end do
804 
805    ! FT projectors to reciprocal space q
806    do iq=1,mqgrid
807      arg = two_pi * qgrid(iq)
808 
809      ! FIXME: add semianalytic form for integral from 0 to first point
810      do ir=1,nr
811        call jbessel(besfact, dummy, dummy2, ll, bessorder, arg*r(ir))
812        work(ir) = drdi(ir) * besfact * proj(ir, iproj) * r(ir) !* r(ir)
813      end do
814      call ctrap (nr, work, one, res)
815      ffspl(iq, 1, iproj) = res
816    end do
817    ABI_FREE(work)
818  end do ! iproj
819 
820  ! add derivative of ffspl(:,1,:) for spline interpolation later
821  ABI_MALLOC(work, (mqgrid))
822  do ipsang=1,lnmax
823    call spline(qgrid,ffspl(:,1,ipsang),mqgrid,zero,zero,ffspl(:,2,ipsang))
824  end do
825  ABI_FREE(work)
826 
827 end subroutine psp11nl