Branch data Line data Source code
1 : : !!****m* ABINIT/m_psp8
2 : : !! NAME
3 : : !! m_psp8
4 : : !!
5 : : !! FUNCTION
6 : : !! Initialize pspcod=8 (pseudopotentials in the format generated by DRH):
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (DRH, XG, AF)
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_psp8
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_splines
28 : :
29 : : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free
30 : : use defs_datatypes, only : nctab_t
31 : : use m_psps, only : nctab_eval_tvalespl
32 : : use m_psptk, only : psp8lo, psp8nl
33 : :
34 : : implicit none
35 : :
36 : : private
37 : : !!***
38 : :
39 : : public :: psp8in
40 : : !!***
41 : :
42 : : contains
43 : : !!***
44 : :
45 : : !!****f* ABINIT/psp8in
46 : : !! NAME
47 : : !! psp8in
48 : : !!
49 : : !! FUNCTION
50 : : !! Initialize pspcod=8 (pseudopotentials in the format generated by DRH):
51 : : !! continue to read the corresponding file, then compute local and non-local potentials.
52 : : !!
53 : : !! INPUTS
54 : : !! lloc=angular momentum choice of local pseudopotential
55 : : !! lmax=value of lmax mentioned at the second line of the psp file
56 : : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
57 : : !! =if useylm=0, max number of (l,n) comp. over all type of psps
58 : : !! lnmax=max. number of (l,n) components over all type of psps
59 : : !! mmax=maximum number of points in real space grid in the psp file
60 : : !! angular momentum of nonlocal pseudopotential
61 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
62 : : !! mpssoang= Maximum number of channels, including those for treating the spin-orbit coupling
63 : : !! when mpspso=1, mpssoang=mpsang
64 : : !! when mpspso=2, mpssoang=2*mpsang-1
65 : : !! mqgrid=dimension of q (or G) grid for arrays.
66 : : !! mqgrid_vl=dimension of q (or G) grid for valence charge (array qgrid_vl)
67 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
68 : : !! qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
69 : : !! qgrid_vl(psps%mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for valence charge
70 : : !! pspso=spin-orbit characteristics, govern the content of ffspl and ekb
71 : : !! if =0: this input requires no spin-orbit characteristics of the psp
72 : : !! if =2: this input requires hgh or psp8/upf2 characteristics of the psp
73 : : !! if =3: this input requires hfn characteristics of the psp
74 : : !! useylm=governs the way the nonlocal operator is to be applied:
75 : : !! 1=using Ylm, 0=using Legendre polynomials
76 : : !! zion=nominal valence of atom as specified in psp file
77 : : !! znucl=nuclear number of atom as specified in psp file
78 : : !!
79 : : !! OUTPUT
80 : : !! ekb(lnmax)=Kleinman-Bylander energy, as read from input file
81 : : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree)
82 : : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
83 : : !! second derivative from spline fit for each angular momentum and
84 : : !! each projector
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 : : !! nproj(mpssoang)=number of projection functions for each angular momentum
88 : : !! qchrg is not used, and could be suppressed later
89 : : !! vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit
90 : : !! xcccrc=XC core correction cutoff radius (bohr)
91 : : !! xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
92 : : !! nctab<nctab_t>=NC tables
93 : : !! %has_tvale=True if the pseudo contains the pseudo valence charge
94 : : !! %tvalespl(mqgrid_vl,2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
95 : : !!
96 : : !! SOURCE
97 : :
98 : 485 : subroutine psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
99 : 485 : mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,&
100 : 485 : useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad)
101 : :
102 : : !Arguments ------------------------------------
103 : : !scalars
104 : : integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid,mqgrid_vl
105 : : integer,intent(in) :: pspso,n1xccc,useylm
106 : : real(dp),intent(in) :: zion,znucl
107 : : real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad
108 : : type(nctab_t),intent(inout) :: nctab
109 : : !arrays
110 : : integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang)
111 : : real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl)
112 : : real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
113 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
114 : :
115 : : !Local variables-------------------------------
116 : : !scalars
117 : : integer :: extension_switch,iln,iln0,pspindex,ipsang,irad,jj,kk,ll,ll_err,llin
118 : : integer :: mm,nn,nso,ir
119 : : real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn
120 : : logical :: has_tvale, debug
121 : : character(len=500) :: msg,errmsg
122 : 485 : type(pawrad_type) :: mesh
123 : : !arrays
124 : 485 : integer, allocatable :: nproj_tmp(:)
125 : 485 : real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),vpspll_irad(:),work_spl(:)
126 : :
127 : : ! ***************************************************************************
128 : :
129 : : !File format of formatted drh psp input, as adapted for use
130 : : !by the ABINIT code (the 3 first lines have already been read in calling -pspatm- routine):
131 : :
132 : : !(1) title (character) line
133 : : !(2) znucl,zion,pspdat
134 : : !(3) pspcod,pspxc,lmax,lloc,mmax,r2well (r2well not used)
135 : : !(4) rchrg,fchrg,qchrg (fchrg /=0 if core charge, qchrg not used)
136 : : !(5) nproj(0:lmax) (several projectors allowed for each l)
137 : : !(6) extension_switch(2) (spin-orbit parameters)
138 : : !Then, for ll=0,lmax :
139 : : !if(nproj(ll)>0)
140 : : !1/<u1|vbkb1>, 1/<u2|vbkb2>, ...
141 : : !for irad=1,mmax : irad, r(irad), vbkb1(irad,ll), vbkb2(irad,ll), ...
142 : : !else if ll=lloc
143 : : !for irad=1,mmax : irad, r(irad), vloc(irad)
144 : : !end if
145 : : !
146 : : !If(lloc>lmax):
147 : : ! for irad=1,mmax : irad, r(irad), vloc(irad)
148 : : !end if
149 : : !
150 : : !vbkb are Bloechl-Kleinman-Bylander projectors,(vpsp(r,ll)-vloc(r))*u(r,ll), unnormalized.
151 : : !Note that an arbitrary local potential is allowed.
152 : : !Set lloc>lmax, and provide projectors for all ll<=lmax
153 : : !
154 : : !Finally, if fchrg>0:
155 : : !
156 : : ! for irad=1,mmax : irad, r(irad), xccc(irad),
157 : : ! xccc'(irac), xccc''(irad), xccc'''(irad), xccc''''(irad)
158 : : !
159 : :
160 : 485 : debug = .False.!; debug = .True.
161 : :
162 : : ! Model core charge for nonlinear core xc correction, and 4 derivatives
163 : :
164 [ + - ]: 485 : read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
165 : 485 : write(msg, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
166 [ + + ]: 1455 : call wrtout([std_out, ab_out], msg)
167 : :
168 [ - + + - : 1455 : ABI_MALLOC(nproj_tmp, (mpssoang))
+ - + - -
+ + - ]
169 [ + + ]: 2110 : nproj_tmp = 0
170 : :
171 [ + - ]: 485 : read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(1:lmax+1)
172 : 485 : write(msg, '(a,5i6)' ) ' nproj',nproj_tmp(1:lmax+1)
173 [ + + ]: 1455 : call wrtout([std_out, ab_out], msg)
174 : :
175 : : !place holder for future implementation of additional optional header
176 : : !lines without invalidating existing psp files
177 : : !Now (12/2014) extended to include spin-orbit projectors
178 : :
179 : : ! The integer labeled "extension switch" on line 6
180 : : ! of the *.psp8 file will be set to 1 (non- or scalar-relativistic)
181 : : ! or 3 (relativistic) to signal to Abinit that the file contains the pseudo valence charge.
182 : :
183 : 485 : has_tvale = .False.
184 [ + - ]: 485 : read (tmp_unit,*, err=10, iomsg=errmsg) extension_switch
185 [ + + ]: 485 : if (any(extension_switch==[2, 3])) then
186 [ + - ]: 103 : read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(lmax+2:2*lmax+1)
187 : 103 : write(msg, '(5x,a,i6)' ) 'spin-orbit psp, extension_switch',extension_switch
188 [ + + ]: 309 : call wrtout([std_out, ab_out], msg)
189 : 103 : write(msg, '(5x,a,5i6)' ) ' nprojso',nproj_tmp(lmax+2:2*lmax+1)
190 [ + + ]: 309 : call wrtout([std_out, ab_out], msg)
191 : 103 : has_tvale = (extension_switch == 3)
192 [ + - ]: 382 : else if (any(extension_switch==[0,1])) then
193 : 382 : write(msg, '(5x,a,i6)' ) 'extension_switch',extension_switch
194 [ + + ]: 1146 : call wrtout([std_out, ab_out], msg)
195 : 382 : has_tvale = (extension_switch == 1)
196 : : else
197 : 0 : write(msg, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,&
198 : 0 : 'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit'
199 : 0 : ABI_ERROR(msg)
200 : : end if
201 : :
202 [ + + ]: 485 : if(lloc<4) then
203 [ - + ]: 1 : if (nproj_tmp(lloc+1)/=0) then
204 : : write(msg, '(a,i4,a,a,i4,5a)' )&
205 : 0 : 'Pseudopotential input file has nproj=',nproj_tmp(lloc+1),ch10,&
206 : 0 : 'for angular momentum',lloc,' which is the local potential.',ch10,&
207 : 0 : 'Should be 0 for the local potential',ch10,&
208 : 0 : 'Action: check your pseudopotential input file.'
209 : 0 : ABI_ERROR(msg)
210 : : end if
211 : : end if
212 : :
213 : : !--------------------------------------------------------------------
214 : :
215 : : !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn
216 : : ! if(pspso==2) then
217 [ + + ]: 485 : if (any(extension_switch == [0,1])) then
218 : : nso=1
219 [ + - ]: 103 : else if (any(extension_switch == [2,3])) then
220 : 103 : nso=2
221 [ + + ]: 103 : if (pspso==0) then
222 : 49 : write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,&
223 : 98 : ' but internal variable pspso is 0'
224 : 49 : ABI_COMMENT(msg)
225 : : end if
226 : : else
227 : 0 : write(msg, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,&
228 : 0 : 'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit'
229 : 0 : ABI_ERROR(msg)
230 : : end if
231 : :
232 [ + + + + ]: 34351 : pspindex=0; iln=0; indlmn=0
233 [ + + ]: 1073 : do nn=1,nso
234 [ + + ]: 2602 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
235 : 1529 : ll = ipsang-(nn-1)*lmax-1
236 [ + + ]: 2117 : if (nproj_tmp(ipsang)>0) then
237 [ + + ]: 4891 : do kk=1,nproj_tmp(ipsang)
238 : 3363 : iln = iln+1
239 [ + + ]: 9100 : do mm=1,2*ll*useylm+1
240 : 4209 : pspindex = pspindex + 1
241 : 4209 : indlmn(1,pspindex) = ll
242 : 4209 : indlmn(2,pspindex) = mm-ll*useylm-1
243 : 4209 : indlmn(3,pspindex) = kk
244 : 4209 : indlmn(4,pspindex) = ll*ll+(1-useylm)*ll+mm
245 : 4209 : indlmn(5,pspindex) = iln
246 : 7572 : indlmn(6,pspindex) = nn
247 : : !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex
248 : : end do
249 : : end do
250 : : end if
251 : : end do
252 : : end do
253 : :
254 : : ! repackage nproj_tmp for proper use by pspatm
255 [ + + ]: 2110 : nproj(:)=0
256 [ + + ]: 1805 : nproj(1:lmax+1)=nproj_tmp(1:lmax+1)
257 [ + + ]: 485 : if(pspso==2) then
258 [ + + ]: 165 : nproj(mpsang+1:mpsang+lmax)=nproj_tmp(lmax+2:2*lmax+1)
259 : : end if
260 : :
261 : : !Can now allocate grids, potentials and projectors
262 [ - + + - : 1455 : ABI_MALLOC(rad,(mmax))
+ - + - -
+ - + ]
263 [ - + + - : 970 : ABI_MALLOC(vloc,(mmax))
- + - + ]
264 [ - + + - : 1940 : ABI_MALLOC(vpspll,(mmax,lnmax))
- + + - +
- + - - +
- + ]
265 [ - + + - : 1455 : ABI_MALLOC(vpspll_irad,(lnmax))
+ - + - -
+ + - ]
266 : :
267 : : !Will now proceed at the reading of pots and projectors
268 : :
269 : : !rad(:)=radial grid r(i)
270 : : !vpspll(:,1),...,vpspll(:,lnmax)=nonlocal projectors
271 : : !vloc(:)=local potential
272 : :
273 : : !Read Vanderbilt-Kleinman-Bylander energies and projectors for each l
274 : : !or read local potential for l=lloc.
275 : : !Also get rad array (actually read more than once)
276 : 1073 : ll_err=0
277 : 1073 : iln0=0
278 [ + + ]: 1073 : do nn=1,nso
279 [ + + ]: 2117 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
280 : 1529 : ll=ipsang-(nn-1)*lmax-1
281 [ + + ]: 2117 : if (nproj_tmp(ipsang)>0) then
282 [ + - ]: 1528 : read(tmp_unit,*, err=10, iomsg=errmsg) llin,ekb(iln0+1:iln0+nproj_tmp(ipsang))
283 [ + - ]: 1528 : if(llin/=ll) then
284 : : ll_err=ipsang
285 : : exit
286 : : end if
287 [ + + ]: 882871 : do irad=1,mmax
288 [ + - ]: 881343 : read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vpspll_irad(1:nproj_tmp(ipsang))
289 [ + + ]: 2833931 : vpspll(irad,iln0+1:iln0+nproj_tmp(ipsang)) = vpspll_irad(1:nproj_tmp(ipsang))
290 : : end do
291 : 1528 : iln0=iln0+nproj_tmp(ipsang)
292 [ + - + - ]: 1 : elseif(ll==lloc .and. nn==1) then
293 [ + - ]: 1 : read(tmp_unit,*, err=10, iomsg=errmsg) llin
294 [ + - ]: 1 : if(llin/=ll) then
295 : : ll_err=ipsang
296 : : exit
297 : : end if
298 [ + + ]: 603 : do irad=1,mmax
299 [ + - ]: 603 : read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad)
300 : : end do
301 : : end if
302 : : end do !ipsang
303 : :
304 : : ! Provision for general local potential /= any angular momentum potential
305 [ + + + + ]: 1073 : if(nn==1 .and. lloc>lmax) then
306 [ + - ]: 484 : read(tmp_unit,*, err=10, iomsg=errmsg) llin
307 [ + - ]: 484 : if(llin==lloc) then
308 [ + + ]: 274997 : do irad=1,mmax
309 [ + - ]: 274997 : read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad)
310 : : end do
311 : : else
312 : 0 : ll_err=lloc+1
313 : 0 : exit
314 : : end if
315 : : end if
316 : : end do !nn
317 : :
318 [ - + ]: 485 : if(ll_err>0) then
319 : : write(msg, '(5a,i4,a,i4,a,a)' )&
320 : 0 : 'Pseudopotential input file does not have angular momenta in order',ch10,&
321 : 0 : 'or has inconsistent general local potential index',ch10,&
322 : 0 : 'Expected',ll_err-1,' , got',ll,ch10,&
323 : 0 : 'Action: check your pseudopotential input file.'
324 : 0 : ABI_ERROR(msg)
325 : : end if
326 : :
327 : : ! Check that rad grid is linear starting at zero
328 : 485 : amesh=rad(2)-rad(1)
329 : 485 : damesh=zero
330 [ + + ]: 274630 : do irad=2,mmax-1
331 : 274630 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
332 : : end do
333 [ + - - + ]: 485 : if(damesh>tol8 .or. rad(1)/=zero) then
334 : : write(msg, '(5a)' )&
335 : 0 : 'Pseudopotential input file requires linear radial mesh',ch10,&
336 : 0 : 'starting at zero.',ch10,&
337 : 0 : 'Action: check your pseudopotential input file.'
338 : 0 : ABI_ERROR(msg)
339 : : end if
340 : :
341 : : !Get core charge function and derivatives, if needed
342 [ + + ]: 485 : if(fchrg>1.0d-15)then
343 : 421 : call psp8cc(mmax, n1xccc, rchrg, xccc1d)
344 : : ! The core charge function for pspcod=8 becomes zero beyond rchrg.
345 : : ! Thus xcccrc must be set equal to rchrg.
346 : 421 : xcccrc=rchrg
347 : : else
348 [ + + + + ]: 495646 : xccc1d(:,:) = zero
349 : 64 : xcccrc = zero
350 : 64 : fchrg = zero
351 : 64 : qchrg = zero
352 : : end if
353 : :
354 : 485 : maxrad = rad(mmax)
355 : :
356 : : !! DEBUG
357 : : !write(std_out,*)' xcccrc = ', xcccrc, rchrg
358 : : !write(std_out,*)
359 : : !write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc
360 : : !do ii = 1, n1xccc
361 : : !write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),&
362 : : ! xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6)
363 : : !enddo
364 : : !write(std_out,*)
365 : : !stop
366 : : !! ENDDEBUG
367 : :
368 : :
369 : : !--------------------------------------------------------------------
370 : : !Carry out calculations for local (lloc) pseudopotential.
371 : : !Obtain Fourier transform (1-d sine transform) to get q^2 V(q).
372 : :
373 : 485 : call psp8lo(amesh, epsatm, mmax, mqgrid, qgrid, vlspl(:,1), rad, vloc, yp1, ypn, zion)
374 : :
375 : : ! Fit spline to q^2 V(q) (Numerical Recipes subroutine)
376 [ - + + - : 1455 : ABI_MALLOC(work_spl,(mqgrid))
+ - + - -
+ - + ]
377 : 485 : call spline(qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
378 [ + + ]: 1474353 : vlspl(:,2)=work_spl(:)
379 [ - + ]: 485 : ABI_FREE(work_spl)
380 : :
381 : : if (debug) then
382 : : write(std_out,*)'# Vlocal psp8 = '
383 : : write(std_out,*)' amesh = ', amesh
384 : : write(std_out,*)' epsatm = ', epsatm
385 : : write(std_out,*)' mmax = ', mmax
386 : : write(std_out,*)' mqgrid = ', mqgrid
387 : : do ir = 1, mqgrid
388 : : write(std_out,*)' qgrid = ', ir, qgrid(ir)
389 : : enddo
390 : : do ir = 1, mqgrid
391 : : write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2)
392 : : enddo
393 : : write(std_out,*)
394 : : do ir = 1, mmax
395 : : write(std_out,*)' rad = ', rad(ir), vloc(ir)
396 : : enddo
397 : : write(std_out,*)
398 : : write(std_out,*)' yp1 = ', yp1
399 : : write(std_out,*)' ypn = ', ypn
400 : : write(std_out,*)' zion = ', zion
401 : : stop
402 : : end if
403 : :
404 : :
405 : : !--------------------------------------------------------------------
406 : : !Take care of non-local part
407 : :
408 : : !Allow for option of no nonlocal corrections (lloc=lmax=0)
409 [ - + - - ]: 485 : if (lloc == 0 .and. lmax == 0) then
410 : 0 : write(msg, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
411 [ # # ]: 0 : call wrtout([std_out, ab_out], msg)
412 : : else
413 : :
414 : : ! Compute Vanderbilt-KB form factors and fit splines
415 : 485 : call psp8nl(amesh, ffspl, indlmn, lmax, lmnmax, lnmax, mmax, mqgrid, qgrid, rad, vpspll)
416 : : end if
417 : :
418 : : !! DEBUG
419 : : ! write(std_out,*)'# KB Projectors = '
420 : : ! write(std_out,*)' amesh = ', amesh
421 : : ! do ir = 1, mqgrid
422 : : ! do il = 1, lnmax
423 : : ! write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il)
424 : : ! enddo
425 : : ! enddo
426 : : ! do il = 1, lmnmax
427 : : ! write(std_out,*)' indlmn = ', il, indlmn(:,il)
428 : : ! enddo
429 : : ! write(std_out,*)' lmax = ', lmax
430 : : ! write(std_out,*)' lmnmax = ', lmnmax
431 : : ! write(std_out,*)' lnmax = ', lnmax
432 : : ! write(std_out,*)' mmax = ', mmax
433 : : ! write(std_out,*)' mqgrid = ', mqgrid
434 : : ! do ir = 1, mqgrid
435 : : ! write(std_out,*)' qgrid = ', ir, qgrid(ir)
436 : : ! enddo
437 : : ! do il = 1, lnmax
438 : : ! write(std_out,*)
439 : : ! write(std_out,*)'# il = ', il
440 : : ! do ir = 1, mmax
441 : : ! write(std_out,*)' rad = ', rad(ir), vpspll(ir,il)
442 : : ! enddo
443 : : ! enddo
444 : : ! stop
445 : : !! ENDDEBUG
446 : :
447 : : ! Read pseudo valence charge in real space on the linear mesh
448 : : ! and transform it to reciprocal space on a regular grid. Use vloc as workspace.
449 [ + + ]: 275600 : vloc(:) = zero
450 [ + + ]: 485 : if (has_tvale) then
451 [ + + ]: 272078 : do irad=1,mmax
452 [ + - ]: 271600 : read(tmp_unit,*, err=10, iomsg=errmsg)jj, rad(irad), vloc(irad)
453 : 272078 : vloc(irad) = vloc(irad) / four_pi
454 : : end do
455 : :
456 : : ! Check that rad grid is linear starting at zero
457 : 478 : amesh = rad(2) - rad(1); damesh = zero
458 [ + + ]: 271122 : do irad=2,mmax-1
459 : 271122 : damesh = max(damesh, abs(rad(irad)+amesh-rad(irad+1)))
460 : : end do
461 : :
462 [ + - - + ]: 478 : if (damesh > tol8 .or. abs(rad(1)) > tol16) then
463 : : write(msg,'(3a)')&
464 : 0 : 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
465 : 0 : 'Action: check your pseudopotential file.'
466 : 0 : ABI_ERROR(msg)
467 : : end if
468 : :
469 : : ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
470 : 478 : call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh)
471 : 478 : call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl)
472 : 478 : call pawrad_free(mesh)
473 : : end if
474 : :
475 [ - + ]: 485 : ABI_FREE(vpspll)
476 [ - + ]: 485 : ABI_FREE(vpspll_irad)
477 [ - + ]: 485 : ABI_FREE(vloc)
478 : 485 : ABI_FREE(rad)
479 [ - + ]: 485 : ABI_FREE(nproj_tmp)
480 : :
481 : 485 : return
482 : :
483 : : ! Handle IO error
484 : : 10 continue
485 : 0 : ABI_ERROR(errmsg)
486 : :
487 [ - + - + : 485 : end subroutine psp8in
- + - + -
+ - + - +
- + ]
488 : : !!***
489 : :
490 : : !!****f* m_psp8/psp8cc
491 : : !! NAME
492 : : !! psp8cc
493 : : !!
494 : : !! FUNCTION
495 : : !! Compute the core charge density, for use in the XC core
496 : : !! correction, following the function definition valid
497 : : !! for format 8 of the pseudopotentials.
498 : : !!
499 : : !! INPUTS
500 : : !! mmax=maximum number of points in real space grid in the psp file
501 : : !! n1xccc=dimension of xccc1d; 0 if no XC core correction is used
502 : : !! rchrg=cut-off radius for the core density
503 : : !!
504 : : !! OUTPUT
505 : : !! xccc1d(n1xccc,6)= 1D core charge function and its four first derivatives
506 : : !!
507 : : !! SOURCE
508 : :
509 : 421 : subroutine psp8cc(mmax, n1xccc, rchrg, xccc1d)
510 : :
511 : : !Arguments ------------------------------------
512 : : !scalars
513 : : integer,intent(in) :: mmax,n1xccc
514 : : real(dp),intent(in) :: rchrg
515 : : !arrays
516 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
517 : :
518 : : !Local variables-------------------------------
519 : : !scalars
520 : : integer :: i1xccc,idum,irad,jj
521 : : real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx
522 : : character(len=500) :: msg,errmsg
523 : : !arrays
524 : : real(dp) :: rscale(5)
525 : 421 : real(dp),allocatable :: ff(:,:),rad(:)
526 : :
527 : : !**********************************************************************
528 : :
529 [ - + + - : 1263 : ABI_MALLOC(ff,(mmax,5))
+ - + - -
+ - + ]
530 [ - + + - : 1263 : ABI_MALLOC(rad,(mmax))
+ - + - -
+ - + ]
531 : :
532 : : pi4i=quarter/pi
533 : : !
534 : : ! Read from pp file the model core charge and its first 4 derivatives
535 : : ! assumed to be on a linear grid starting at zero.
536 : : ! The input functions contain the 4pi factor, and must be rescaled.
537 : :
538 [ + + ]: 250936 : do irad=1,mmax
539 [ + - ]: 250936 : read(tmp_unit,*, err=10, iomsg=errmsg) idum,rad(irad),(ff(irad,jj),jj=1,5)
540 : : end do
541 : :
542 : : ! Check that rad grid is linear starting at zero
543 : 421 : amesh=rad(2)-rad(1)
544 : 421 : damesh=zero
545 [ + + ]: 250094 : do irad=2,mmax-1
546 : 250094 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
547 : : end do
548 : :
549 [ + - - + ]: 421 : if(damesh>tol8 .or. rad(1)/=zero) then
550 : : write(msg, '(5a)' )&
551 : 0 : 'Pseudopotential input file requires linear radial mesh',ch10,&
552 : 0 : 'starting at zero.',ch10,&
553 : 0 : 'Action: check your pseudopotential input file.'
554 : 0 : ABI_ERROR(msg)
555 : : end if
556 : :
557 : : ! Check that input rchrg is consistent with last grid point
558 [ - + ]: 421 : if(rchrg>rad(mmax)) then
559 : : write(msg, '(5a)' )&
560 : 0 : 'Pseudopotential input file core charge mesh',ch10,&
561 : 0 : 'is inconsistent with rchrg in header.',ch10,&
562 : 0 : 'Action: check your pseudopotential input file.'
563 : 0 : ABI_ERROR(msg)
564 : : end if
565 : :
566 : : !Factors for unit range scaling
567 [ + + ]: 2526 : do jj = 1, 5
568 : 2526 : rscale(jj)=rchrg**(jj-1)
569 : : end do
570 : :
571 : : !Generate uniform mesh xx in the box cut by rchrg
572 : : !and interpolate the core charge and derivatives
573 : : !Cubic polynomial interpolation is used which is consistent
574 : : !with the original interpolation of these functions from
575 : : !a log grid to the input linear grid.
576 : :
577 : 421 : dri=1.d0/amesh
578 [ + + ]: 1053342 : do i1xccc=1,n1xccc
579 : 1052921 : xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
580 : :
581 : : ! index to find bracketing input mesh points
582 : 1052921 : irad = int(dri * xx) + 1
583 : 1052921 : irad = max(irad,2)
584 : 1052921 : irad = min(irad,mmax-2)
585 : : ! interpolation coefficients
586 : 1052921 : xp = dri * (xx - rad(irad))
587 : 1052921 : xpp1 = xp + one
588 : 1052921 : xpm1 = xp - one
589 : 1052921 : xpm2 = xp - two
590 : 1052921 : c1 = -xp * xpm1 * xpm2 * sixth
591 : 1052921 : c2 = xpp1 * xpm1 * xpm2 * half
592 : 1052921 : c3 = - xp * xpp1 * xpm2 * half
593 : 1052921 : c4 = xp * xpp1 * xpm1 * sixth
594 : : ! Now do the interpolation on all derivatives for this grid point
595 : : ! Include 1/4pi normalization and unit range scaling
596 [ + + ]: 6317947 : do jj=1,5
597 : : tff = c1 * ff(irad - 1, jj) &
598 : : & + c2 * ff(irad , jj) &
599 : : & + c3 * ff(irad + 1, jj) &
600 : 5264605 : & + c4 * ff(irad + 2, jj)
601 : 6317526 : xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff
602 : : end do
603 : : end do
604 : :
605 : : !5th derivative is apparently not in use, so set to zero
606 [ + + ]: 1053342 : xccc1d(:,6)=zero
607 : :
608 [ - + ]: 421 : ABI_FREE(ff)
609 : 421 : ABI_FREE(rad)
610 : :
611 : : return
612 : :
613 : : ! Handle IO error
614 : : 10 continue
615 : 0 : ABI_ERROR(errmsg)
616 : :
617 [ # # # # ]: 0 : end subroutine psp8cc
618 : : !!***
619 : :
620 : : end module m_psp8
621 : : !!***
|
