Line data Source code
1 : !!****m* ABINIT/m_psp9
2 : !! NAME
3 : !! m_psp9
4 : !!
5 : !! FUNCTION
6 : !! Initialize pspcod=9 (pseudopotentials from the PSML XML format):
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (JJ, MVer, YP)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_psp9
23 :
24 : use defs_basis
25 : use m_splines
26 : use m_errors
27 : use m_abicore
28 : #if defined HAVE_LIBPSML
29 : use m_psml
30 : #endif
31 :
32 : use defs_datatypes, only : nctab_t
33 : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free
34 : use m_psps, only : nctab_eval_tvalespl, nctab_eval_tvaletauspl
35 : use m_psptk, only : psp8lo, psp8nl
36 :
37 : implicit none
38 :
39 : private
40 : !!***
41 :
42 : public :: psp9in
43 : !!***
44 :
45 : contains
46 : !!***
47 :
48 : !!****f* ABINIT/psp9in
49 : !! NAME
50 : !! psp9in
51 : !!
52 : !! FUNCTION
53 : !! Initialize pspcod=9 (pseudopotentials from the PSML XML format):
54 : !! continue to read the corresponding file, then compute the
55 : !! local and non-local potentials.
56 : !!
57 : !! INPUTS
58 : !! filpsp=filename of the PSML pseudopotential
59 : !! lloc=angular momentum choice of local pseudopotential
60 : !! lmax=value of lmax mentioned at the second line of the psp file
61 : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
62 : !! =if useylm=0, max number of (l,n) comp. over all type of psps
63 : !! lnmax=max. number of (l,n) components over all type of psps
64 : !! mmax=maximum number of points in real space grid in the psp file
65 : !! angular momentum of nonlocal pseudopotential
66 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
67 : !! mpssoang= 2*maximum angular momentum for nonlocal pseudopotentials - 1
68 : !! mqgrid=dimension of q (or G) grid for arrays.
69 : !! mqgrid_vl=dimension of q (or G) grid for valence charge (array qgrid_vl)
70 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
71 : !! qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
72 : !! qgrid_vl(psps%mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for valence charge
73 : !! pspso=spin-orbit characteristics, govern the content of ffspl and ekb
74 : !! if =0 : this input requires NO spin-orbit characteristics of the psp
75 : !! if =2 : this input requires HGH or psp8 characteristics of the psp
76 : !! if =3 : this input requires HFN characteristics of the psp
77 : !! useylm=governs the way the nonlocal operator is to be applied:
78 : !! 1=using Ylm, 0=using Legendre polynomials
79 : !! zion=nominal valence of atom as specified in psp file
80 : !! znucl=nuclear number of atom as specified in psp file
81 : !!
82 : !! OUTPUT
83 : !! ekb(lnmax)=Kleinman-Bylander energy,
84 : !! {{\ \begin{equation}
85 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
86 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
87 : !! \end{equation} }}
88 : !! for each (l,n)
89 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree)
90 : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
91 : !! second derivative from spline fit for each angular momentum and
92 : !! each projector
93 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
94 : !! or i=lmn (if useylm=1)
95 : !! nproj(mpssoang)=number of projection functions for each angular momentum
96 : !! qchrg is not used, and could be suppressed later
97 : !! vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit
98 : !! xcccrc=XC core correction cutoff radius (bohr)
99 : !! xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
100 : !! nctab<nctab_t>=NC tables
101 : !! %has_tvale=True if the pseudo contains the pseudo valence charge
102 : !! %tvalespl(mqgrid_vl,2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
103 : !!
104 : !! SOURCE
105 :
106 15 : subroutine psp9in(filpsp,ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
107 15 : & mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,&
108 15 : & useylm,vlspl,xcccrc,xccc1d,xcctau1d,zion,znucl,nctab,maxrad)
109 :
110 : !Arguments ------------------------------------
111 : !scalars
112 : integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mpsang,mpssoang,mqgrid,mqgrid_vl
113 : integer,intent(in) :: pspso,n1xccc,useylm
114 : integer,intent(out) :: mmax
115 : real(dp),intent(in) :: zion,znucl
116 : real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad
117 : type(nctab_t),intent(inout) :: nctab
118 : character(len=fnlen),intent(in) :: filpsp
119 : !arrays
120 : integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang)
121 : real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl)
122 : real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
123 : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
124 : real(dp),intent(inout) :: xcctau1d(n1xccc,6)
125 :
126 : !Local variables-------------------------------
127 : !scalars
128 : #if defined HAVE_LIBPSML
129 : integer :: iln,pspindex,ipsang,irad,kk,ll
130 : integer :: mm,nn,nso,ii,ir,il
131 : integer :: nshells
132 : integer :: iproj,irelt,nders
133 : integer :: np_dn, np_lj, np_nr, np_so, np_sr, np_up, val_l, val_n
134 : real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn,zval
135 : real(dp) :: rmax,rmatch,z,chgvps
136 : real(dp) :: val_occ
137 : real(dp) :: rchrg_tau
138 : logical :: has_nlcc,has_spin,has_metagga
139 : logical :: has_tvale,oncvpsp
140 : character(len=500) :: message
141 : character(len=30) :: creator
142 : character(len=7), parameter :: oncvpsp_name = "ONCVPSP"
143 15 : type(pawrad_type) :: mesh
144 : #endif
145 : !arrays
146 : #if defined HAVE_LIBPSML
147 15 : integer, allocatable :: idx_so(:),idx_sr(:)
148 15 : real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_spl(:)
149 15 : real(dp),allocatable :: ff_cc(:),ff_tau(:)
150 15 : type(ps_t) :: psxml
151 : #endif
152 :
153 : ! ***************************************************************************
154 :
155 : #if defined HAVE_LIBPSML
156 :
157 15 : call ps_destroy(psxml)
158 15 : call psml_reader(filpsp,psxml,debug=.true.)
159 :
160 : !Identify the atomic code that generated the pseudopotential
161 15 : call ps_Provenance_Get(psxml, 1, creator=creator)
162 : !Check whether the pseudopotential has been created with ONCVPSP,
163 : !Don Hamann's code
164 : oncvpsp = (trim(creator(1:7)) .eq. trim(oncvpsp_name))
165 : !DEBUG
166 : !write(std_out,*)' psp9in : creator : ', creator
167 : !write(std_out,*)' psp9in : oncvpsp : ', oncvpsp
168 : !ENDDEBUG
169 :
170 : ! SIESTA's ATOM uses spherical harmonics, while ONCVPSP uses Legendre
171 : ! polynomials, which means we have to check the consistency of input variables
172 : ! wrt the pseudos
173 : !
174 : ! Note: commented because NC pseudos do not have non-diagonal terms
175 : !
176 : ! if ( oncvpsp ) then
177 : ! if ( useylm /= 0 ) then
178 : ! write(message,'(3a)') "ONCVPSP pseudos use Legendre polynomials but we use spherical harmonics", &
179 : !& ch10, "ACTION: set useylm to 0 in your input file"
180 : ! ABI_ERROR(message)
181 : ! endif
182 : ! else
183 : ! if ( useylm == 0 ) then
184 : ! write(message,'(3a)') "ATOM pseudos use spherical harmonics but we use Legendre polynomials", &
185 : !& ch10, "ACTION: set useylm to 1 in your input file"
186 : ! ABI_ERROR(message)
187 : ! endif
188 : ! endif
189 :
190 : ! The atomic number is a real number instead of a simple integer
191 : ! z (in Abinit), atomic-number in the header of the PSML file.
192 : ! z = ps_AtomicNumber(psxml)
193 : !
194 : ! The difference between the number of protons in the nucleus and the
195 : ! sum of the populations of the core shells is the effective atomic number
196 : ! of the pseudo-atom, Zval (in Abinit), z-pseudo in the header of the
197 : ! PSML file.
198 : ! zval = ps_Zpseudo(psxml)
199 :
200 : has_metagga = .false.
201 : #if defined HAVE_LIBPSML_METAGGA
202 : call ps_PseudoAtomSpec_Get(psxml, &
203 : & atomic_number=z, z_pseudo=zval, &
204 15 : & spin_dft=has_spin, core_corrections=has_nlcc, meta_gga=has_metagga)
205 : #else
206 : call ps_PseudoAtomSpec_Get(psxml, &
207 : & atomic_number=z, z_pseudo=zval, &
208 : & spin_dft=has_spin, core_corrections=has_nlcc)
209 : #endif
210 :
211 : !---
212 :
213 : !Feb 2015: shifted to Hamann grid for convenience - libpsml interpolates anyway
214 : !
215 : ! The following lines are taken from the oncvpsp.f90 subroutine of the oncvpsp
216 : ! code implemented by D. Hamann
217 : ! The atomic number of the element is read from the header of the XML file
218 : ! Logarithmic grid defined by Hamann in oncvpsp code
219 : ! z = psxml%header%z
220 : ! amesh = 1.012d0
221 : ! al = dlog(amesh)
222 : ! rr1 = .0005d0/z
223 : ! mmax = dlog(45.0d0 /rr1)/al
224 : !
225 : ! ABI_MALLOC( rad,(mmax) )
226 : !
227 : ! do ir = 1, mmax
228 : ! rad(ir) = rr1 * dexp(al*(ir-1))
229 : ! end do
230 :
231 : !Determine the maximum number of points in the grid ---
232 15 : rmax = 6.0_dp
233 15 : amesh = 0.01_dp
234 15 : mmax = int(rmax/amesh)
235 : ! if(mod(mmax,2) .eq. 0) mmax = mmax + 1
236 :
237 : !Print core charge info, for compatibility with psp8
238 15 : rchrg = zero
239 15 : fchrg = zero
240 15 : if (has_nlcc) then
241 9 : rchrg = amesh * (mmax - 2)
242 : ! PSML does not store fchrg for now but we know we have core corrections,
243 : ! then let's set it arbitrarily to 1.0
244 9 : fchrg = one
245 : else
246 6 : write(message, '(a)' ) '- psp9in: No XC core correction.'
247 6 : call wrtout(std_out,message,'COLL')
248 : end if
249 15 : write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,zero,'rchrg,fchrg,qchrg'
250 15 : call wrtout(ab_out,message,'COLL')
251 15 : call wrtout(std_out,message,'COLL')
252 :
253 : !Do we have a valence charge?
254 15 : call ps_ValenceConfiguration_Get(psxml, nshells=nshells)
255 15 : has_tvale = (nshells > 0)
256 :
257 : ! Compute the valence charge of the reference configuration used to
258 : ! generate the pseudopotential
259 15 : chgvps = 0.0_dp
260 : ! Loop on all the shells included in the valence
261 49 : do il = 1, nshells
262 : ! Sum the corresponding occupation of each shell
263 : ! FIXME: What if there is spin?
264 34 : call ps_ValenceShell_Get(psxml, il, n=val_n, l=val_l, occupation=val_occ)
265 34 : chgvps = chgvps + val_occ
266 34 : write(std_out,*)' psp9in : n, l, occupation = ', &
267 83 : & val_n, val_l, val_occ
268 : end do
269 :
270 : !DEBUG
271 : !write(std_out,*)' psp9in : atomic number'
272 : !write(std_out,*)' psp9in : z = ', z
273 : !write(std_out,*)' psp9in : valence charge of the reference configuration'
274 : !write(std_out,*)' psp9in : chgvps = ', chgvps
275 : !write(std_out,*)' psp9in : nominal valence charge'
276 : !write(std_out,*)' psp9in : zval = ', zval
277 : !write(std_out,*)' psp9in : mqgrid_vl = ', mqgrid_vl
278 : !write(std_out,*)' psp9in : parameters to define the points of the grid'
279 : !write(std_out,*)' psp9in : amesh = ', amesh
280 : !write(std_out,*)' psp9in : rmax = ', rmax
281 : !write(std_out,*)' psp9in : mmax = ', mmax
282 : !ENDDEBUG
283 :
284 : ! TODO: should be simple to average these and get difference for SREL+SOC,
285 : ! but also the Ekb etc...
286 15 : call ps_NonlocalProjectors_Filter(psxml, set=SET_DOWN, number=np_dn)
287 15 : call ps_NonlocalProjectors_Filter(psxml, set=SET_LJ, number=np_lj)
288 15 : call ps_NonlocalProjectors_Filter(psxml, set=SET_NONREL, number=np_nr)
289 15 : call ps_NonlocalProjectors_Filter(psxml, set=SET_SO, number=np_so)
290 15 : call ps_NonlocalProjectors_Filter(psxml, set=SET_SREL, number=np_sr)
291 15 : call ps_NonlocalProjectors_Filter(psxml, set=SET_UP, number=np_up)
292 15 : if (np_lj > 0) then
293 0 : message = 'For the moment LJ format projectors are not supported; SREL + SO is the internal abinit format'
294 0 : ABI_BUG(message)
295 : end if
296 :
297 15 : if (np_up > 0 .or. np_dn > 0) then
298 0 : write (message,'(3a)') 'For the moment separate spin up and down format projectors are not supported;',ch10,&
299 0 : & ' spin average is the internal abinit format'
300 0 : ABI_BUG(message)
301 : end if
302 :
303 : !--------------------------------------------------------------------
304 :
305 : !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn
306 15 : if(pspso==2) then
307 : nso=2
308 : else
309 15 : nso=1
310 : end if
311 :
312 : !Find the number of projectors per angular momentum shell
313 68 : nproj(:)=0
314 15 : if (np_nr > 0) then
315 6 : call ps_NonlocalProjectors_Filter(psxml, set=SET_NONREL, indexes=idx_sr)
316 32 : do iproj = 1, np_nr
317 26 : call ps_Projector_Get(psxml, idx_sr(iproj), l=il)
318 32 : nproj(il+1) = nproj(il+1) + 1
319 : end do
320 : else
321 9 : if (np_sr > 0) then
322 9 : call ps_NonlocalProjectors_Filter(psxml, set=SET_SREL, indexes=idx_sr)
323 60 : do iproj = 1, np_sr
324 51 : call ps_Projector_Get(psxml, idx_sr(iproj), l=il)
325 60 : nproj(il+1) = nproj(il+1) + 1
326 : end do
327 : else ! this should not happen
328 0 : ABI_BUG('Your psml potential should have either scalar- or non- relativistic projectors')
329 : end if
330 : end if
331 :
332 15 : write(message, '(a,5i6)' ) ' nproj',nproj(1:lmax+1)
333 15 : call wrtout(ab_out,message,'COLL')
334 15 : call wrtout(std_out, message,'COLL')
335 :
336 15 : irelt = 0
337 15 : if (nso == 2) then
338 0 : call ps_NonlocalProjectors_Filter(psxml, set=SET_SO, indexes=idx_so)
339 0 : do iproj = 1, np_so
340 0 : call ps_Projector_Get(psxml, idx_so(iproj), l=il)
341 0 : nproj(il+lmax+2) = nproj(il+lmax+2) + 1
342 0 : irelt = 1
343 : end do
344 : end if
345 :
346 554 : pspindex=0;iln=0;indlmn(:,:)=0
347 30 : do nn=1,nso
348 83 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
349 53 : ll=ipsang-(nn-1)*lmax-1
350 68 : if (nproj(ipsang)>0) then
351 130 : do kk=1,nproj(ipsang)
352 77 : iln=iln+1
353 207 : do mm=1,2*ll*useylm+1
354 77 : pspindex=pspindex+1
355 77 : indlmn(1,pspindex)=ll ! l angular momentum channel
356 77 : indlmn(2,pspindex)=mm-ll*useylm-1 ! hash of position in m
357 77 : indlmn(3,pspindex)=kk ! index of projector
358 77 : indlmn(4,pspindex)=ll*ll+(1-useylm)*ll+mm ! hash of position in l(l+1) array
359 77 : indlmn(5,pspindex)=iln ! absolute index of l, n disregarding m values
360 154 : indlmn(6,pspindex)=nn ! spin orbit index!!! NOT the n shell index
361 : end do
362 : end do
363 : end if
364 : end do
365 : end do
366 :
367 : ! Determine whether the atomic calculation to generate the pseudopotential
368 : ! is relativistic or not
369 :
370 : !DEBUG
371 : !write(std_out,*)' psp9in : pseudopotential generation relativity ', ps_Relativity(psxml)
372 : !write(std_out,*)' psp9in : SOC pseudopotential? (1=yes, 0 =no) '
373 : !write(std_out,*)' psp9in : irelt = ', irelt
374 : !write(ab_out,*)' psp9in : irelt = ', irelt
375 : !ENDDEBUG
376 :
377 : !Can now allocate grids, potentials and projectors
378 45 : ABI_MALLOC(rad,(mmax))
379 30 : ABI_MALLOC(vloc,(mmax))
380 60 : ABI_MALLOC(vpspll,(mmax,lnmax))
381 :
382 : !Feb 2015: shifted to Hamann grid for convenience - libpsml interpolates anyway
383 9015 : do ir=1,mmax
384 9015 : rad(ir) = amesh * (ir - 1)
385 : end do
386 : !! DEBUG
387 : ! do ir = 2, mmax
388 : ! write(std_out,'(i5,f20.12)')ir, rad(ir)
389 : ! end do
390 : !! ENDDEBUG
391 : !---
392 15 : write(message, '(a,i5,es16.6,es16.6)')' psp9in : mmax, amesh, rad(mmax) = ', mmax, amesh, rad(mmax)
393 15 : call wrtout(ab_out,message,'COLL')
394 15 : call wrtout(std_out,message,'COLL')
395 :
396 : !Check that rad grid is linear starting at zero
397 15 : amesh=rad(2)-rad(1)
398 15 : damesh=zero
399 8985 : do irad=2,mmax-1
400 8985 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
401 : end do
402 15 : if(damesh>tol8 .or. rad(1)/=zero) then
403 : write(message, '(5a)' )&
404 0 : & 'Pseudopotential input file requires linear radial mesh',ch10,&
405 0 : & 'starting at zero.',ch10,&
406 0 : & 'Action: check your pseudopotential input file.'
407 0 : ABI_ERROR(message)
408 : end if
409 :
410 : !Take care of the non-linear core corrections
411 : !----------------------------------------------------------------------------
412 : ! xcccrc : XC core correction cutoff radius (bohr)
413 : ! It is defined as the radius where the pseudo-core
414 : ! charge density becomes zero
415 : ! (here we have set up a tolerance of 1.d-12).
416 :
417 15 : rmatch = zero
418 15 : nders = 0
419 15 : maxrad = zero
420 15 : if (has_nlcc) then
421 :
422 : ! In Abinit, at least for the Troullier-Martins pseudopotential,
423 : ! the pseudocore charge density and its derivatives (xccc1d)
424 : ! are introduced in a linear grid.
425 : ! This grid is normalized, so the radial coordinates run between
426 : ! from 0 and 1 (from 0 to xcccrc, where xcccrc is the radius
427 : ! where the pseudo-core becomes zero).
428 :
429 9 : call ps_CoreCharge_get(psxml, rc=rmatch, nderivs=nders)
430 : write (message,'(1X,A,A,5X,A,1X,F8.3,A,5X,A,I8,A)') &
431 9 : & "Reading pseudocore charge",ch10, &
432 9 : & "- matching radius:",rmatch,ch10, &
433 18 : & "- number of continuous derivatives",nders,ch10
434 9 : call wrtout(std_out,message,'COLL')
435 :
436 : !Get core charge function and derivatives, if needed
437 9 : if(fchrg>1.0d-15)then
438 : ! Evaluate core charge on radial grid
439 27 : ABI_MALLOC(ff_cc,(mmax))
440 5409 : do ii=1,mmax
441 5409 : ff_cc(ii) = ps_CoreCharge_Value(psxml,rad(ii))
442 : end do
443 9 : rchrg = zero
444 9 : call psp9cc(mmax,n1xccc,rad,ff_cc,rchrg,xccc1d)
445 9 : ABI_FREE(ff_cc)
446 : ! The core charge function for pspcod=9
447 : ! becomes zero beyond rchrg. Thus xcccrc must be set
448 : ! equal to rchrg.
449 9 : xcccrc=rchrg
450 :
451 : ! Get core kinetic energy density if meta-GGA pseudopotential
452 : #if defined HAVE_LIBPSML_METAGGA
453 9 : if (has_metagga) then
454 3 : ABI_MALLOC(ff_tau,(mmax))
455 601 : do ii=1,mmax
456 601 : ff_tau(ii) = ps_CoreKineticDensity_Value(psxml,rad(ii))
457 : end do
458 1 : rchrg_tau = rchrg
459 1 : call psp9cc(mmax,n1xccc,rad,ff_tau,rchrg_tau,xcctau1d,apply_pi4i=.false.)
460 1 : ABI_FREE(ff_tau)
461 1 : write (message,'(1X,A)') "Reading pseudocore kinetic energy density for meta-GGA"
462 1 : call wrtout(std_out,message,'COLL')
463 : end if
464 : #endif
465 :
466 : else
467 0 : xccc1d(:,:) = zero
468 0 : xcctau1d(:,:) = zero
469 0 : xcccrc = zero
470 0 : fchrg = zero
471 0 : qchrg = zero
472 : end if
473 :
474 9 : maxrad = rad(mmax)
475 :
476 : end if ! has_nlcc
477 :
478 : !! DEBUG
479 : ! write(std_out,*)' xcccrc = ', xcccrc, rchrg
480 : ! write(std_out,*)
481 : ! write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc
482 : ! do ii = 1, n1xccc
483 : ! write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),&
484 : ! & xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6)
485 : ! enddo
486 : ! write(std_out,*)
487 : ! stop
488 : !! ENDDEBUG
489 :
490 :
491 : !--------------------------------------------------------------------
492 : !Carry out calculations for local (lloc) pseudopotential.
493 : !Obtain Fourier transform (1-d sine transform)
494 : !to get q^2 V(q).
495 :
496 : !Read and process vlocal:
497 : !The local potential is given by a <radfunc> element under the <local-potential>
498 : !element.
499 : !After reading, this is a copy of the treatment to the
500 : !local part carry out in psp8
501 : !i.e. (as in Hamann pseudopotential)
502 : !
503 : !Read the local component of the pseudopotential
504 9015 : vloc = zero
505 9015 : do ir = 1, mmax
506 9015 : vloc(ir) = ps_LocalPotential_Value(psxml, rad(ir))
507 : end do
508 :
509 : call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,&
510 15 : & vlspl(:,1),rad,vloc,yp1,ypn,zion)
511 :
512 : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
513 45 : ABI_MALLOC(work_spl,(mqgrid))
514 15 : call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
515 45030 : vlspl(:,2)=work_spl(:)
516 15 : ABI_FREE(work_spl)
517 :
518 : !! DEBUG
519 : ! write(std_out,*)'# Vlocal = '
520 : ! write(std_out,*)' amesh = ', amesh
521 : ! write(std_out,*)' epsatm = ', epsatm
522 : ! write(std_out,*)' mmax = ', mmax
523 : ! write(std_out,*)' mqgrid = ', mqgrid
524 : ! do ir = 1, mqgrid
525 : ! write(std_out,*)' qgrid = ', ir, qgrid(ir)
526 : ! enddo
527 : ! do ir = 1, mqgrid
528 : ! write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2)
529 : ! enddo
530 : ! write(std_out,*)
531 : ! do ir = 1, mmax
532 : ! write(std_out,*)' rad = ', rad(ir), vloc(ir)
533 : ! enddo
534 : ! write(std_out,*)
535 : ! write(std_out,*)' yp1 = ', yp1
536 : ! write(std_out,*)' ypn = ', ypn
537 : ! write(std_out,*)' zion = ', zion
538 : ! stop
539 : !! ENDDEBUG
540 :
541 :
542 : !--------------------------------------------------------------------
543 : !Take care of non-local part
544 :
545 : !Zero out all Kleinman-Bylander energies to initialize
546 92 : do ii = 1, lmnmax ! loop over all possible projectors
547 92 : if (indlmn(6,ii) == 1) then
548 77 : call ps_Projector_Get(psxml, idx_sr(indlmn(5,ii)), ekb=ekb(indlmn(5,ii)))
549 0 : else if (indlmn(6,ii) == 2) then
550 0 : call ps_Projector_Get(psxml, idx_so(indlmn(5,ii)), ekb=ekb(indlmn(5,ii)))
551 : end if
552 : end do
553 :
554 : !Read the KB projectors from the PSML file
555 : !Note than in the PSML file the radial part of the projector is stored,
556 : !while Abinit expects the radial part times the radii.
557 : !We have to multiply by r after reading it.
558 : !Note than in Hamann's format (psp8), Abinit directly reads r * radial_part_KB
559 46292 : vpspll = zero
560 92 : do ii = 1, lmnmax
561 92 : if (indlmn(6,ii) == 1) then
562 46277 : do ir = 1, mmax
563 46200 : vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_sr(indlmn(5,ii)), rad(ir))
564 46277 : vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii))
565 : end do
566 0 : else if (indlmn(6,ii) == 2) then
567 0 : do ir = 1, mmax
568 0 : vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_so(indlmn(5,ii)), rad(ir))
569 0 : vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii))
570 : end do
571 : end if
572 : end do
573 :
574 : !Allow for option of no nonlocal corrections (lloc=lmax=0)
575 15 : if (lloc==0.and.lmax==0) then
576 0 : write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
577 0 : call wrtout(ab_out,message,'COLL')
578 0 : call wrtout(std_out,message,'COLL')
579 : else
580 :
581 : ! ----------------------------------------------------------------------
582 : ! Compute Vanderbilt-KB form factors and fit splines
583 :
584 : call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,&
585 15 : & mqgrid,qgrid,rad,vpspll)
586 :
587 : end if
588 :
589 : !! DEBUG
590 : ! write(std_out,*)'# KB Projectors = '
591 : ! write(std_out,*)' amesh = ', amesh
592 : ! do ir = 1, mqgrid
593 : ! do il = 1, lnmax
594 : ! write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il)
595 : ! enddo
596 : ! enddo
597 : ! do il = 1, lmnmax
598 : ! write(std_out,*)' indlmn = ', il, indlmn(:,il)
599 : ! enddo
600 : ! write(std_out,*)' lmax = ', lmax
601 : ! write(std_out,*)' lmnmax = ', lmnmax
602 : ! write(std_out,*)' lnmax = ', lnmax
603 : ! write(std_out,*)' mmax = ', mmax
604 : ! write(std_out,*)' mqgrid = ', mqgrid
605 : ! do ir = 1, mqgrid
606 : ! write(std_out,*)' qgrid = ', ir, qgrid(ir)
607 : ! enddo
608 : ! do il = 1, lnmax
609 : ! write(std_out,*)
610 : ! write(std_out,*)'# il = ', il
611 : ! do ir = 1, mmax
612 : ! write(std_out,*)' rad = ', rad(ir), vpspll(ir,il)
613 : ! enddo
614 : ! enddo
615 : ! stop
616 : !! ENDDEBUG
617 :
618 : ! Read pseudo valence charge in real space on the linear mesh
619 : ! and transform it to reciprocal space on a regular grid. Use vloc as workspace.
620 9015 : vloc(:) = zero
621 15 : if (has_tvale) then
622 9015 : do irad=1,mmax
623 9000 : vloc(irad) = ps_ValenceCharge_Value(psxml,rad(irad))
624 9015 : vloc(irad) = vloc(irad) / four_pi
625 : end do
626 :
627 : !! DEBUG
628 : ! do irad = 1, mmax
629 : ! write(std_out,*)' Valence Charge = ', rad(irad), vloc(irad)
630 : ! enddo
631 : ! stop
632 : !! ENDDEBUG
633 :
634 :
635 : ! Check that rad grid is linear starting at zero
636 : amesh=rad(2)-rad(1)
637 15 : damesh=zero
638 8985 : do irad=2,mmax-1
639 8985 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
640 : end do
641 15 : if(damesh>tol8 .or. rad(1)/=zero) then
642 : write(message, '(5a)' )&
643 0 : & 'Pseudopotential input file requires linear radial mesh',ch10,&
644 0 : & 'starting at zero.',ch10,&
645 0 : & 'Action: check your pseudopotential input file.'
646 0 : ABI_ERROR(message)
647 : end if
648 :
649 : ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
650 15 : call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh)
651 15 : call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl)
652 15 : call pawrad_free(mesh)
653 : end if
654 :
655 : ! Read pseudo valence kinetic energy density if meta-GGA pseudopotential
656 : #if defined HAVE_LIBPSML_METAGGA
657 15 : if (has_metagga) then
658 601 : vloc(:) = zero
659 601 : do irad=1,mmax
660 601 : vloc(irad) = ps_ValenceKineticDensity_Value(psxml,rad(irad))
661 : ! NB: no 4pi division for tau quantities
662 : end do
663 :
664 1 : call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh)
665 1 : call nctab_eval_tvaletauspl(nctab, mesh, vloc, mqgrid_vl, qgrid_vl)
666 1 : call pawrad_free(mesh)
667 1 : write (message,'(1X,A)') "Reading pseudo valence kinetic energy density for meta-GGA"
668 1 : call wrtout(std_out,message,'COLL')
669 : end if
670 : #endif
671 :
672 15 : ABI_FREE(vpspll)
673 15 : ABI_FREE(vloc)
674 15 : ABI_FREE(rad)
675 15 : if (allocated(idx_sr)) then
676 15 : ABI_FREE_NOCOUNT(idx_sr)
677 : end if
678 15 : if (allocated(idx_so)) then
679 0 : ABI_FREE_NOCOUNT(idx_so)
680 : end if
681 :
682 15 : call ps_destroy(psxml)
683 :
684 : !--------------------------------------------------------------------
685 :
686 : #else
687 : ABI_UNUSED(mpsang)
688 : ABI_UNUSED(pspso)
689 : ABI_UNUSED(qgrid_vl)
690 : ABI_UNUSED(nctab%mqgrid_vl)
691 : !Initialize some arguments, for portability at compile time
692 : indlmn=0 ; mmax=0 ; nproj=0
693 : ekb=zero ; epsatm=zero ; ffspl=zero ; qchrg=zero ; vlspl=zero ; xcccrc=zero ; xccc1d=zero ; xcctau1d=zero
694 :
695 : if(.false.)write(std_out,*)filpsp ! Just to keep filpsp when HAVE_LIBPSML is false
696 : if(.false.)write(std_out,*)lloc ! Just to keep lloc when HAVE_LIBPSML is false
697 : if(.false.)write(std_out,*)lmax ! Just to keep lmax when HAVE_LIBPSML is false
698 : if(.false.)write(std_out,*)mpsang ! Just to keep mpsang when HAVE_LIBPSML is false
699 : if(.false.)write(std_out,*)pspso ! Just to keep pspso when HAVE_LIBPSML is false
700 : if(.false.)write(std_out,*)qgrid ! Just to keep qgrid when HAVE_LIBPSML is false
701 : if(.false.)write(std_out,*)qgrid_vl ! Just to keep qgrid_vl when HAVE_LIBPSML is false
702 : if(.false.)write(std_out,*)useylm ! Just to keep useylm when HAVE_LIBPSML is false
703 : if(.false.)write(std_out,*)zion ! Just to keep zion when HAVE_LIBPSML is false
704 : if(.false.)write(std_out,*)znucl ! Just to keep znucl when HAVE_LIBPSML is false
705 : #endif
706 :
707 45 : end subroutine psp9in
708 : !!***
709 :
710 : !!****f* ABINIT/psp9cc
711 : !! NAME
712 : !! psp9cc
713 : !!
714 : !! FUNCTION
715 : !! Given a radial function on a linear grid (e.g. core charge density
716 : !! or core kinetic energy density), compute its first 4 derivatives
717 : !! and interpolate onto a normalized uniform mesh for use in ABINIT.
718 : !!
719 : !! INPUTS
720 : !! mmax=maximum number of points in real space grid
721 : !! n1xccc=dimension of result1d ; 0 if no XC core correction is used
722 : !! rad(mmax)=radial grid points
723 : !! ff_values(mmax)=function values on the radial grid
724 : !! apply_pi4i=if .true. (default), multiply by 1/(4*pi); set to .false.
725 : !! for kinetic energy density quantities which have no 4pi factor
726 : !!
727 : !! OUTPUT
728 : !! rchrg=cut-off radius (bohr). If zero on input, determined from the data;
729 : !! if positive on input, used as-is (e.g. reusing core charge cutoff for tau).
730 : !! result1d(n1xccc,6)= 1D function and its four first derivatives on normalized grid
731 : !!
732 : !! SOURCE
733 :
734 10 : subroutine psp9cc(mmax,n1xccc,rad,ff_values,rchrg,result1d,apply_pi4i)
735 :
736 : !Arguments ------------------------------------
737 : !scalars
738 : integer,intent(in) :: mmax,n1xccc
739 : real(dp),intent(inout) :: rchrg
740 : logical,intent(in),optional :: apply_pi4i
741 : !arrays
742 : real(dp),intent(in) :: rad(mmax)
743 : real(dp),intent(in) :: ff_values(mmax)
744 : real(dp),intent(inout) :: result1d(n1xccc,6)
745 :
746 : !Local variables-------------------------------
747 : !scalars
748 : integer :: i1xccc,idum,irad,jj
749 : real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,normfact,tff,xp,xpm1,xpm2,xpp1,xx,twelvth
750 : logical :: do_pi4i
751 : character(len=500) :: message
752 : !arrays
753 : integer :: iwork(8)
754 : real(dp) :: rscale(5),dpoly(6,6),vpoly(6)
755 10 : real(dp),allocatable :: ff(:,:)
756 :
757 : !**********************************************************************
758 :
759 10 : do_pi4i = .true.
760 10 : if (present(apply_pi4i)) do_pi4i = apply_pi4i
761 :
762 : !Check that rad grid is linear starting at zero
763 10 : amesh=rad(2)-rad(1)
764 10 : damesh=zero
765 5990 : do irad=2,mmax-1
766 5990 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
767 : end do
768 :
769 10 : if(damesh>tol8 .or. rad(1)/=zero) then
770 : write(message, '(5a)' )&
771 0 : & 'Pseudopotential input file requires linear radial mesh',ch10,&
772 0 : & 'starting at zero.',ch10,&
773 0 : & 'Action: check your pseudopotential input file.'
774 0 : ABI_ERROR(message)
775 : end if
776 :
777 30 : ABI_MALLOC(ff,(mmax,5))
778 :
779 10 : dri = one / amesh
780 10 : pi4i = quarter / pi
781 10 : twelvth = one / 12.0_dp
782 :
783 : !Normalization factor: 1/(4*pi) for densities, 1 for kinetic energy densities
784 11 : normfact = merge(pi4i, one, do_pi4i)
785 :
786 : !Store the function values
787 30060 : ff(:,:) = zero
788 6010 : do jj=1,mmax
789 6010 : ff(jj,1) = ff_values(jj)
790 : end do
791 :
792 : !Calculate 4 first derivatives with 5-point stencil, except borders
793 5970 : do irad=3,mmax-2
794 : ff(irad,2) = (-ff(irad+2,1) + 8.0d0*ff(irad+1,1) - &
795 5960 : & 8.0d0*ff(irad-1,1) + ff(irad-2,1)) * twelvth * dri
796 : ff(irad,3) = (-ff(irad+2,1) + 16.0d0*ff(irad+1,1) - 30.0d0*ff(irad,1) + &
797 5960 : & 16.0d0*ff(irad-1,1) - ff(irad-2,1)) * twelvth * dri * dri
798 : ff(irad,4) = (ff(irad+2,1) - 2.0d0*ff(irad+1,1) + &
799 5960 : & 2.0d0*ff(irad-1,1) - ff(irad-2,1)) * half * dri * dri * dri
800 : ff(irad,5) = (ff(irad+2,1) - 4.0d0*ff(irad+1,1) + 6.0d0*ff(irad,1) - &
801 5970 : & 4.0d0*ff(irad-1,1) + ff(irad-2,1)) * dri * dri * dri * dri
802 : end do
803 :
804 : !Add border near zero using polynomial fit
805 10 : dpoly(:,:) = zero
806 70 : dpoly(:,1) = one
807 10 : vpoly(:) = zero
808 10 : vpoly(1) = ff(1,1)
809 60 : do irad=2,6
810 350 : do jj=1,6
811 350 : dpoly(irad,jj) = rad(irad)**(jj-1)
812 : end do
813 60 : vpoly(irad) = ff(irad,1)
814 : end do
815 10 : call dgesv(6,1,dpoly,6,iwork,vpoly,6,idum)
816 :
817 30 : do irad=1,2
818 : ff(irad,2) = &
819 : & vpoly(2) + 2.0d0*vpoly(3)*rad(irad) + &
820 : & 3.0d0*vpoly(4)*rad(irad)*rad(irad) + &
821 : & 4.0d0*vpoly(5)*rad(irad)*rad(irad)*rad(irad) + &
822 20 : & 5.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)*rad(irad)
823 : ff(irad,3) = &
824 : & 2.0d0*vpoly(3)*rad(irad) + &
825 : & 6.0d0*vpoly(4)*rad(irad) + &
826 : & 12.0d0*vpoly(5)*rad(irad)*rad(irad) + &
827 20 : & 20.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)
828 : ff(irad,4) = &
829 : & 6.0d0*vpoly(4) + &
830 : & 24.0d0*vpoly(5)*rad(irad) + &
831 20 : & 60.0d0*vpoly(6)*rad(irad)*rad(irad)
832 : ff(irad,5) = 24.0d0*vpoly(5) + &
833 30 : & 120.0d0*vpoly(6)*rad(irad)
834 : end do
835 :
836 : !Make linear approximation for the tail near mmax
837 30 : do irad=1,2
838 20 : ff(mmax-2+irad,2) = ff(mmax-2,2) + irad * (ff(mmax-2,2) - ff(mmax-3,2))
839 20 : ff(mmax-2+irad,3) = ff(mmax-2,3) + irad * (ff(mmax-2,3) - ff(mmax-3,3))
840 20 : ff(mmax-2+irad,4) = ff(mmax-2,4) + irad * (ff(mmax-2,4) - ff(mmax-3,4))
841 30 : ff(mmax-2+irad,5) = ff(mmax-2,5) + irad * (ff(mmax-2,5) - ff(mmax-3,5))
842 : end do
843 :
844 : !Determine rchrg where the function becomes 0, unless already provided
845 10 : if (rchrg <= zero) then
846 9 : rchrg = zero
847 1010 : do jj=mmax,1,-1
848 1010 : if (ff(jj,1) > tol13) then
849 9 : rchrg=rad(jj)
850 9 : exit
851 : end if
852 : end do
853 : end if
854 :
855 : !Check that rchrg is consistent with last grid point
856 10 : if(rchrg>rad(mmax)) then
857 : write(message, '(5a)' )&
858 0 : & 'Pseudopotential input file core charge mesh',ch10,&
859 0 : & 'is inconsistent with rchrg in header.',ch10,&
860 0 : & 'Action: check your pseudopotential input file.'
861 0 : ABI_ERROR(message)
862 : end if
863 :
864 : !Factors for unit range scaling
865 60 : do jj = 1, 5
866 60 : rscale(jj)=rchrg**(jj-1)
867 : end do
868 :
869 : !Generate uniform mesh xx in the box cut by rchrg
870 : !and interpolate the function and derivatives
871 : !Cubic polynomial interpolation is used which is consistent
872 : !with the original interpolation of these functions from
873 : !a log grid to the input linear grid.
874 :
875 25020 : dri=1.d0/amesh
876 25020 : do i1xccc=1,n1xccc
877 25010 : xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
878 :
879 : ! index to find bracketing input mesh points
880 25010 : irad = int(dri * xx) + 1
881 25010 : irad = max(irad,2)
882 25010 : irad = min(irad,mmax-2)
883 : ! interpolation coefficients
884 25010 : xp = dri * (xx - rad(irad))
885 25010 : xpp1 = xp + one
886 25010 : xpm1 = xp - one
887 25010 : xpm2 = xp - two
888 25010 : c1 = -xp * xpm1 * xpm2 * sixth
889 25010 : c2 = xpp1 * xpm1 * xpm2 * half
890 25010 : c3 = - xp * xpp1 * xpm2 * half
891 25010 : c4 = xp * xpp1 * xpm1 * sixth
892 : ! Now do the interpolation on all derivatives for this grid point
893 : ! Include normalization and unit range scaling
894 150070 : do jj=1,5
895 : tff = c1 * ff(irad - 1, jj) &
896 : & + c2 * ff(irad , jj) &
897 : & + c3 * ff(irad + 1, jj) &
898 125050 : & + c4 * ff(irad + 2, jj)
899 150060 : result1d(i1xccc,jj)=normfact*rscale(jj)*tff
900 : end do
901 : end do
902 :
903 : !5th derivative is apparently not in use, so set to zero
904 25020 : result1d(:,6)=zero
905 :
906 10 : ABI_FREE(ff)
907 :
908 10 : end subroutine psp9cc
909 : !!***
910 :
911 : end module m_psp9
912 : !!***
|
