Branch data 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 : : !!***
|
