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