Line data Source code
1 : !!****m* ABINIT/m_upf2abinit
2 : !! NAME
3 : !! m_upf2abinit
4 : !!
5 : !! FUNCTION
6 : !! Procedures to read NC pseudos in UPF1/UPF2 format and convert data into
7 : !! the internal ABINIT representation.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2009-2026 ABINIT group (MJV, MG, DRH)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
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
43 : !!***
44 :
45 : public :: upf1_to_abinit
46 : public :: upf2_to_abinit
47 : !!***
48 :
49 : contains
50 : !!***
51 :
52 : !!****f* ABINIT/upf1_to_abinit
53 : !! NAME
54 : !! upf1_to_abinit
55 : !!
56 : !! FUNCTION
57 : !! This routine wraps a call to a PWSCF module, which reads in
58 : !! a UPF1 (PWSCF / Espresso) format pseudopotential, then transfers
59 : !! data to abinit internal variables.
60 : !! "UPF1 PWSCF format" (pspcod=11)
61 : !!
62 : !! INPUTS
63 : !! filpsp = name of file with UPF data
64 : !! psps = sturcture with global dimension data for pseudopotentials, header info ...
65 : !! used contents:
66 : !! psps%lmnmax
67 : !! psps%mqgrid_ff
68 : !! psps%mqgrid_vl
69 : !! psps%dimekb
70 : !! psps%n1xccc
71 : !! psps%qgrid_ff
72 : !! psps%qgrid_vl
73 : !!
74 : !! OUTPUT
75 : !! pspxc = index of xc functional for this pseudo
76 : !! lmax_ = maximal angular momentum
77 : !! lloc = local component chosen for pseudopotential
78 : !! mmax = maximum number of points in real space radial grid
79 : !! znucl = charge of species nucleus
80 : !! zion = valence charge
81 : !! epsatm = integral of local potential - coulomb potential of zion
82 : !! xcccrc = radius for non linear core correction
83 : !! ekb(dimekb)= Kleinman Bylander energies, see pspatm.F90
84 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
85 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
86 : !! or i=lmn (if useylm=1)
87 : !! vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
88 : !! ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
89 : !! second derivative from spline fit for each angular momentum and
90 : !! each projector; if any, spin-orbit components begin at l=mpsang+1
91 : !! nproj= number of projectors for each channel
92 : !! xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives,
93 : !! from psp file (used in NC only)
94 : !!
95 : !! SOURCE
96 :
97 1 : subroutine upf1_to_abinit(filpsp, znucl, zion, pspxc, lmax_, lloc, mmax, &
98 1 : 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 1 : logical, allocatable :: found_l(:)
125 1 : real(dp), allocatable :: work_spl(:), ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:)
126 :
127 : ! *********************************************************************
128 :
129 : ! ######### in module pseudo: ############
130 : !
131 : ! only npsx = 1 is used here
132 : ! grids are allocated for much larger fixed length (ndm=2000)
133 : ! number of species (6) and projectors (8) as well...
134 : !
135 : ! psd(npsx) = specied string
136 : ! pseudotype = uspp / nc string
137 : ! dft(npsx) = exchange correlation string (20 chars)
138 : ! lmax(npsx) = maximum l channel
139 : ! mesh(npsx) = number of points for local pot
140 : ! nbeta(npsx) = number of projectors (beta functions for uspp)
141 : ! nlcc(npsx) = flag for presence of NL core correction
142 : ! zp(npsx) = valence ionic charge
143 : ! r(ndm,npsx) = radial mesh
144 : ! rab(ndm,npsx) = dr / di for radial mesh
145 : ! rho_atc(ndm,npsx) = NLCC pseudocharge density
146 : ! rho_at(ndm,npsx) = 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 1 : if (open_file (filpsp,msg,newunit=iunit,status='old',form='formatted') /= 0) then
156 0 : ABI_ERROR(msg)
157 : end if
158 :
159 : ! read in psp data to static data in pseudo module, for ipsx == 1
160 1 : call read_pseudo(1,iunit)
161 1 : close (iunit)
162 :
163 : ! convert from Rydberg to Ha units
164 12007 : vloc0 = half * vloc0
165 439 : dion = half * dion
166 :
167 : ! if upf file is a USPP one, stop
168 1 : if (pseudotype == 'US') then
169 0 : 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 1 : call upfxc2abi(dft(1), pspxc)
176 1 : 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 3 : ABI_MALLOC(found_l, (0:lmax_))
181 5 : found_l = .true.
182 5 : do ll = 0, lmax_
183 11 : if (any(lll(1:nbeta(1),1) == ll)) found_l(ll) = .false.
184 : end do
185 :
186 5 : if (count(found_l) /= 1) then
187 0 : lloc = -1
188 : else
189 3 : do ll = 0, lmax_
190 3 : if (found_l(ll)) then
191 1 : lloc = ll
192 1 : exit
193 : end if
194 : end do
195 : end if
196 :
197 1 : ABI_FREE(found_l)
198 : !FIXME: do something about lloc == -1
199 :
200 1 : call atomdata_from_symbol(atom,psd(1))
201 1 : znucl = atom%znucl
202 1 : zion = zp(1)
203 1 : mmax = mesh(1)
204 :
205 : call psp11lo(rab(1:mmax,1), epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl, &
206 1 : 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 3 : ABI_MALLOC(work_spl, (psps%mqgrid_vl))
211 :
212 1 : call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl)
213 :
214 3002 : vlspl(:,2) = work_spl(:)
215 1 : ABI_FREE(work_spl)
216 :
217 : ! this has to do the FT of the projectors to reciprocal space
218 : ! allocate proj to avoid temporary copy.
219 4 : ABI_MALLOC(proj, (mmax,1:nbeta(1)))
220 1805 : proj = betar(1:mmax,1:nbeta(1), 1)
221 :
222 : call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, &
223 : nbeta(1), proj, lll(1:nbeta(1),1), ikk2(1:nbeta(1),1), &
224 1 : psps%qgrid_ff, r(1:mmax,1), rab(1:mmax,1), psps%useylm)
225 :
226 1 : ABI_FREE(proj)
227 :
228 5 : nproj_l = 0
229 4 : do iproj = 1, nbeta(1)
230 3 : ll = lll(iproj,1)
231 4 : nproj_l(ll+1) = nproj_l(ll+1) + 1
232 : end do
233 :
234 : ! shape = dimekb vs. shape = n_proj
235 4 : do ll = 1, nbeta(1)
236 4 : ekb(ll) = dion(ll,ll,1)
237 : end do
238 :
239 15013 : xcccrc = zero; xccc1d = zero
240 : ! if we find a core density, do something about it
241 : ! rho_atc contains the nlcc density
242 : ! rho_at contains the total density
243 1 : if (nlcc(1)) then
244 3 : ABI_MALLOC(ff, (mmax))
245 2 : ABI_MALLOC(ff1, (mmax))
246 2 : ABI_MALLOC(ff2, (mmax))
247 601 : ff(1:mmax) = rho_atc(1:mmax,1) ! model core charge without derivative factor
248 :
249 601 : ff1 = zero
250 1 : call nderiv(one,ff,ff1,mmax,1) ! first derivative
251 601 : ff1(1:mmax) = ff1(1:mmax) / rab(1:mmax,1)
252 1 : call smooth(ff1, mmax, 15) ! run 15 iterations of smoothing
253 :
254 601 : ff2 = zero
255 1 : call nderiv(one, ff1, ff2, mmax, 1) ! second derivative
256 601 : ff2(1:mmax) = ff2(1:mmax) / rab(1:mmax,1)
257 1 : call smooth(ff2, mmax, 15) ! run 15 iterations of smoothing?
258 :
259 : ! determine a good rchrg = xcccrc
260 147 : do ir = mmax, 1, -1
261 147 : if (abs(ff(ir)) > 1.e-6) then
262 1 : xcccrc = r(ir,1)
263 1 : exit
264 : end if
265 : end do
266 3 : ABI_MALLOC(rad_cc, (mmax))
267 602 : rad_cc = r(1:mmax,1)
268 1 : rad_cc(1) = zero ! force this so that the core charge covers whole spline interval.
269 :
270 1 : call cc_derivatives(rad_cc,ff,ff1,ff2,mmax,psps%n1xccc,xcccrc,xccc1d)
271 :
272 1 : ABI_FREE(rad_cc)
273 1 : ABI_FREE(ff)
274 1 : ABI_FREE(ff1)
275 1 : ABI_FREE(ff2)
276 : end if ! nlcc present
277 :
278 2 : end subroutine upf1_to_abinit
279 : !!***
280 :
281 : !!****f* ABINIT/upf2_to_abinit
282 : !! NAME
283 : !! upf2_to_abinit
284 : !!
285 : !! FUNCTION
286 : !! This routine wraps a call to a PWSCF module, which reads in
287 : !! a UPF2 (PWSCF / Espresso) format pseudopotential, then transfers
288 : !! data to abinit internal variables.
289 : !! "UPF2 PWSCF format" (pspcod=12)
290 : !!
291 : !! INPUTS
292 : !! filpsp = name of file with UPF2 data
293 : !! vloc_rcut= Real-space cutoff for local part
294 : !! psps = sturcture with global dimension data for pseudopotentials, header info ...
295 : !! used contents:
296 : !! psps%lmnmax
297 : !! psps%mqgrid_ff
298 : !! psps%mqgrid_vl
299 : !! psps%dimekb
300 : !! psps%n1xccc
301 : !! psps%qgrid_ff
302 : !! psps%qgrid_vl
303 : !!
304 : !! OUTPUT
305 : !! pspxc = index of xc functional for this pseudo
306 : !! lmax_ = maximal angular momentum
307 : !! lloc = local component chosen for pseudopotential
308 : !! mmax = maximum number of points in real space radial grid
309 : !! znucl = charge of species nucleus
310 : !! zion = valence charge
311 : !! epsatm = integral of local potential - coulomb potential of zion
312 : !! xcccrc = radius for non linear core correction
313 : !! ekb(dimekb)= Kleinman Bylander energies, see pspatm.F90
314 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
315 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
316 : !! or i=lmn (if useylm=1)
317 : !! vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
318 : !! ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
319 : !! second derivative from spline fit for each angular momentum and
320 : !! each projector; if any, spin-orbit components begin at l=mpsang+1
321 : !! nproj= number of projectors for each channel
322 : !! xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives,
323 : !! from psp file (used in NC only)
324 : !! xcctau1d(n1xccc*(1-usepaw),6)=1D core charge kinetic energy density function and five derivatives,
325 : !! from psp file (used in NC only)
326 : !! nctab<nctab_t>=NC tables
327 : !! %has_tvale=True if the pseudo contains the pseudo valence charge
328 : !! %tvalespl(mqgrid_vl,2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
329 : !!
330 : !! SOURCE
331 :
332 10 : subroutine upf2_to_abinit(ipsp, filpsp, vloc_rcut, znucl, zion, pspxc, lmax, lloc, mmax, &
333 10 : psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj_l, vlspl, xccc1d, xcctau1d, nctab, maxrad)
334 :
335 : use pseudo_types, only : pseudo_upf, deallocate_pseudo_upf !, pseudo_config
336 : use read_upf_new_module, only : read_upf_new
337 : use defs_datatypes, only : nctab_t
338 : use m_psps, only : nctab_eval_tvalespl
339 : use m_pspheads, only : upfdft_to_ixc, upf2_jl2srso
340 :
341 : !Arguments -------------------------------
342 : integer,intent(in) :: ipsp
343 : character(len=fnlen), intent(in) :: filpsp
344 : real(dp),intent(in) :: vloc_rcut
345 : type(pseudopotential_type),intent(in) :: psps
346 : type(nctab_t),intent(inout) :: nctab
347 : integer, intent(out) :: pspxc, lmax, lloc, mmax
348 : real(dp), intent(out) :: znucl, zion, epsatm, xcccrc, maxrad
349 : !arrays
350 : integer, intent(out) :: indlmn(6,psps%lmnmax), nproj_l(psps%mpssoang)
351 : real(dp), intent(inout) :: ekb(psps%dimekb)
352 : real(dp), intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax)
353 : real(dp), intent(out) :: vlspl(psps%mqgrid_vl,2)
354 : real(dp), intent(inout) :: xccc1d(psps%n1xccc,6)
355 : real(dp), intent(inout) :: xcctau1d(psps%n1xccc,6)
356 :
357 : !Local variables -------------------------
358 : integer :: ierr, ir, irad, iprj, il, ll, smooth_niter, nso, nn, iln, kk, mm, pspindex, iwfc !, iq
359 : integer :: atmwfc_lmax, mmax_cut
360 : !real(dp),parameter :: vloc_rcut = 10.0_dp ! QE Value
361 : !real(dp),parameter :: vloc_rcut = 6.0_dp ! PseudoDojo value used to generate psp8 files
362 : real(dp) :: yp1, ypn, amesh, damesh, intg
363 : character(len=500) :: msg
364 : logical :: linear_mesh, debug
365 10 : type(pseudo_upf) :: upf
366 : type(atomdata_t) :: atom
367 10 : type(pawrad_type) :: mesh
368 : integer :: my_nproj_l(0:3), my_nprojso_l(1:3), units(2)
369 : !integer :: nproj_tmp(psps%mpssoang)
370 10 : integer,allocatable :: awfc_indlmn(:,:)
371 10 : logical,allocatable :: found_l(:)
372 10 : real(dp),allocatable :: work_spl(:), ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:), chi_tmp(:)
373 10 : real(dp),allocatable :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:)
374 :
375 : ! *********************************************************************
376 :
377 30 : units = [std_out, ab_out]
378 :
379 : ! See also https://github.com/QEF/qeschemas/blob/master/UPF/qe_pp-0.99.xsd
380 : ! and https://github.com/QEF/qeschemas/files/9497267/pp.md
381 10 : call read_upf_new(filpsp, upf, ierr)
382 10 : ABI_CHECK(ierr == 0, sjoin("read_upf_new returned ierr:", itoa(ierr)))
383 :
384 10 : call atomdata_from_symbol(atom, upf%psd)
385 10 : znucl = atom%znucl
386 10 : zion = upf%zp
387 10 : mmax = upf%mesh
388 10 : maxrad = upf%rmax
389 :
390 10 : ABI_CHECK(upfdft_to_ixc(upf%dft, pspxc, msg) == 0, msg)
391 10 : lmax = upf%lmax
392 :
393 : ! Write some description of file
394 10 : write(msg, '(3(a,1x))' ) '-',trim(upf%psd), trim(upf%generated)
395 10 : call wrtout(units, msg)
396 10 : write(msg,'(a,f9.5,f10.5,2x,a,t47,a)')'-',znucl,zion,trim(upf%date),'znucl, zion, pspdat'
397 10 : call wrtout(units, msg)
398 10 : write(msg, '(5(i0,1x),t47,a)' ) 12, pspxc, lmax, upf%lloc, mmax,'pspcod,pspxc,lmax,lloc,mmax'
399 10 : call wrtout(units, msg)
400 :
401 : ! Check that rad grid is linear starting at zero
402 10 : linear_mesh = .True.
403 10 : amesh = upf%r(2) - upf%r(1); damesh = zero
404 14456 : do irad=2,mmax-1
405 14456 : damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1)))
406 : end do
407 10 : linear_mesh = damesh < tol8
408 :
409 10 : if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then
410 : write(msg,'(3a)')&
411 0 : 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
412 0 : 'Action: check your pseudopotential file.'
413 0 : ABI_ERROR(msg)
414 : end if
415 :
416 : ! Check if the local component is one of the angular momentum channels
417 : ! effectively if one of the ll is absent from the NL projectors
418 30 : ABI_MALLOC(found_l, (0:lmax))
419 40 : found_l = .true.
420 40 : do ll=0,lmax
421 108 : if (any(upf%lll(1:upf%nbeta) == ll)) found_l(ll) = .false.
422 : end do
423 :
424 40 : if (count(found_l) /= 1) then
425 10 : lloc = -1
426 : else
427 0 : do ll=0,lmax
428 0 : if (found_l(ll)) then
429 0 : lloc = ll
430 0 : exit
431 : end if
432 : end do
433 : end if
434 10 : ABI_FREE(found_l)
435 : !FIXME: do something about lloc == -1
436 :
437 : ! convert vloc from Rydberg to Ha
438 14476 : upf%vloc = half * upf%vloc
439 :
440 : ! =================================================
441 : ! This comment is from q-e/Modules/read_pseudo.F90
442 : ! =================================================
443 :
444 : ! the radial grid is defined up to r(mesh) but we introduce
445 : ! an auxiliary variable msh to limit the grid up to rcut=10 a.u.
446 : ! This is used to cut off the numerical noise arising from the
447 : ! large-r tail in cases like the integration of V_loc-Z/r
448 : !
449 : ! NB: In QE, the default value for vloc_rcut is 10 Bohr.
450 10 : call wrtout(std_out, sjoin(" Cutting radial-mesh for vloc using vloc_rcut:", ftoa(vloc_rcut), "(Bohr)"))
451 10 : mmax_cut = mmax
452 8508 : do ir=1,upf%mesh
453 8508 : if (upf%r(ir) > vloc_rcut) then
454 : mmax_cut = ir
455 : ! msh is forced to be odd for simpson integration (maybe obsolete?)
456 6 : mmax_cut = 2 * ( (mmax_cut + 1) / 2) - 1
457 6 : exit
458 : end if
459 : end do
460 :
461 : !write(std_out,*)" UPF file with mmax: ", mmax, " with r_max:", upf%r(mmax)
462 : !write(std_out,*)" Using mmax_cut: ", mmax_cut, " with r_cut:", upf%r(mmax_cut)
463 :
464 : ! Note mmax_cut here.
465 10 : if (linear_mesh) then
466 : call psp8lo(amesh, epsatm, mmax_cut, psps%mqgrid_vl, psps%qgrid_vl, &
467 10 : vlspl(:,1), upf%r, upf%vloc, yp1, ypn, zion)
468 : else
469 : call psp11lo(upf%rab, epsatm, mmax_cut, psps%mqgrid_vl, psps%qgrid_vl,&
470 0 : vlspl(:,1), upf%r, upf%vloc, yp1, ypn, zion)
471 : end if
472 :
473 : ! Fit spline to q^2 V(q) (Numerical Recipes subroutine)
474 30 : ABI_MALLOC(work_spl, (psps%mqgrid_vl))
475 :
476 10 : call spline(psps%qgrid_vl, vlspl(:,1), psps%mqgrid_vl, yp1, ypn, work_spl)
477 :
478 30020 : vlspl(:,2) = work_spl(:)
479 10 : ABI_FREE(work_spl)
480 :
481 10 : debug = .False.!; debug = .True.
482 : if (debug) then
483 : write(std_out,*)'# Vlocal upf = '
484 : write(std_out,*)' amesh = ', amesh
485 : write(std_out,*)' epsatm = ', epsatm
486 : write(std_out,*)' mmax = ', mmax
487 : write(std_out,*)' mqgrid = ', psps%mqgrid_vl
488 : do ir = 1, psps%mqgrid_vl
489 : write(std_out,*)' qgrid = ', ir, psps%qgrid_vl(ir)
490 : enddo
491 : do ir = 1, psps%mqgrid_vl
492 : write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2)
493 : enddo
494 : write(std_out,*)
495 : do ir = 1, mmax
496 : write(std_out,*)' rad = ', upf%r(ir), upf%vloc(ir)
497 : enddo
498 : write(std_out,*)
499 : write(std_out,*)' yp1 = ', yp1
500 : write(std_out,*)' ypn = ', ypn
501 : write(std_out,*)' zion = ', zion
502 : !stop
503 : end if
504 :
505 48 : nproj_l = 0
506 :
507 10 : if (.not. upf%has_so) then
508 40 : do iprj=1,upf%nbeta
509 34 : ll = upf%lll(iprj)
510 40 : nproj_l(ll+1) = nproj_l(ll+1) + 1
511 : end do
512 6 : write(msg, '(a,*(i6))' ) ' nproj',nproj_l
513 6 : call wrtout(units, msg)
514 :
515 : ! shape = dimekb vs. shape = n_proj
516 : ! convert from Rydberg to Ha
517 40 : do ll=1,upf%nbeta
518 40 : ekb(ll) = upf%dion(ll,ll) * half
519 : end do
520 :
521 : ! this has to do the FT of the projectors to reciprocal space
522 : ! allocate proj to avoid temporary copy.
523 24 : ABI_MALLOC(proj, (mmax,1:upf%nbeta))
524 45562 : proj = upf%beta(1:mmax,1:upf%nbeta)
525 :
526 : call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, &
527 : upf%nbeta, proj, upf%lll(1:upf%nbeta), upf%kbeta(1:upf%nbeta), &
528 6 : psps%qgrid_ff, upf%r(1:mmax), upf%rab(1:mmax), psps%useylm)
529 :
530 : ! This to reproduce psp8in version with linear meshes.
531 : ! Compute Vanderbilt-KB form factors and fit splines
532 : call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, &
533 6 : psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj)
534 :
535 6 : ABI_FREE(proj)
536 :
537 : else
538 4 : call upf2_jl2srso(upf, my_nproj_l, my_nprojso_l, vsr, esr, vso, eso)
539 :
540 4 : write(msg, '(a,*(i6))' ) ' nproj',my_nproj_l
541 4 : call wrtout(units, msg)
542 4 : write(msg, '(5x,a)' ) "spin-orbit psp"
543 4 : call wrtout(units, msg)
544 4 : write(msg, '(5x,a,*(i6))' ) ' nprojso',my_nprojso_l
545 4 : call wrtout(units, msg)
546 :
547 98596 : ABI_CALLOC(proj, (mmax, psps%lnmax))
548 424 : pspindex = 0; iln=0; indlmn(:,:)=0
549 4 : nso = 2
550 :
551 4 : if (psps%pspso(ipsp) == 0) then
552 2 : write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,&
553 4 : ' but internal variable pspso is 0'
554 2 : ABI_COMMENT(msg)
555 2 : nso = 1
556 : end if
557 :
558 10 : do nn=1,nso
559 : !do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
560 : ! ll = ipsang-(nn-1)*lmax-1
561 28 : do il=1,lmax+1
562 18 : ll = il - 1
563 24 : if (nn == 1) then
564 : ! SR part
565 42 : do iprj=1,my_nproj_l(ll)
566 30 : iln = iln + 1
567 30 : ekb(iln) = esr(iprj, il)
568 48722 : proj(:,iln) = vsr(:, iprj, il)
569 30 : nproj_l(il) = my_nproj_l(ll)
570 30 : kk = iprj
571 72 : do mm=1,2*ll*psps%useylm+1
572 30 : pspindex = pspindex + 1
573 30 : indlmn(1,pspindex) = ll
574 30 : indlmn(2,pspindex) = mm-ll*psps%useylm-1
575 30 : indlmn(3,pspindex) = kk
576 30 : indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm
577 30 : indlmn(5,pspindex) = iln
578 60 : indlmn(6,pspindex) = nn
579 : !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex
580 : end do
581 : end do
582 :
583 : else
584 : ! SOC part
585 6 : if (ll == 0) cycle
586 18 : do iprj=1,my_nprojso_l(ll)
587 14 : iln = iln + 1
588 14 : ekb(iln) = eso(iprj, il)
589 23002 : proj(:,iln) = vso(:,iprj,il)
590 : ! Note ll in nproj_l i.e. the s channel in the SOC part is not included in nproj.
591 14 : nproj_l(ll + psps%mpsang) = my_nprojso_l(ll)
592 14 : kk = iprj !+ my_nproj_l(ll)
593 34 : do mm=1,2*ll*psps%useylm+1
594 14 : pspindex = pspindex + 1
595 14 : indlmn(1,pspindex) = ll
596 14 : indlmn(2,pspindex) = mm-ll*psps%useylm-1
597 14 : indlmn(3,pspindex) = kk
598 14 : indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm
599 14 : indlmn(5,pspindex) = iln
600 28 : indlmn(6,pspindex) = nn
601 : !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex
602 : end do
603 : end do
604 : end if
605 :
606 : end do ! il
607 : end do ! nn
608 :
609 : ! This to reproduce psp8in version with linear meshes.
610 : ! Compute Vanderbilt-KB form factors and fit splines
611 : call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, &
612 4 : psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj)
613 :
614 4 : ABI_FREE(proj)
615 4 : ABI_FREE(vsr)
616 4 : ABI_FREE(esr)
617 4 : ABI_FREE(vso)
618 4 : ABI_FREE(eso)
619 : !ABI_WARNING("upf2_to_abinit: UPF2 with SOC")
620 : end if
621 :
622 : ! In Abinit, at least for the Troullier-Martins pseudopotential,
623 : ! the pseudocore charge density and its derivatives (xccc1d)
624 : ! are introduced in a linear grid.
625 : ! This grid is normalized, so the radial coordinates run between
626 : ! from 0 and 1 (from 0 to xcccrc, where xcccrc is the radius
627 : ! where the pseudo-core becomes zero).
628 :
629 300250 : xcccrc = zero; xccc1d = zero; xcctau1d = zero
630 :
631 10 : if (upf%nlcc) then
632 30 : ABI_MALLOC(ff, (mmax))
633 20 : ABI_MALLOC(ff1, (mmax))
634 20 : ABI_MALLOC(ff2, (mmax))
635 : ! model core charge without derivative factor
636 14476 : ff(1:mmax) = upf%rho_atc(1:mmax)
637 : !smooth_niter = 15 ! run 15 iterations of smoothing?
638 10 : smooth_niter = 0 ! Don't smooth core charges to be consistent with the treatment done in psp8in
639 :
640 14476 : ff1 = zero
641 10 : call nderiv(one, ff, ff1, mmax, 1) ! first derivative
642 14476 : ff1(1:mmax) = ff1(1:mmax) / upf%rab(1:mmax)
643 10 : call smooth(ff1, mmax, smooth_niter)
644 :
645 14476 : ff2 = zero
646 10 : call nderiv(one, ff1, ff2, mmax, 1) ! second derivative
647 14476 : ff2(1:mmax) = ff2(1:mmax) / upf%rab(1:mmax)
648 10 : call smooth(ff2, mmax, smooth_niter)
649 :
650 : ! determine a good rchrg = xcccrc
651 10728 : do ir = mmax, 1, -1
652 : !if (abs(ff(ir)) > tol6) then
653 10728 : if (abs(ff(ir)) > tol20) then
654 10 : xcccrc = upf%r(ir); exit
655 : end if
656 : end do
657 : !xcccrc = upf%r(mmax)
658 :
659 30 : ABI_MALLOC(rad_cc, (mmax))
660 14486 : rad_cc = upf%r(1:mmax)
661 10 : rad_cc(1) = zero ! force this so that the core charge covers whole spline interval.
662 :
663 10 : call cc_derivatives(rad_cc, ff, ff1, ff2, mmax, psps%n1xccc, xcccrc, xccc1d)
664 : !call psp8cc(mmax, psps%n1xccc, xcccrc, xccc1d)
665 :
666 10 : ABI_FREE(ff)
667 10 : ABI_FREE(ff1)
668 10 : ABI_FREE(ff2)
669 :
670 : ! Add corresponding splines for tau_mod, and flag to test for its presence before calculation
671 : ! NB: there are no r**2 or 4 pi factors in the tau quantities from Don Hamann in oncvpsp metagga
672 30 : ABI_MALLOC(ff, (mmax))
673 20 : ABI_MALLOC(ff1, (mmax))
674 20 : ABI_MALLOC(ff2, (mmax))
675 14476 : ff(1:mmax) = upf%tau_mod(1:mmax) ! model core charge without derivative factor
676 : !smooth_niter = 15 ! run 15 iterations of smoothing?
677 : smooth_niter = 0 ! Don't smooth core charges to be consistent with the treatment done in psp8in
678 :
679 14476 : ff1 = zero
680 10 : call nderiv(one, ff, ff1, mmax, 1) ! first derivative
681 14476 : ff1(1:mmax) = ff1(1:mmax) / upf%rab(1:mmax)
682 10 : call smooth(ff1, mmax, smooth_niter)
683 :
684 14476 : ff2 = zero
685 10 : call nderiv(one, ff1, ff2, mmax, 1) ! second derivative
686 14476 : ff2(1:mmax) = ff2(1:mmax) / upf%rab(1:mmax)
687 10 : call smooth(ff2, mmax, smooth_niter)
688 :
689 : ! use same xcccrc and rad_cc as for density above ??
690 :
691 10 : call cc_derivatives(rad_cc, ff, ff1, ff2, mmax, psps%n1xccc, xcccrc, xcctau1d)
692 :
693 10 : ABI_FREE(ff)
694 10 : ABI_FREE(ff1)
695 10 : ABI_FREE(ff2)
696 :
697 10 : ABI_FREE(rad_cc)
698 :
699 : end if ! nlcc present
700 :
701 : ! Read pseudo valence charge in real space on the linear mesh
702 : ! and transform it to reciprocal space on a regular grid
703 : ! TODO: Spline input data on linear mesh if not linear
704 30 : ABI_MALLOC(ff, (mmax))
705 14486 : ff = upf%rho_at(1:mmax) / four_pi
706 43408 : where (abs(upf%r) > tol16)
707 : ff = ff / upf%r(1:mmax) ** 2
708 : else where
709 : ff = zero
710 : end where
711 :
712 : ! TODO: use the tau_at pseudo atomic kinetic energy density as well, in upf%tau_at
713 : ! would also need to be splined as below for the atomic pseudo charge
714 : ! NB: there are no r**2 or 4 pi factors in the tau quantities from Don Hamann in oncvpsp metagga
715 : ! TODO: as well, nctab contents with the form factors of xccc1d and xcctau1d could be calculated here instead of in m_pspini
716 :
717 : ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
718 10 : call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh)
719 10 : call nctab_eval_tvalespl(nctab, zion, mesh, ff, psps%mqgrid_vl, psps%qgrid_vl)
720 :
721 10 : nctab%num_tphi = upf%nwfc
722 10 : if (upf%nwfc > 0) then
723 : ! Store atomic wavefunctions and metadata in nctab, then compute form factors for spline.
724 : ! NB: lchi is the radial part of the KS equation, multiplied by r.
725 30 : ABI_MALLOC(awfc_indlmn, (6, upf%nwfc))
726 248 : awfc_indlmn = huge(1)
727 10 : atmwfc_lmax = -1
728 30 : ABI_MALLOC(chi_tmp, (upf%mesh))
729 :
730 44 : do iwfc=1, upf%nwfc
731 : !print *, "label: ", upf%els(iwfc); print *, "n", upf%nchi(iwfc)
732 : !print *, "nn", upf%nn(iwfc); print *, "j", upf%jchi(iwfc)
733 : !print *, "l", upf%lchi(iwfc); print *, "occ:", upf%oc(iwfc)
734 34 : atmwfc_lmax = max(atmwfc_lmax, upf%lchi(iwfc))
735 51736 : chi_tmp = upf%chi(:, iwfc) ** 2; call simp_gen(intg, chi_tmp, mesh)
736 : !write(std_out, *)" wavefunction (before rescaling) integrates to: ",intg
737 51702 : upf%chi(:, iwfc) = upf%chi(:, iwfc) / sqrt(intg)
738 51736 : chi_tmp = upf%chi(:, iwfc) ** 2; call simp_gen(intg, chi_tmp, mesh)
739 : !write(std_out, *)" wavefunction (after rescaling) integrates to: ",intg
740 :
741 : ! NB: we only need ll (1), and iln (5) in psp8nl
742 34 : awfc_indlmn(1, iwfc) = upf%lchi(iwfc)
743 78 : awfc_indlmn(5, iwfc) = iwfc
744 : end do
745 10 : ABI_FREE(chi_tmp)
746 :
747 : ! All this sfree/remalloc stuff is for handling memory in multi dataset mode!
748 40 : ABI_REMALLOC(nctab%tphi_qspl, (psps%mqgrid_ff, 2, upf%nwfc))
749 10 : ABI_SFREE(nctab%tphi_n)
750 10 : ABI_SFREE(nctab%tphi_l)
751 10 : ABI_SFREE(nctab%tphi_occ)
752 10 : ABI_SFREE(nctab%tphi_jtot)
753 10 : call alloc_copy(upf%nchi, nctab%tphi_n)
754 10 : call alloc_copy(upf%lchi, nctab%tphi_l)
755 10 : call alloc_copy(upf%oc, nctab%tphi_occ)
756 10 : nctab%has_jtot = upf%has_so
757 10 : if (upf%has_so) call alloc_copy(upf%jchi, nctab%tphi_jtot)
758 :
759 : !call psp8nl(amesh, nctab%tphi_qspl, awfc_indlmn, atmwfc_lmax, upf%nwfc, upf%nwfc, mmax, &
760 : ! psps%mqgrid_ff, psps%qgrid_ff, upf%r, upf%chi)
761 :
762 : !do iwfc=1, upf%nwfc
763 : ! do iq=1,psps%mqgrid_ff
764 : ! write(555, *) nctab%tphi_qspl(iq,:,iwfc)
765 : ! end do
766 : !end do
767 :
768 : !do iwfc=1, upf%nwfc
769 : ! where (abs(upf%r) > tol16)
770 : ! upf%chi(:,iwfc) = upf%chi(:,iwfc) / upf%r
771 : ! else where
772 : ! upf%chi(:,iwfc) = zero
773 : ! end where
774 : !end do
775 :
776 10 : call pawpsp_nl(nctab%tphi_qspl, awfc_indlmn, upf%nwfc, upf%nwfc, psps%mqgrid_ff, psps%qgrid_ff, mesh, upf%chi)
777 :
778 : !do iwfc=1, upf%nwfc
779 : ! do iq=1,psps%mqgrid_ff
780 : ! write(666, *) nctab%tphi_qspl(iq,:,iwfc)
781 : ! end do
782 : !end do
783 :
784 10 : ABI_FREE(awfc_indlmn)
785 : !stop "nwfc"
786 : end if ! upf%nwfc > 0
787 :
788 10 : ABI_FREE(ff)
789 10 : call pawrad_free(mesh)
790 10 : call deallocate_pseudo_upf(upf)
791 :
792 30 : end subroutine upf2_to_abinit
793 : !!***
794 :
795 : !!****f* m_upf2abinit/psp11nl
796 : !! NAME
797 : !! psp11nl
798 : !!
799 : !! FUNCTION
800 : !! Fourier transform the real space UPF projector functions to reciprocal space
801 : !!
802 : !! INPUTS
803 : !! lmax=maximum ang momentum for which nonlocal form factor is desired.
804 : !! Usually lmax=1, sometimes = 0 (e.g. for oxygen); lmax <= 2 allowed.
805 : !! mmax=number of radial grid points for atomic grid
806 : !! lnmax= maximum index for all l channel projectors, dimension of ffspl
807 : !! lmnmax= maximum index for all projectors, dimension of indlmn
808 : !! mqgrid=number of grid points for q grid
809 : !! n_proj = total number of NL projectors read in
810 : !! proj = projector functions times r, on a real space grid
811 : !! proj_l = angular momentum channel for each projector
812 : !! proj_nr = max number of r-points used for each projector
813 : !! qgrid(mqgrid)=q-values at which form factors are returned
814 : !! r(mmax)=radial grid values
815 : !! drdi=derivative of grid point wrt index
816 : !! useylm = input to use m dependency of NL part, or only Legendre polynomials
817 : !!
818 : !! OUTPUT
819 : !! ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and
820 : !! second derivative from spline fit for each angular momentum
821 : !! indlmn = indexing of each projector, for n, l, m, s, ln, lmn (see pspatm.F90)
822 : !!
823 : !! SOURCE
824 :
825 7 : subroutine psp11nl(ffspl,indlmn, mmax, lnmax, lmnmax, mqgrid, n_proj, &
826 7 : proj, proj_l, proj_nr, qgrid, r, drdi, useylm)
827 :
828 : !Arguments ------------------------------------
829 : !scalars
830 : integer,intent(in) :: mmax, lnmax, lmnmax, mqgrid, useylm, n_proj
831 : !arrays
832 : integer, intent(in) :: proj_l(n_proj)
833 : integer, intent(in) :: proj_nr(n_proj)
834 : integer, intent(out) :: indlmn(6,lmnmax)
835 : real(dp),intent(in) :: r(mmax)
836 : real(dp),intent(in) :: drdi(mmax)
837 : real(dp),intent(in) :: proj(mmax,n_proj)
838 : real(dp),intent(in) :: qgrid(mqgrid)
839 : real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
840 :
841 : !Local variables-------------------------------
842 : !scalars
843 : integer,parameter :: bessorder0 = 0 ! never calculate derivatives of bessel functions
844 : integer :: iproj, nr, ll, llold, ipsang, i_indlmn
845 : integer :: iproj_1l, ir, iq, mm
846 : real(dp) :: res, arg, besfact, dummy, dummy2
847 7 : real(dp), allocatable :: work(:)
848 : character(len=500) :: msg
849 :
850 : !*************************************************************************
851 :
852 240662 : ffspl = zero; indlmn = 0; i_indlmn = 0
853 : llold = -1; iproj_1l = 1
854 :
855 : ! loop over all projectors
856 44 : do iproj=1,n_proj
857 37 : if (iproj > lmnmax) then
858 0 : write(msg,'(a,2i0)') ' Too many projectors found. n_proj, lmnmax = ',n_proj, lmnmax
859 0 : ABI_ERROR(msg)
860 : end if
861 :
862 37 : nr = proj_nr(iproj)
863 111 : ABI_MALLOC(work, (nr))
864 37 : ll = proj_l(iproj)
865 37 : if (ll < llold) then
866 0 : ABI_ERROR('UPF projectors are not in order of increasing ll')
867 37 : else if (ll == llold) then
868 16 : iproj_1l = iproj_1l + 1
869 : else
870 : iproj_1l = 1
871 : llold = ll
872 : end if
873 :
874 : ! determine indlmn for this projector (keep in UPF order and enforce that they are in increasing ll)
875 : ! indlmn(6,lmnmax,ntypat)
876 : ! For each type of psp,
877 : ! array giving l,m,n,lm,ln,spin for i=ln (if useylm=0)
878 : ! or i=lmn (if useylm=1)
879 : ! NB: spin is used for NC pseudos with SOC term: 1 if scalar term (spin diagonal), 2 if SOC term.
880 94 : do mm=1,2*ll*useylm+1
881 57 : i_indlmn = i_indlmn + 1
882 57 : indlmn(1, i_indlmn) = ll
883 57 : indlmn(2, i_indlmn) = mm-ll*useylm-1
884 57 : indlmn(3, i_indlmn) = iproj_1l
885 57 : indlmn(4, i_indlmn) = ll*ll+(1-useylm)*ll+mm
886 57 : indlmn(5, i_indlmn) = iproj
887 94 : indlmn(6, i_indlmn) = 1 !spin? FIXME: to get j for relativistic cases
888 : end do
889 :
890 : ! FT projectors to reciprocal space q
891 111074 : do iq=1,mqgrid
892 111037 : arg = two_pi * qgrid(iq)
893 :
894 : ! FIXME: add semianalytic form for integral from 0 to first point
895 24551181 : do ir=1,nr
896 24440144 : call jbessel(besfact, dummy, dummy2, ll, bessorder0, arg*r(ir))
897 24551181 : work(ir) = drdi(ir) * besfact * proj(ir, iproj) * r(ir) !* r(ir)
898 : end do
899 111037 : call ctrap (nr, work, one, res)
900 111074 : ffspl(iq, 1, iproj) = res
901 : end do
902 44 : ABI_FREE(work)
903 : end do ! iproj
904 :
905 : ! add derivative of ffspl(:,1,:) for spline interpolation later
906 21 : ABI_MALLOC(work, (mqgrid))
907 47 : do ipsang=1,lnmax
908 47 : call spline(qgrid,ffspl(:,1,ipsang),mqgrid,zero,zero,ffspl(:,2,ipsang))
909 : end do
910 7 : ABI_FREE(work)
911 :
912 10 : end subroutine psp11nl
913 : !!***
914 :
915 : !!****f* ABINIT/psp11lo
916 : !! NAME
917 : !! psp11lo
918 : !!
919 : !! FUNCTION
920 : !! Compute sine transform to transform from V(r) to q^2 V(q).
921 : !! Computes integrals on logarithmic grid using related uniform
922 : !! grid in exponent and corrected trapezoidal integration.
923 : !! Generalized from psp5lo for non-log grids using dr/di weights.
924 : !!
925 : !! INPUTS
926 : !! drdi=derivative of radial grid wrt index
927 : !! mmax=number of radial r grid points
928 : !! mqgrid=number of grid points in q from 0 to qmax.
929 : !! qgrid(mqgrid)=q grid values (bohr**-1).
930 : !! rad(mmax)=r grid values (bohr).
931 : !! vloc(mmax)=V(r) on radial grid.
932 : !! zion=nominal valence charge of atom.
933 : !!
934 : !! OUTPUT
935 : !! epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
936 : !!{{\\ \begin{equation}
937 : !! q2vq(mqgrid)
938 : !! =q^2 V(q)
939 : !! = -\frac{Zv}{\pi}
940 : !! + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
941 : !!\end{equation} }}
942 : !! yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax (needed for spline fitter).
943 : !!
944 : !! SOURCE
945 :
946 1 : subroutine psp11lo(drdi,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion)
947 :
948 : !Arguments----------------------------------------------------------
949 : !scalars
950 : integer,intent(in) :: mmax,mqgrid
951 : real(dp),intent(in) :: zion
952 : real(dp),intent(out) :: epsatm,yp1,ypn
953 : !arrays
954 : real(dp),intent(in) :: drdi(mmax)
955 : real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
956 : real(dp),intent(out) :: q2vq(mqgrid)
957 :
958 : !Local variables-------------------------------
959 : !scalars
960 : integer :: iq,ir
961 : real(dp),parameter :: scale=10.0d0
962 : real(dp) :: arg,result_ctrap,test,ztor1
963 : !arrays
964 1 : real(dp),allocatable :: work(:)
965 :
966 : ! *************************************************************************
967 :
968 3 : ABI_MALLOC(work,(mmax))
969 :
970 : ! Do q=0 separately (compute epsatm)
971 : ! Do integral from 0 to r1
972 1 : ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
973 :
974 : ! Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
975 : ! with extra factor of drdi to convert to uniform grid
976 601 : do ir = 1, mmax
977 : ! First handle tail region
978 600 : test=vloc(ir)+zion/rad(ir)
979 : ! Ignore small contributions, or impose a cut-off in the case
980 : ! the pseudopotential data are in single precision.
981 : ! (it is indeed expected that vloc is very close to zero beyond 20,
982 : ! so a value larger than 2.0d-8 is considered anomalous)
983 600 : if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
984 70 : work(ir)=zero
985 : else
986 530 : work(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
987 : end if
988 601 : work(ir)=work(ir)*drdi(ir)
989 : end do
990 :
991 : ! Do integral from r(1) to r(max)
992 1 : call ctrap(mmax,work,one,result_ctrap)
993 1 : epsatm=4.d0*pi*(result_ctrap+ztor1)
994 :
995 1 : q2vq(1)=-zion/pi
996 :
997 : ! Loop over q values
998 3001 : do iq=2,mqgrid
999 3000 : arg=2.d0*pi*qgrid(iq)
1000 : ! ztor1=$ -Zv/\pi + 2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
1001 3000 : ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion) * cos(arg*rad(1)) )/pi
1002 :
1003 : ! set up integrand
1004 1803000 : do ir=1,mmax
1005 : !test=vloc(ir)+zion/rad(ir)
1006 : !Ignore contributions within decade of machine precision (suppressed ...)
1007 : !if ((scale+abs(test)).eq.scale) then
1008 : !work(ir)=zero
1009 : !else
1010 1800000 : work(ir)=sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
1011 : !end if
1012 1803000 : work(ir)=work(ir)*drdi(ir)
1013 : end do
1014 : ! do integral from r(1) to r(mmax)
1015 3000 : call ctrap(mmax,work,one,result_ctrap)
1016 :
1017 : ! store q^2 v(q)
1018 : ! FIXME: I only see one factor q, not q^2, but the same is done in other pspXlo.F90
1019 3001 : q2vq(iq)=ztor1+2.d0*qgrid(iq)*result_ctrap
1020 :
1021 : end do
1022 :
1023 : ! Compute derivatives of q^2 v(q) at ends of interval
1024 1 : yp1=0.0d0
1025 : !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
1026 : !integral from 0 to r1
1027 1 : arg=2.0d0*pi*qgrid(mqgrid)
1028 1 : ztor1=zion*rad(1)*sin(arg*rad(1))
1029 1 : 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))
1030 : !integral from r(mmax) to infinity is overkill; ignore
1031 : !set up integrand
1032 601 : do ir=1,mmax
1033 : !test=vloc(ir)+zion/rad(ir)
1034 : !Ignore contributions within decade of machine precision (supressed ...)
1035 : !if ((scale+abs(test)).eq.scale) then
1036 : !work(ir)=0.0d0
1037 : !else
1038 600 : work(ir)=(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * (rad(ir)*vloc(ir)+zion)
1039 : !end if
1040 601 : work(ir)=work(ir)*drdi(ir)
1041 : end do
1042 :
1043 1 : call ctrap(mmax,work,one,result_ctrap)
1044 1 : ypn=2.0d0 * (ztor1 + result_ctrap)
1045 1 : ABI_FREE(work)
1046 :
1047 1 : end subroutine psp11lo
1048 : !!***
1049 :
1050 : end module m_upf2abinit
1051 : !!***
|
