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 479 : subroutine psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
99 479 : mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,&
100 479 : 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 479 : type(pawrad_type) :: mesh
123 : !arrays
124 479 : integer, allocatable :: nproj_tmp(:)
125 479 : 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 479 : debug = .False.!; debug = .True.
161 :
162 : ! Model core charge for nonlinear core xc correction, and 4 derivatives
163 :
164 479 : read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
165 479 : write(msg, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
166 1437 : call wrtout([std_out, ab_out], msg)
167 :
168 1437 : ABI_MALLOC(nproj_tmp, (mpssoang))
169 2080 : nproj_tmp = 0
170 :
171 479 : read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(1:lmax+1)
172 479 : write(msg, '(a,5i6)' ) ' nproj',nproj_tmp(1:lmax+1)
173 1437 : 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 479 : has_tvale = .False.
184 479 : read (tmp_unit,*, err=10, iomsg=errmsg) extension_switch
185 479 : if (any(extension_switch==[2, 3])) then
186 100 : read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(lmax+2:2*lmax+1)
187 100 : write(msg, '(5x,a,i6)' ) 'spin-orbit psp, extension_switch',extension_switch
188 300 : call wrtout([std_out, ab_out], msg)
189 100 : write(msg, '(5x,a,5i6)' ) ' nprojso',nproj_tmp(lmax+2:2*lmax+1)
190 300 : call wrtout([std_out, ab_out], msg)
191 100 : has_tvale = (extension_switch == 3)
192 379 : else if (any(extension_switch==[0,1])) then
193 379 : write(msg, '(5x,a,i6)' ) 'extension_switch',extension_switch
194 1137 : call wrtout([std_out, ab_out], msg)
195 379 : 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 479 : 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 479 : if (any(extension_switch == [0,1])) then
218 : nso=1
219 100 : else if (any(extension_switch == [2,3])) then
220 100 : nso=2
221 100 : if (pspso==0) then
222 48 : write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,&
223 96 : ' but internal variable pspso is 0'
224 48 : 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 33820 : pspindex=0; iln=0; indlmn=0
233 1058 : do nn=1,nso
234 2563 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
235 1505 : ll = ipsang-(nn-1)*lmax-1
236 2084 : if (nproj_tmp(ipsang)>0) then
237 4804 : do kk=1,nproj_tmp(ipsang)
238 3300 : iln = iln+1
239 8938 : do mm=1,2*ll*useylm+1
240 4134 : pspindex = pspindex + 1
241 4134 : indlmn(1,pspindex) = ll
242 4134 : indlmn(2,pspindex) = mm-ll*useylm-1
243 4134 : indlmn(3,pspindex) = kk
244 4134 : indlmn(4,pspindex) = ll*ll+(1-useylm)*ll+mm
245 4134 : indlmn(5,pspindex) = iln
246 7434 : 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 2080 : nproj(:)=0
256 1781 : nproj(1:lmax+1)=nproj_tmp(1:lmax+1)
257 479 : if(pspso==2) then
258 159 : 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 1437 : ABI_MALLOC(rad,(mmax))
263 958 : ABI_MALLOC(vloc,(mmax))
264 1916 : ABI_MALLOC(vpspll,(mmax,lnmax))
265 1437 : 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 1058 : ll_err=0
277 1058 : iln0=0
278 1058 : do nn=1,nso
279 2084 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
280 1505 : ll=ipsang-(nn-1)*lmax-1
281 2084 : if (nproj_tmp(ipsang)>0) then
282 1504 : read(tmp_unit,*, err=10, iomsg=errmsg) llin,ekb(iln0+1:iln0+nproj_tmp(ipsang))
283 1504 : if(llin/=ll) then
284 : ll_err=ipsang
285 : exit
286 : end if
287 868447 : do irad=1,mmax
288 866943 : read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vpspll_irad(1:nproj_tmp(ipsang))
289 2781707 : vpspll(irad,iln0+1:iln0+nproj_tmp(ipsang)) = vpspll_irad(1:nproj_tmp(ipsang))
290 : end do
291 1504 : 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 1058 : if(nn==1 .and. lloc>lmax) then
306 478 : read(tmp_unit,*, err=10, iomsg=errmsg) llin
307 478 : if(llin==lloc) then
308 271391 : do irad=1,mmax
309 271391 : 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 479 : 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 479 : amesh=rad(2)-rad(1)
329 479 : damesh=zero
330 271036 : do irad=2,mmax-1
331 271036 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
332 : end do
333 479 : 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 479 : if(fchrg>1.0d-15)then
343 415 : 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 415 : 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 479 : 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 479 : 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 1437 : ABI_MALLOC(work_spl,(mqgrid))
377 479 : call spline(qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
378 1456341 : vlspl(:,2)=work_spl(:)
379 479 : 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 479 : 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 479 : 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 271994 : vloc(:) = zero
450 479 : if (has_tvale) then
451 268472 : do irad=1,mmax
452 268000 : read(tmp_unit,*, err=10, iomsg=errmsg)jj, rad(irad), vloc(irad)
453 268472 : vloc(irad) = vloc(irad) / four_pi
454 : end do
455 :
456 : ! Check that rad grid is linear starting at zero
457 472 : amesh = rad(2) - rad(1); damesh = zero
458 267528 : do irad=2,mmax-1
459 267528 : damesh = max(damesh, abs(rad(irad)+amesh-rad(irad+1)))
460 : end do
461 :
462 472 : 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 472 : call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh)
471 472 : call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl)
472 472 : call pawrad_free(mesh)
473 : end if
474 :
475 479 : ABI_FREE(vpspll)
476 479 : ABI_FREE(vpspll_irad)
477 479 : ABI_FREE(vloc)
478 479 : ABI_FREE(rad)
479 479 : ABI_FREE(nproj_tmp)
480 :
481 479 : return
482 :
483 : ! Handle IO error
484 : 10 continue
485 0 : ABI_ERROR(errmsg)
486 :
487 479 : 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 415 : 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 415 : real(dp),allocatable :: ff(:,:),rad(:)
526 :
527 : !**********************************************************************
528 :
529 1245 : ABI_MALLOC(ff,(mmax,5))
530 1245 : 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 247330 : do irad=1,mmax
539 247330 : 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 415 : amesh=rad(2)-rad(1)
544 415 : damesh=zero
545 246500 : do irad=2,mmax-1
546 246500 : damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
547 : end do
548 :
549 415 : 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 415 : 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 2490 : do jj = 1, 5
568 2490 : 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 415 : dri=1.d0/amesh
578 1038330 : do i1xccc=1,n1xccc
579 1037915 : xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
580 :
581 : ! index to find bracketing input mesh points
582 1037915 : irad = int(dri * xx) + 1
583 1037915 : irad = max(irad,2)
584 1037915 : irad = min(irad,mmax-2)
585 : ! interpolation coefficients
586 1037915 : xp = dri * (xx - rad(irad))
587 1037915 : xpp1 = xp + one
588 1037915 : xpm1 = xp - one
589 1037915 : xpm2 = xp - two
590 1037915 : c1 = -xp * xpm1 * xpm2 * sixth
591 1037915 : c2 = xpp1 * xpm1 * xpm2 * half
592 1037915 : c3 = - xp * xpp1 * xpm2 * half
593 1037915 : 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 6227905 : do jj=1,5
597 : tff = c1 * ff(irad - 1, jj) &
598 : & + c2 * ff(irad , jj) &
599 : & + c3 * ff(irad + 1, jj) &
600 5189575 : & + c4 * ff(irad + 2, jj)
601 6227490 : 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 1038330 : xccc1d(:,6)=zero
607 :
608 415 : ABI_FREE(ff)
609 415 : 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 : !!***
|
