Line data Source code
1 : !!****m* ABINIT/m_psp_hgh
2 : !! NAME
3 : !! m_psp_hgh
4 : !!
5 : !! FUNCTION
6 : !! Initialize pspcod=2, 3, 10 pseudopotentials (GTH)
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MT, FD, PT)
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_psp_hgh
23 :
24 : use defs_basis
25 : use m_splines
26 : use m_abicore
27 : use m_errors
28 : use m_dtset
29 :
30 : use defs_datatypes, only : pseudopotential_type
31 : use m_special_funcs, only : abi_derfc
32 : #if defined HAVE_BIGDFT
33 : use BigDFT_API, only: atomic_info
34 : #endif
35 : use m_wvl_descr_psp, only : wvl_descr_psp_fill
36 :
37 : implicit none
38 :
39 : private
40 : !!***
41 :
42 : public :: psp2in
43 : public :: psp3in
44 : public :: psp10in
45 : !!***
46 :
47 : contains
48 : !!***
49 :
50 : !!****f* m_psp_hgh/psp2in
51 : !! NAME
52 : !! psp2in
53 : !!
54 : !! FUNCTION
55 : !! Initialize pspcod=2 pseudopotentials (GTH format):
56 : !! continue to read the file, then compute the corresponding
57 : !! local and non-local potentials.
58 : !!
59 : !! INPUTS
60 : !! dtset <type(dataset_type)>=all input variables in this dataset
61 : !! ipsp=id in the array of the pseudo-potential.
62 : !! lmax=value of lmax mentioned at the second line of the psp file
63 : !! zion=nominal valence of atom as specified in psp file
64 : !!
65 : !! OUTPUT
66 : !! ekb(lnmax)=Kleinman-Bylander energy,
67 : !! {{\ \begin{equation}
68 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
69 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
70 : !! \end{equation} }}
71 : !! for each (l,n)
72 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+Zv/r) dr]$ (hartree)
73 : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
74 : !! second derivative from spline fit for each angular momentum and
75 : !! each projector
76 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
77 : !! or i=lmn (if useylm=1)
78 : !! nproj(mpsang)=number of projection functions for each angular momentum
79 : !! vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
80 : !! dvlspl(mqgrid_vl,2)=dVloc(r)/dr and second derivatives from spline fit (only
81 : !! allocated if vlspl_recipSpace is false.
82 : !!
83 : !! SIDE EFFECTS
84 : !! Input/output
85 : !! lmax : at input =value of lmax mentioned at the second line of the psp file
86 : !! at output= 1
87 : !! psps <type(pseudopotential_type)>=at output, values depending on the read
88 : !! pseudo are set.
89 : !! | lmnmax(IN)=if useylm=1, max number of (l,m,n) comp. over all type of psps
90 : !! | =if useylm=0, max number of (l,n) comp. over all type of psps
91 : !! | lnmax(IN)=max. number of (l,n) components over all type of psps
92 : !! | angular momentum of nonlocal pseudopotential
93 : !! | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
94 : !! | mqgrid_ff(IN)=dimension of q (or G) grid for nl form factors (array ffspl)
95 : !! | mqgrid_vl(IN)=dimension of q (or G) grid or r grid (if vlspl_recipSpace = .false.)
96 : !! | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
97 : !! | qgrid_vl(mqgrid_vl)(IN)=values of q on grid from 0 to qmax (bohr^-1) for Vloc
98 : !! | if vlspl_recipSpace is .true. else values of r on grid from
99 : !! | 0 to 2pi / qmax * mqgrid_ff (bohr).
100 : !! | useylm(IN)=governs the way the nonlocal operator is to be applied:
101 : !! | 1=using Ylm, 0=using Legendre polynomials
102 : !! | vlspl_recipSpace(IN)=.true. if pseudo are expressed in reciprocal space.
103 : !! | gth_params(OUT)=store GTH coefficients and parameters.
104 : !!
105 : !! SOURCE
106 :
107 328 : subroutine psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion)
108 :
109 : !Arguments ------------------------------------
110 : !scalars
111 : integer,intent(in) :: ipsp,lmax
112 : real(dp),intent(in) :: zion
113 : real(dp),intent(out) :: epsatm
114 : type(dataset_type),intent(in) :: dtset
115 : type(pseudopotential_type),intent(inout) :: psps
116 : !arrays
117 : integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpsang)
118 : real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2),ekb(psps%lnmax)
119 : real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i
120 : real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2)
121 :
122 : !Local variables-------------------------------
123 : !scalars
124 : integer :: iln,index,ipsang,kk,ll,mm
125 : real(dp) :: cc1,cc2,cc3,cc4,h1p,h1s,h2s,rloc,rrp,rrs
126 : real(dp) :: yp1,ypn
127 : character(len=500) :: msg,errmsg
128 : !arrays
129 328 : real(dp),allocatable :: work_space(:),work_spl(:)
130 328 : real(dp),allocatable :: dvloc(:)
131 :
132 : ! ***************************************************************************
133 :
134 : !Set various terms to 0 in case not defined below
135 : !GTH values
136 328 : rloc=0.d0
137 328 : cc1=0.d0
138 328 : cc2=0.d0
139 328 : cc3=0.d0
140 328 : cc4=0.d0
141 328 : rrs=0.d0
142 328 : h1s=0.d0
143 328 : h2s=0.d0
144 328 : rrp=0.d0
145 328 : h1p=0.d0
146 771 : nproj(1:psps%mpsang)=0
147 :
148 : !Read and write different lines of the pseudopotential file
149 328 : read (tmp_unit,*, err=10, iomsg=errmsg) rloc,cc1,cc2,cc3,cc4
150 328 : write(msg, '(a,f12.7)' ) ' rloc=',rloc
151 328 : call wrtout(ab_out,msg,'COLL')
152 328 : call wrtout(std_out, msg,'COLL')
153 328 : write(msg, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' cc1=',cc1,'; cc2=',cc2,'; cc3=',cc3,'; cc4=',cc4
154 328 : call wrtout(ab_out,msg,'COLL')
155 328 : call wrtout(std_out, msg,'COLL')
156 :
157 328 : read (tmp_unit,*, err=10, iomsg=errmsg) rrs,h1s,h2s
158 328 : write(msg, '(a,f12.7,a,f12.7,a,f12.7)' )' rrs=',rrs,'; h1s=',h1s,'; h2s=',h2s
159 328 : call wrtout(ab_out,msg,'COLL')
160 328 : call wrtout(std_out, msg,'COLL')
161 :
162 328 : read (tmp_unit,*, err=10, iomsg=errmsg) rrp,h1p
163 328 : write(msg, '(a,f12.7,a,f12.7)' )' rrp=',rrp,'; h1p=',h1p
164 328 : call wrtout(ab_out,msg,'COLL')
165 328 : call wrtout(std_out, msg,'COLL')
166 :
167 : !Store the coefficients.
168 328 : psps%gth_params%set(ipsp) = .true.
169 2624 : psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, 0.d0, 0.d0 /)
170 2624 : psps%gth_params%psppar(1, :, ipsp) = (/ rrs, h1s, h2s, 0.d0, 0.d0, 0.d0, 0.d0 /)
171 2624 : psps%gth_params%psppar(2, :, ipsp) = (/ rrp, h1p, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /)
172 328 : if (dtset%usewvl == 1) then
173 0 : call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
174 : end if
175 :
176 328 : if (abs(h1s)>1.d-08) nproj(1)=1
177 328 : if (abs(h2s)>1.d-08) nproj(1)=2
178 :
179 328 : if (abs(h1p)>1.d-08) then
180 104 : if(psps%mpsang<2)then
181 : write(msg, '(a,es12.4,a,a,a,i2,a)' )&
182 0 : & 'With non-zero h1p (=',h1p,'), mpsang should be at least 2,',ch10,&
183 0 : & 'while mpsang=',psps%mpsang,'.'
184 0 : ABI_ERROR(msg)
185 : end if
186 104 : nproj(2)=1
187 104 : if (lmax<1) then
188 : write(msg, '(a,i5,a,e12.4,a,a,a,a)' )&
189 0 : & 'Input lmax=',lmax,' disagree with input h1p=',h1p,'.',&
190 0 : & 'Your pseudopotential is incoherent.',ch10,&
191 0 : & 'Action: correct your pseudopotential file.'
192 0 : ABI_ERROR(msg)
193 : end if
194 : end if
195 :
196 : !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn
197 4178 : index=0;iln=0;indlmn(:,:)=0
198 760 : do ipsang=1,lmax+1
199 760 : if(nproj(ipsang)>0)then
200 209 : ll=ipsang-1
201 522 : do kk=1,nproj(ipsang)
202 313 : iln=iln+1
203 835 : do mm=1,2*ll*psps%useylm+1
204 313 : index=index+1
205 313 : indlmn(1,index)=ll
206 313 : indlmn(2,index)=mm-ll*psps%useylm-1
207 313 : indlmn(3,index)=kk
208 313 : indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
209 313 : indlmn(5,index)=iln
210 626 : indlmn(6,index)=1
211 : end do
212 : end do
213 : end if
214 : end do
215 :
216 : !First, the local potential --
217 : !compute q^2V(q) or V(r)
218 : !MJV NOTE: psp2lo should never be called with dvspl unallocated, which
219 : !is possible unless .not.psps%vlspl_recipSpace
220 984 : ABI_MALLOC(dvloc,(psps%mqgrid_vl))
221 : call psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,psps%mqgrid_vl,psps%qgrid_vl,&
222 328 : & vlspl(:,1),rloc,psps%vlspl_recipSpace,yp1,ypn,zion)
223 :
224 : !Fit spline to (q^2)V(q) or V(r)
225 984 : ABI_MALLOC(work_space,(psps%mqgrid_vl))
226 656 : ABI_MALLOC(work_spl,(psps%mqgrid_vl))
227 328 : call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl)
228 1078984 : vlspl(:,2)=work_spl(:)
229 328 : if (.not.psps%vlspl_recipSpace) then
230 0 : dvlspl(:,1) = dvloc
231 0 : call spline (psps%qgrid_vl,dvlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl)
232 0 : dvlspl(:,2)=work_spl(:)
233 : end if
234 :
235 328 : ABI_FREE(work_space)
236 328 : ABI_FREE(work_spl)
237 328 : ABI_FREE(dvloc)
238 :
239 :
240 : !Second, compute KB energies and form factors and fit splines
241 878 : ekb(:)=0.0d0
242 : !First check if any nonlocal projectors are being used
243 760 : if (maxval(nproj(1:lmax+1))>0) then
244 105 : call psp2nl(ekb,ffspl,h1p,h1s,h2s,psps%lnmax,psps%mqgrid_ff,psps%qgrid_ff,rrp,rrs)
245 : end if
246 :
247 : return
248 :
249 : ! Handle IO error
250 : 10 continue
251 0 : ABI_ERROR(errmsg)
252 :
253 : end subroutine psp2in
254 : !!***
255 :
256 : !!****f* m_psp_hgh/psp2nl
257 : !! NAME
258 : !! psp2nl
259 : !!
260 : !! FUNCTION
261 : !! Goedecker-Teter-Hutter nonlocal pseudopotential (from preprint of 1996).
262 : !! Uses Gaussians for fully nonlocal form, analytic expressions.
263 : !!
264 : !! INPUTS
265 : !! h1p=factor defining strength of 1st projector for l=1 channel
266 : !! h1s=factor defining strength of 1st projector for l=0 channel
267 : !! h2s=factor defining strength of 2nd projector for l=0 channel
268 : !! lnmax=max. number of (l,n) components over all type of psps
269 : !! mqgrid=number of grid points for qgrid
270 : !! qgrid(mqgrid)=array of |G| values
271 : !! rrp=core radius for p channel (bohr)
272 : !! rrs=core radius for s channel (bohr)
273 : !!
274 : !! OUTPUT
275 : !! ekb(lnmax)=Kleinman-Bylander energy
276 : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
277 : !! second derivative from spline fit for each angular momentum
278 : !! and each projector
279 : !!
280 : !! SOURCE
281 :
282 105 : subroutine psp2nl(ekb,ffspl,h1p,h1s,h2s,lnmax,mqgrid,qgrid,rrp,rrs)
283 :
284 : !Arguments ------------------------------------
285 : !scalars
286 : integer,intent(in) :: lnmax,mqgrid
287 : real(dp),intent(in) :: h1p,h1s,h2s,rrp,rrs
288 : !arrays
289 : real(dp),intent(in) :: qgrid(mqgrid)
290 : real(dp),intent(inout) :: ekb(lnmax),ffspl(mqgrid,2,lnmax) !vz_i
291 :
292 : !Local variables-------------------------------
293 : !scalars
294 : integer :: iln,iqgrid
295 : real(dp) :: qmax,yp1,ypn
296 : !arrays
297 105 : real(dp),allocatable :: work(:)
298 :
299 : ! *************************************************************************
300 :
301 315 : ABI_MALLOC(work,(mqgrid))
302 :
303 : !Kleinman-Bylander energies ekb were set to zero in calling program
304 :
305 : !Compute KB energies
306 105 : iln=0
307 105 : if (abs(h1s)>1.d-12) then
308 105 : iln=iln+1
309 105 : ekb(iln)=h1s*32.d0*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2)
310 : end if
311 105 : if (abs(h2s)>1.d-12) then
312 104 : iln=iln+1
313 104 : ekb(iln) =h2s*(128.d0/15.d0)*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2)
314 : end if
315 105 : if (abs(h1p)>1.d-12) then
316 104 : iln=iln+1
317 104 : ekb(iln)=h1p*(64.d0/3.d0)*rrp**5*(pi**(2.5d0)/(4.d0*pi)**2)
318 : end if
319 :
320 : !Compute KB form factor
321 105 : iln=0
322 :
323 : !l=0 first projector
324 105 : if (abs(h1s)>1.d-12) then
325 316974 : iln=iln+1
326 316974 : do iqgrid=1,mqgrid
327 316974 : ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
328 : end do
329 : ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
330 105 : yp1=0.d0
331 105 : qmax=qgrid(mqgrid)
332 105 : ypn=-4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2)
333 : ! Fit spline to get second derivatives by spline fit
334 105 : call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
335 : ! else
336 : ! or else put first projector nonlocal correction at l=0 to 0
337 : ! ffspl(:,:,iln)=0.0d0
338 : end if
339 :
340 : !l=0 second projector
341 105 : if (abs(h2s)>1.d-12) then
342 104 : iln=iln+1
343 313972 : do iqgrid=1,mqgrid
344 : ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
345 313972 : & (3.d0-(two_pi*qgrid(iqgrid)*rrs)**2)
346 : end do
347 : ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
348 104 : yp1=0.d0
349 104 : qmax=qgrid(mqgrid)
350 : ypn=4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
351 104 : & (-5.d0+(two_pi*qmax*rrs)**2)
352 : ! Fit spline to get second derivatives by spline fit
353 104 : call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
354 : ! else if(mproj>=2)then
355 : ! or else put second projector nonlocal correction at l=0 to 0
356 : ! ffspl(:,:,iln)=0.0d0
357 : end if
358 :
359 : !l=1 first projector
360 105 : if (abs(h1p)>1.d-12) then
361 104 : iln=iln+1
362 313972 : do iqgrid=1,mqgrid
363 : ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
364 313972 : & (two_pi*qgrid(iqgrid))
365 : end do
366 : ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
367 104 : yp1=two_pi
368 104 : qmax=qgrid(mqgrid)
369 104 : ypn=-two_pi*((two_pi*qmax*rrp)**2-1.d0) * exp(-0.5d0*(two_pi*qmax*rrp)**2)
370 : ! Fit spline to get second derivatives by spline fit
371 104 : call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
372 : ! else if(mpsang>=2)then
373 : ! or else put first projector l=1 nonlocal correction to 0
374 : ! ffspl(:,:,iln)=0.0d0
375 : end if
376 :
377 105 : ABI_FREE(work)
378 :
379 105 : end subroutine psp2nl
380 : !!***
381 :
382 : !!****f* ABINIT/psp2lo
383 : !! NAME
384 : !! psp2lo
385 : !!
386 : !! FUNCTION
387 : !! Treat local part of Goedecker-Teter-Hutter pseudopotentials (pspcod=2),
388 : !! as well as Hartwigsen-Goedecker-Hutter pseudopotentials (pspcod=3)
389 : !!
390 : !! INPUTS
391 : !! cc1,2,3,4=parameters from analytic pseudopotential form
392 : !! mqgrid=number of grid points in q from 0 to qmax.
393 : !! qgrid(mqgrid)=values of q (or G) on grid from 0 to qmax (bohr^-1)
394 : !! if vlspl_recipSpace is .true. else values of r on grid from
395 : !! 0 to 2pi / qmax * mqgrid_ff (bohr).
396 : !! rloc=local pseudopotential core radius (bohr)
397 : !! vlspl_recipSpace= .true. if computation of vlspl is done in reciprocal space
398 : !! zion=valence charge of atom
399 : !! parameters for local potential: rloc,c1,c2,c3,c4
400 : !!
401 : !! OUTPUT
402 : !! dvloc(mqgrid)=dVloc(r)/dr (only allocated if vlspl_recipSpace is false).
403 : !! epsatm=$4\pi\int[r^2 (v(r)+\frac{Zv}{r} dr]$
404 : !!{{\ \begin{eqnarray}
405 : !! q2vq(mqgrid)&=&q^2 v(q) \nonumber \\
406 : !! &=&-Zv/\pi
407 : !! +q^2 4\pi\int[(\frac{\sin(2\pi qr)}{2\pi qr})(r^2 v(r)+r Zv)dr]\nonumber\\
408 : !! &=&\exp(-K^2*rloc^2/2) \nonumber \\
409 : !! && *(-\frac{zion}{\pi}+(\frac{K^2*rloc^3}{\sqrt{2*\pi}}*
410 : !! (c1+c2*(3-(rloc*K)^2) \nonumber \\
411 : !! && +c3*(15-10(rloc*K)^2+(rloc*K)^4) \nonumber \\
412 : !! && +c4*(105-105*(rloc*K)^2+21*(rloc*K)^4-(rloc*K)^6)) \nonumber
413 : !!\end{eqnarray} }}
414 : !! for GTH vloc with $K=(2\pi q)$.
415 : !! yp1,ypn=derivative of q^2 v(q) wrt q at q=0 and q=qmax
416 : !! (needed for spline fitter).
417 : !!
418 : !! SOURCE
419 :
420 658 : subroutine psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,mqgrid,qgrid,q2vq,&
421 : & rloc,vlspl_recipSpace,yp1,ypn,zion)
422 :
423 : !Arguments ------------------------------------
424 : !scalars
425 : integer,intent(in) :: mqgrid
426 : real(dp),intent(in) :: cc1,cc2,cc3,cc4,rloc,zion
427 : real(dp),intent(out) :: epsatm,yp1,ypn
428 : logical,intent(in) :: vlspl_recipSpace
429 : !arrays
430 : real(dp),intent(in) :: qgrid(mqgrid)
431 : real(dp),intent(out) :: dvloc(mqgrid),q2vq(mqgrid)
432 :
433 : !Local variables-------------------------------
434 : !scalars
435 : integer :: iqgrid
436 : real(dp) :: erfValue,gaussValue,polyValue,qmax,rq,rq2
437 : character(len=500) :: msg
438 :
439 : ! *************************************************************************
440 :
441 : !Compute epsatm = lim(q->0) [Vloc(q) + zion/(Pi*q^2)]
442 : epsatm=2.d0*pi*rloc**2*zion+(2.d0*pi)**(1.5d0)*rloc**3*&
443 658 : & (cc1+3.d0*cc2+15.d0*cc3+105.d0*cc4)
444 :
445 : !If vlspl_recipSpace is .true., we compute V(q)*q^2 in reciprocal space,
446 : !else we compute V(r) in real space.
447 658 : if (vlspl_recipSpace) then
448 658 : write(msg, '(a)' ) '- Local part computed in reciprocal space.'
449 658 : call wrtout(ab_out,msg,'COLL')
450 658 : call wrtout(std_out, msg,'COLL')
451 :
452 : ! d(q^2*V(q))/d(q) at q=0 and q=qmax
453 658 : qmax=qgrid(mqgrid)
454 658 : rq2=(2.d0*pi*qmax*rloc)**2
455 658 : yp1=0.d0
456 : ypn= (2.d0*pi*qmax*rloc**2)*exp(-0.5d0*rq2)* &
457 : & (2.d0*zion + sqrt(2.d0*pi)*rloc*&
458 : & (cc1*(2.d0-rq2) + cc2*(6.d0-7.d0*rq2+rq2**2) +&
459 : & cc3*(30.d0-55.d0*rq2+16.d0*rq2**2-rq2**3) +&
460 658 : & cc4*(210.d0-525.d0*rq2+231.d0*rq2**2-29.d0*rq2**3+rq2**4)))
461 : ! ypn has been tested against Maple-derived expression.
462 :
463 : ! Compute q^2*vloc(q) on uniform grid
464 2073916 : do iqgrid=1,mqgrid
465 2073258 : rq2=(2.d0*pi*qgrid(iqgrid)*rloc)**2
466 : q2vq(iqgrid)=exp(-0.5d0*rq2)*(-zion/pi+rq2*(rloc/sqrt(2.d0*pi)) *&
467 : & ( cc1 + cc2*(3.d0-rq2) + cc3*(15.d0-10.d0*rq2+rq2**2) +&
468 2073916 : & cc4*(105.d0-rq2*(105.d0-rq2*(21.d0-rq2))) ))
469 : end do
470 : else
471 0 : write(msg, '(a)' ) '- Local part computed in real space.'
472 0 : call wrtout(ab_out,msg,'COLL')
473 0 : call wrtout(std_out, msg,'COLL')
474 :
475 : ! Compute derivatives for splines computations
476 0 : yp1 = 0.d0
477 0 : rq2 = (qgrid(mqgrid) / rloc) ** 2
478 0 : erfValue = abi_derfc(sqrt(0.5d0 * rq2))
479 0 : ypn = - 2.0d0 * zion / sqrt(2.d0 * pi) / qgrid(mqgrid) / rloc
480 0 : ypn = ypn - rq2 * (cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3) / qgrid(mqgrid)
481 0 : ypn = ypn + (2.d0 * cc2 * rq2 + 4.d0 * cc3 * rq2 ** 2 + 6.d0 * cc4 * rq2 ** 3) / qgrid(mqgrid)
482 0 : ypn = ypn * exp(-0.5d0 * rq2)
483 0 : ypn = ypn + zion / qgrid(mqgrid) ** 2 * erfValue
484 : ! Note that ypn has been calculated on a full-proof a4 paper sheet.
485 :
486 : ! Compute local potential and its first derivatives.
487 0 : do iqgrid = 1, mqgrid, 1
488 0 : rq2 = (qgrid(iqgrid) / rloc) ** 2
489 : ! Compute erf() part
490 : ! Case r = 0
491 0 : gaussValue = exp(-0.5d0 * rq2)
492 0 : if (qgrid(iqgrid) == 0.d0) then
493 0 : q2vq(iqgrid) = -zion / rloc * sqrt(2.d0 / pi)
494 0 : dvloc(iqgrid) = 0.d0
495 : else
496 0 : erfValue = abi_derfc(sqrt(0.5d0 * rq2))
497 0 : q2vq(iqgrid) = -zion / qgrid(iqgrid) * (1.0d0 - erfValue)
498 : dvloc(iqgrid) = - sqrt(2.d0 / pi) * zion * gaussValue / (qgrid(iqgrid) * rloc) - &
499 0 : & q2vq(iqgrid) / qgrid(iqgrid)
500 : end if
501 : ! Add the gaussian part
502 0 : polyValue = cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3
503 0 : q2vq(iqgrid) = q2vq(iqgrid) + gaussValue * polyValue
504 0 : rq = qgrid(iqgrid) / rloc
505 : dvloc(iqgrid) = dvloc(iqgrid) - qgrid(iqgrid) / rloc ** 2 * gaussValue * polyValue + &
506 : & gaussValue * (2.0d0 * cc2 * rq / rloc + 3.0d0 * cc3 * rq ** 3 / rloc + &
507 0 : & 6.0d0 * cc4 * rq ** 5 / rloc)
508 : end do
509 :
510 : write(msg, '(a,f12.7,a,a,f12.7,a,a,a,f12.7)' ) &
511 0 : & ' | dr spline step is : ', qgrid(2), ch10, &
512 0 : & ' | r > ', qgrid(mqgrid) ,' is set to 0.', ch10, &
513 0 : & ' | last non-nul potential value is : ', q2vq(mqgrid)
514 0 : call wrtout(ab_out,msg,'COLL')
515 0 : call wrtout(std_out, msg,'COLL')
516 : end if
517 :
518 658 : end subroutine psp2lo
519 : !!***
520 :
521 : !!****f* ABINIT/psp3in
522 : !! NAME
523 : !! psp3in
524 : !!
525 : !! FUNCTION
526 : !! Initialize pspcod=3 pseudopotentials (HGH psps PRB58,3641(1998) [[cite:Hartwigsen1998]]):
527 : !! continue to read the file, then compute the corresponding
528 : !! local and non-local potentials.
529 : !!
530 : !! INPUTS
531 : !! dtset <type(dataset_type)>=all input variables in this dataset
532 : !! pspso=spin-orbit characteristics, govern the content of ffspl and ekb
533 : !! if =0 : this input requires NO spin-orbit characteristics of the psp
534 : !! if =2 : this input requires HGH characteristics of the psp
535 : !! if =3 : this input requires HFN characteristics of the psp
536 : !! ipsp=id in the array of the pseudo-potential.
537 : !! zion=nominal valence of atom as specified in psp file
538 : !!
539 : !! OUTPUT
540 : !! ekb(lnmax)=Kleinman-Bylander energy,
541 : !! {{\ \begin{equation}
542 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
543 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
544 : !! \end{equation} }}
545 : !! for each (l,n)
546 : !! if any, spin-orbit components begin at l=mpsang+1
547 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$ (hartree)
548 : !! ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
549 : !! second derivative from spline fit for each angular momentum and
550 : !! each projector; if any, spin-orbit components begin at l=mpsang+1
551 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
552 : !! or i=lmn (if useylm=1)
553 : !! nproj(mpssoang)=number of projection functions for each angular momentum
554 : !! vlspl(mqgrid_ff,2)=q^2 Vloc(q) and second derivatives from spline fit
555 : !!
556 : !! SIDE EFFECTS
557 : !! Input/output
558 : !! lmax : at input =value of lmax mentioned at the second line of the psp file
559 : !! at output= 1
560 : !! psps <type(pseudopotential_type)>=at output, values depending on the read
561 : !! pseudo are set.
562 : !! | lmnmax(IN)=if useylm=1, max number of (l,m,n) comp. over all type of psps
563 : !! | =if useylm=0, max number of (l,n) comp. over all type of psps
564 : !! | lnmax(IN)=max. number of (l,n) components over all type of psps
565 : !! | angular momentum of nonlocal pseudopotential
566 : !! | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
567 : !! | mpssoang(IN)= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
568 : !! | mqgrid_ff(IN)=dimension of q (or G) grid for arrays.
569 : !! | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
570 : !! | useylm(IN)=governs the way the nonlocal operator is to be applied:
571 : !! | 1=using Ylm, 0=using Legendre polynomials
572 : !!
573 : !! SOURCE
574 :
575 318 : subroutine psp3in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion)
576 :
577 : !Arguments ------------------------------------
578 : !scalars
579 : integer,intent(in) :: ipsp,pspso
580 : integer,intent(inout) :: lmax
581 : real(dp),intent(in) :: zion
582 : real(dp),intent(out) :: epsatm
583 : type(dataset_type),intent(in) :: dtset
584 : type(pseudopotential_type),intent(inout) :: psps
585 : !arrays
586 : integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang)
587 : real(dp),intent(inout) :: ekb(psps%lnmax),ffspl(psps%mqgrid_ff,2,psps%lnmax)!vz_i
588 : real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2)
589 :
590 : !Local variables-------------------------------
591 : !scalars
592 : integer :: iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nso
593 : real(dp) :: cc1,cc2,cc3,cc4,h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s
594 : real(dp) :: k11d,k11f,k11p,k22d,k22p,k33d,k33p,rloc
595 : real(dp) :: rrd,rrf,rrp,rrs,yp1,ypn
596 : character(len=500) :: msg,errmsg
597 : !arrays
598 318 : real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:)
599 318 : real(dp),allocatable :: ffspl_sr(:,:,:,:),work_space(:),work_spl(:)
600 :
601 : ! ***************************************************************************
602 :
603 : !Set various terms to 0 in case not defined below
604 : !HGH values
605 318 : rloc=zero ; rrs=zero ; h11p=zero ; k33p=zero ; k11d=zero;
606 318 : cc1=zero ; h11s=zero ; h22p=zero ; rrd=zero ; k22d=zero;
607 318 : cc2=zero ; h22s=zero ; h33p=zero ; h11d=zero ; k33d=zero;
608 318 : cc3=zero ; h33s=zero ; k11p=zero ; h22d=zero ; h11f=zero;
609 318 : cc4=zero ; rrp=zero ; k22p=zero ; h33d=zero ; k11f=zero;
610 318 : rrf=zero
611 1802 : nproj(1:psps%mpssoang)=0
612 :
613 : !Read and write different lines of the pseudopotential file
614 :
615 318 : read (tmp_unit,*,err=10,iomsg=errmsg) rloc,cc1,cc2,cc3,cc4
616 318 : write(msg, '(a,f12.7)' ) ' rloc=',rloc
617 318 : call wrtout(ab_out,msg,'COLL')
618 318 : call wrtout(std_out, msg,'COLL')
619 318 : write(msg, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' cc1 =',cc1,'; cc2 =',cc2,'; cc3 =',cc3,'; cc4 =',cc4
620 318 : call wrtout(ab_out,msg,'COLL')
621 318 : call wrtout(std_out, msg,'COLL')
622 :
623 : !For the time being, the s state line must be present and is read,
624 : !even for local pseudopotentials (zero must appear)
625 318 : read (tmp_unit,*,err=10,iomsg=errmsg) rrs,h11s,h22s,h33s
626 318 : write(msg, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrs =',rrs,'; h11s=',h11s,'; h22s=',h22s,'; h33s=',h33s
627 318 : call wrtout(ab_out,msg,'COLL')
628 318 : call wrtout(std_out, msg,'COLL')
629 :
630 318 : if (lmax > 0) then
631 :
632 298 : read (tmp_unit,*,err=10,iomsg=errmsg) rrp,h11p,h22p,h33p
633 298 : write(msg, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrp =',rrp,'; h11p=',h11p,'; h22p=',h22p,'; h33p=',h33p
634 298 : call wrtout(ab_out,msg,'COLL')
635 298 : call wrtout(std_out, msg,'COLL')
636 :
637 298 : read (tmp_unit,*,err=10,iomsg=errmsg) k11p,k22p,k33p
638 298 : write(msg, '(a,f12.7,a,f12.7,a,f12.7)' )' k11p=',k11p,'; k22p=',k22p,'; k33p=',k33p
639 298 : call wrtout(ab_out,msg,'COLL')
640 298 : call wrtout(std_out, msg,'COLL')
641 :
642 : end if
643 :
644 318 : if (lmax > 1) then
645 250 : read (tmp_unit,*,err=10,iomsg=errmsg) rrd,h11d,h22d,h33d
646 250 : write(msg, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrd =',rrd,'; h11d=',h11d,'; h22d=',h22d,'; h33d=',h33d
647 250 : call wrtout(ab_out,msg,'COLL')
648 250 : call wrtout(std_out, msg,'COLL')
649 :
650 250 : read (tmp_unit,*,err=10,iomsg=errmsg) k11d,k22d,k33d
651 250 : write(msg, '(a,f12.7,a,f12.7,a,f12.7)' )' k11d=',k11d,'; k22d=',k22d,'; k33d=',k33d
652 250 : call wrtout(ab_out,msg,'COLL')
653 250 : call wrtout(std_out, msg,'COLL')
654 : end if
655 :
656 318 : if (lmax > 2) then
657 2 : read (tmp_unit,*,err=10,iomsg=errmsg) rrf,h11f
658 2 : write(msg, '(a,f12.7,a,f12.7)' )' rrf =',rrf,'; h11f=',h11f
659 2 : call wrtout(ab_out,msg,'COLL')
660 2 : call wrtout(std_out, msg,'COLL')
661 :
662 2 : read (tmp_unit,*,err=10,iomsg=errmsg) k11f
663 2 : write(msg, '(a,f12.7)' )' k11f=',k11f
664 2 : call wrtout(ab_out,msg,'COLL')
665 2 : call wrtout(std_out, msg,'COLL')
666 : end if
667 :
668 318 : if (abs(h11s)>1.d-08) nproj(1)=1
669 318 : if (abs(h22s)>1.d-08) nproj(1)=2
670 318 : if (abs(h33s)>1.d-08) nproj(1)=3
671 :
672 318 : if (abs(h11p)>1.d-08) then
673 273 : nproj(2)=1
674 273 : if (lmax<1) then
675 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
676 0 : & ' psp3in : COMMENT -',ch10,&
677 0 : & ' input lmax=',lmax,' does not agree with input h11p=',h11p,ch10,&
678 0 : & ' setting lmax to 1'
679 0 : call wrtout(ab_out,msg,'COLL')
680 0 : call wrtout(std_out, msg,'COLL')
681 0 : lmax=1
682 : end if
683 : end if
684 :
685 318 : if (abs(h22p)>1.d-08) then
686 252 : nproj(2)=2
687 252 : if (lmax<1) then
688 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
689 0 : & ' psp3in : COMMENT -',ch10,&
690 0 : & ' input lmax=',lmax,' does not agree with input h22p=',h22p,ch10,&
691 0 : & ' setting lmax to 1'
692 0 : call wrtout(ab_out,msg,'COLL')
693 0 : call wrtout(std_out, msg,'COLL')
694 0 : lmax=1
695 : end if
696 : end if
697 :
698 318 : if (abs(h33p)>1.d-08) then
699 2 : nproj(2)=3
700 2 : if (lmax<1) then
701 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
702 0 : & ' psp3in : COMMENT -',ch10,&
703 0 : & ' input lmax=',lmax,' does not agree with input h33p=',h33p,ch10,&
704 0 : & ' setting lmax to 1'
705 0 : call wrtout(ab_out,msg,'COLL')
706 0 : call wrtout(std_out, msg,'COLL')
707 0 : lmax=1
708 : end if
709 : end if
710 :
711 318 : if (abs(h11d)>1.d-08) then
712 250 : nproj(3)=1
713 250 : if (lmax<2) then
714 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
715 0 : & ' psp3in : COMMENT -',ch10,&
716 0 : & ' input lmax=',lmax,' does not agree with input h11d=',h11d,ch10,&
717 0 : & ' setting lmax to 2'
718 0 : call wrtout(ab_out,msg,'COLL')
719 0 : call wrtout(std_out, msg,'COLL')
720 0 : lmax=2
721 : end if
722 : end if
723 :
724 318 : if (abs(h22d)>1.d-08) then
725 6 : nproj(3)=2
726 6 : if (lmax<2) then
727 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
728 0 : & ' psp3in : COMMENT -',ch10,&
729 0 : & ' input lmax=',lmax,' does not agree with input h22d=',h22d,ch10,&
730 0 : & ' setting lmax to 2'
731 0 : call wrtout(ab_out,msg,'COLL')
732 0 : call wrtout(std_out, msg,'COLL')
733 0 : lmax=2
734 : end if
735 : end if
736 :
737 318 : if (abs(h33d)>1.d-08) then
738 0 : nproj(3)=3
739 0 : if (lmax<2) then
740 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
741 0 : & ' psp3in : COMMENT -',ch10,&
742 0 : & ' input lmax=',lmax,' does not agree with input h33d=',h33d,ch10,&
743 0 : & ' setting lmax to 2'
744 0 : call wrtout(ab_out,msg,'COLL')
745 0 : call wrtout(std_out, msg,'COLL')
746 0 : lmax=2
747 : end if
748 : end if
749 :
750 318 : if (abs(h11f)>1.d-08) then
751 2 : nproj(4)=1
752 2 : if (lmax<3) then
753 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
754 0 : & ' psp3in : COMMENT -',ch10,&
755 0 : & ' input lmax=',lmax,' does not agree with input h11f=',h11f,ch10,&
756 0 : & ' setting lmax to 3'
757 0 : call wrtout(ab_out,msg,'COLL')
758 0 : call wrtout(std_out, msg,'COLL')
759 0 : lmax=3
760 : end if
761 : end if
762 :
763 318 : if(pspso/=0) then
764 :
765 28 : if (abs(k11p)>1.d-08) then
766 28 : nproj(psps%mpsang+1)=1
767 28 : if (lmax<1) then
768 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
769 0 : & ' psp3in : COMMENT -',ch10,&
770 0 : & ' input lmax=',lmax,' does not agree with input k11p=',k11p,ch10,&
771 0 : & ' setting lmax to 1'
772 0 : call wrtout(ab_out,msg,'COLL')
773 0 : call wrtout(std_out, msg,'COLL')
774 0 : lmax=1
775 : end if
776 : end if
777 :
778 28 : if (abs(k22p)>1.d-08) then
779 28 : nproj(psps%mpsang+1)=2
780 28 : if (lmax<1) then
781 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
782 0 : & ' psp3in : COMMENT -',ch10,&
783 0 : & ' input lmax=',lmax,' does not agree with input k22p=',k22p,ch10,&
784 0 : & ' setting lmax to 1'
785 0 : call wrtout(ab_out,msg,'COLL')
786 0 : call wrtout(std_out, msg,'COLL')
787 0 : lmax=1
788 : end if
789 : end if
790 :
791 :
792 28 : if (abs(k33p)>1.d-08) then
793 2 : nproj(psps%mpsang+1)=3
794 2 : if (lmax<1) then
795 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
796 0 : & ' psp3in : COMMENT -',ch10,&
797 0 : & ' input lmax=',lmax,' does not agree with input k33p=',k33p,ch10,&
798 0 : & ' setting lmax to 1'
799 0 : call wrtout(ab_out,msg,'COLL')
800 0 : call wrtout(std_out, msg,'COLL')
801 0 : lmax=1
802 : end if
803 : end if
804 :
805 28 : if (abs(k11d)>1.d-08) then
806 26 : nproj(psps%mpsang+2)=1
807 26 : if (lmax<2) then
808 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
809 0 : & ' psp3in : COMMENT -',ch10,&
810 0 : & ' input lmax=',lmax,' does not agree with input k11d=',k11d,ch10,&
811 0 : & ' setting lmax to 2'
812 0 : call wrtout(ab_out,msg,'COLL')
813 0 : call wrtout(std_out, msg,'COLL')
814 0 : lmax=2
815 : end if
816 : end if
817 :
818 28 : if (abs(k22d)>1.d-08) then
819 4 : nproj(psps%mpsang+2)=2
820 4 : if (lmax<2) then
821 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
822 0 : & ' psp3in : COMMENT -',ch10,&
823 0 : & ' input lmax=',lmax,' does not agree with input k22d=',k22d,ch10,&
824 0 : & ' setting lmax to 2'
825 0 : call wrtout(ab_out,msg,'COLL')
826 0 : call wrtout(std_out, msg,'COLL')
827 0 : lmax=2
828 : end if
829 : end if
830 :
831 28 : if (abs(k33d)>1.d-08) then
832 0 : nproj(psps%mpsang+2)=3
833 0 : if (lmax<2) then
834 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
835 0 : & ' psp3in : COMMENT -',ch10,&
836 0 : & ' input lmax=',lmax,' does not agree with input k33d=',k33d,ch10,&
837 0 : & ' setting lmax to 2'
838 0 : call wrtout(ab_out,msg,'COLL')
839 0 : call wrtout(std_out, msg,'COLL')
840 0 : lmax=2
841 : end if
842 : end if
843 :
844 28 : if (abs(k11f)>1.d-08) then
845 0 : nproj(psps%mpsang+3)=1
846 0 : if (lmax<3) then
847 0 : write(msg, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
848 0 : & ' psp3in : COMMENT -',ch10,&
849 0 : & ' input lmax=',lmax,' does not agree with input k11f=',k11f,ch10,&
850 0 : & ' setting lmax to 3'
851 0 : call wrtout(ab_out,msg,'COLL')
852 0 : call wrtout(std_out, msg,'COLL')
853 0 : lmax=3
854 : end if
855 : end if
856 :
857 : end if
858 :
859 : !Store the coefficients.
860 318 : psps%gth_params%set(ipsp) = .true.
861 2544 : psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, zero, zero /)
862 2544 : psps%gth_params%psppar(1, :, ipsp) = (/ rrs, h11s, h22s, h33s, zero, zero, zero /)
863 2544 : psps%gth_params%psppar(2, :, ipsp) = (/ rrp, h11p, h22p, h33p, zero, zero, zero /)
864 2544 : psps%gth_params%psppar(3, :, ipsp) = (/ rrd, h11d, h22d, h33d, zero, zero, zero /)
865 2544 : psps%gth_params%psppar(4, :, ipsp) = (/ rrf, h11f, zero, zero, zero, zero, zero /)
866 :
867 : !Store the k coefficients
868 1272 : psps%gth_params%psp_k_par(1, :, ipsp) = (/ zero, zero, zero /)
869 1272 : psps%gth_params%psp_k_par(2, :, ipsp) = (/ k11p, k22p, k33p /)
870 1272 : psps%gth_params%psp_k_par(3, :, ipsp) = (/ k11d, k22d, k33d /)
871 1272 : psps%gth_params%psp_k_par(4, :, ipsp) = (/ k11f, zero, zero /)
872 :
873 : !Additionnal wavelet parameters
874 318 : if (dtset%usewvl == 1) then
875 0 : call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
876 : end if
877 :
878 : !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn
879 318 : nso=1;if(pspso/=0) nso=2
880 19295 : index=0;iln=0;indlmn(:,:)=0
881 664 : do nn=1,nso
882 1586 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
883 1268 : if (nproj(ipsang)>0) then
884 875 : ll=ipsang-(nn-1)*lmax-1
885 2553 : do kk=1,nproj(ipsang)
886 1678 : iln=iln+1
887 4249 : do mm=1,2*ll*psps%useylm+1
888 1696 : index=index+1
889 1696 : indlmn(1,index)=ll
890 1696 : indlmn(2,index)=mm-ll*psps%useylm-1
891 1696 : indlmn(3,index)=kk
892 1696 : indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
893 1696 : indlmn(5,index)=iln
894 3374 : indlmn(6,index)=nn
895 : end do
896 : end do
897 : end if
898 : end do
899 : end do
900 :
901 954 : ABI_MALLOC(dvlspl,(psps%mqgrid_ff))
902 : !First, the local potential -- compute on q grid and fit spline
903 318 : call psp2lo(cc1,cc2,cc3,cc4,dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,vlspl(:,1),rloc,.true.,yp1,ypn,zion)
904 318 : ABI_FREE(dvlspl)
905 :
906 : !DEBUG
907 : !write(std_out,*)' psp3in : after psp2lo '
908 : !ENDDEBUG
909 :
910 : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
911 954 : ABI_MALLOC(work_space,(psps%mqgrid_ff))
912 636 : ABI_MALLOC(work_spl,(psps%mqgrid_ff))
913 318 : call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl)
914 958908 : vlspl(:,2)=work_spl(:)
915 318 : ABI_FREE(work_space)
916 318 : ABI_FREE(work_spl)
917 :
918 : !Second, compute KB energies and form factors and fit splines
919 2997 : ekb(:)=zero
920 :
921 : !Check if any nonlocal projectors are being used
922 1802 : mproj=maxval(nproj)
923 318 : if (mproj>0) then
924 :
925 1192 : ABI_MALLOC(ekb_sr,(psps%mpsang,mproj))
926 1490 : ABI_MALLOC(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj))
927 894 : ABI_MALLOC(ekb_so,(psps%mpsang,mproj))
928 1192 : ABI_MALLOC(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj))
929 :
930 : call psp3nl(ekb_sr,ffspl_sr,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
931 298 : & h33d,h11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs)
932 298 : if(pspso/=0) then
933 : call psp3nl(ekb_so,ffspl_so,zero,zero,zero,k11p,k22p,k33p,k11d,&
934 28 : & k22d,k33d,k11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs)
935 : end if
936 :
937 :
938 : ! Convert ekb and ffspl
939 298 : iln0=0
940 2975 : do jj=1,psps%lmnmax
941 2677 : iln=indlmn(5,jj)
942 2975 : if (iln>iln0) then
943 1678 : iln0=iln
944 1678 : if (indlmn(6,jj)<=1) then
945 1594 : ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj))
946 9619058 : ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj))
947 : else
948 84 : ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj))
949 504420 : ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj))
950 : end if
951 : end if
952 : end do
953 :
954 298 : ABI_FREE(ekb_sr)
955 298 : ABI_FREE(ffspl_sr)
956 298 : ABI_FREE(ekb_so)
957 298 : ABI_FREE(ffspl_so)
958 : end if
959 :
960 : return
961 :
962 : ! Handle IO error
963 : 10 continue
964 0 : ABI_ERROR(errmsg)
965 :
966 : end subroutine psp3in
967 : !!***
968 :
969 : !!****f* m_psp_hgh/psp3nl
970 : !! NAME
971 : !! psp3nl
972 : !!
973 : !! FUNCTION
974 : !! Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998).
975 : !! Uses Gaussians for fully nonlocal form, analytic expressions.
976 : !!
977 : !! INPUTS
978 : !! h11s=factor defining strength of 1st projector for l=0 channel
979 : !! h22s=factor defining strength of 2nd projector for l=0 channel
980 : !! h33s=factor defining strength of 3rd projector for l=0 channel
981 : !! h11p=factor defining strength of 1st projector for l=1 channel
982 : !! h22p=factor defining strength of 2nd projector for l=1 channel
983 : !! h33p=factor defining strength of 2nd projector for l=1 channel
984 : !! h11d=factor defining strength of 1st projector for l=2 channel
985 : !! h22d=factor defining strength of 2nd projector for l=2 channel
986 : !! h33d=factor defining strength of 2nd projector for l=2 channel
987 : !! h11f=factor defining strength of 1st projector for l=3 channel
988 : !! mproj=maximum number of projectors in any channel
989 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
990 : !! mqgrid=number of grid points for qgrid
991 : !! qgrid(mqgrid)=array of |G| values
992 : !! rrd=core radius for d channel (bohr)
993 : !! rrf=core radius for f channel (bohr)
994 : !! rrp=core radius for p channel (bohr)
995 : !! rrs=core radius for s channel (bohr)
996 : !!
997 : !! OUTPUT
998 : !! ekb(mpsang,mproj)=Kleinman-Bylander energies
999 : !! ffspl(mqgrid,2,mpssang,mproj)=Kleinman-Bylander form factor f_l(q) and
1000 : !! second derivative from spline fit for each angular momentum and
1001 : !! each projectors
1002 : !!
1003 : !! SOURCE
1004 :
1005 326 : subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
1006 326 : & h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs)
1007 :
1008 : !Arguments ------------------------------------
1009 : !scalars
1010 : integer,intent(in) :: mproj,mpsang,mqgrid
1011 : real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd
1012 : real(dp),intent(in) :: rrf,rrp,rrs
1013 : !arrays
1014 : real(dp),intent(in) :: qgrid(mqgrid)
1015 : real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
1016 :
1017 : !Local variables-------------------------------
1018 : !scalars
1019 : integer :: info,iproj,iqgrid,ldz,mu,nproj,nu
1020 : real(dp) :: qmax
1021 : character(len=500) :: msg
1022 : character :: jobz,uplo
1023 : !arrays
1024 : real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
1025 326 : real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
1026 :
1027 : ! *************************************************************************
1028 :
1029 1630 : ABI_MALLOC(ppspl,(mqgrid,2,mpsang,mproj))
1030 978 : ABI_MALLOC(work,(mqgrid))
1031 :
1032 326 : qmax=qgrid(mqgrid)
1033 326 : jobz='v'
1034 326 : uplo='u'
1035 3863 : ekb(:,:)=0.0d0
1036 :
1037 : !---------------------------------------------------------------
1038 : !Treat s channel
1039 :
1040 326 : nproj=0
1041 326 : ap(:,:)=0.0d0
1042 : !If there is at least one s-projector
1043 326 : if ( abs(h11s) >= 1.0d-8 ) then
1044 298 : nproj=1 ; ldz=1 ; ap(1,1)=h11s
1045 : end if
1046 326 : nproj=1
1047 : !If there is a second projector
1048 326 : if ( abs(h22s) >= 1.0d-8 ) then
1049 266 : nproj=2 ; ldz=2 ; ap(1,3)=h22s
1050 266 : ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s
1051 : end if
1052 : !If there is a third projector
1053 326 : if ( abs(h33s) >= 1.0d-8 ) then
1054 245 : nproj=3 ; ldz=3 ; ap(1,6)=h33s
1055 245 : ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s
1056 245 : ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s
1057 : end if
1058 :
1059 326 : if(nproj/=0)then
1060 :
1061 1304 : ABI_MALLOC(uu,(nproj,nproj))
1062 1304 : ABI_MALLOC(zz,(2,nproj,nproj))
1063 :
1064 326 : if (nproj > 1) then
1065 266 : call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1066 3332 : uu(:,:)=zz(1,:,:)
1067 : else
1068 60 : ww(1)=h11s
1069 60 : uu(1,1)=1.0d0
1070 : end if
1071 :
1072 : ! Initialization of ekb, and spline fitting
1073 1163 : do iproj=1,nproj
1074 837 : ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2
1075 837 : if(iproj==1)then
1076 982886 : do iqgrid=1,mqgrid
1077 982886 : ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
1078 : end do
1079 326 : yp1j(1)=0.d0
1080 326 : ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2)
1081 511 : else if(iproj==2)then
1082 802766 : do iqgrid=1,mqgrid
1083 : ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) &
1084 : & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) &
1085 802766 : & *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 )
1086 : end do
1087 266 : yp1j(2)=0.0d0
1088 : ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax &
1089 266 : & *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2)
1090 245 : else if(iproj==3)then
1091 739104 : do iqgrid=1,mqgrid
1092 : ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
1093 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
1094 : & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + &
1095 739104 : & (two_pi*qgrid(iqgrid)*rrs)**4)
1096 : end do
1097 245 : yp1j(3)=0.0d0
1098 : ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
1099 245 : & (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4)
1100 : end if
1101 : call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
1102 1163 : & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
1103 : end do
1104 :
1105 : ! Linear combination using the eigenvectors
1106 5362935 : ffspl(:,:,1,:)=0.0d0
1107 1163 : do mu=1,nproj
1108 3512 : do nu=1,nproj
1109 7087541 : do iqgrid=1,mqgrid
1110 : ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) &
1111 21255414 : & +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu)
1112 : end do
1113 : end do
1114 : end do
1115 :
1116 326 : ABI_FREE(uu)
1117 326 : ABI_FREE(zz)
1118 : end if ! End condition on nproj(/=0)
1119 :
1120 : !--------------------------------------------------------------------
1121 : !Now treat p channel
1122 :
1123 326 : nproj=0
1124 326 : ap(:,:)=0.0d0
1125 : !If there is at least one projector
1126 326 : if ( abs(h11p) >= 1.0d-8 ) then
1127 301 : nproj=1 ; ldz=1 ; ap(1,1)=h11p
1128 : end if
1129 : !If there is a second projector
1130 326 : if ( abs(h22p) >= 1.0d-8 ) then
1131 280 : nproj=2 ; ldz=2 ; ap(1,3)=h22p
1132 280 : ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p
1133 : end if
1134 : !If there is a third projector
1135 326 : if ( abs(h33p) >= 1.0d-8 ) then
1136 4 : nproj=3 ; ldz=3 ; ap(1,6)=h33p
1137 4 : ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p
1138 4 : ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p
1139 : end if
1140 :
1141 326 : if(nproj/=0)then
1142 :
1143 1204 : ABI_MALLOC(uu,(nproj,nproj))
1144 1204 : ABI_MALLOC(zz,(2,nproj,nproj))
1145 :
1146 301 : if (nproj > 1) then
1147 280 : call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1148 1984 : uu(:,:)=zz(1,:,:)
1149 : else
1150 21 : ww(1)=h11p
1151 21 : uu(1,1)=1.0d0
1152 : end if
1153 :
1154 : ! Initialization of ekb, and spline fitting
1155 886 : do iproj=1,nproj
1156 585 : ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2
1157 585 : if(iproj==1)then
1158 907836 : do iqgrid=1,mqgrid
1159 : ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
1160 907836 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid))
1161 : end do
1162 301 : yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
1163 : ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*&
1164 301 : & (1.0d0/sqrt(3.0d0))
1165 284 : else if(iproj==2)then
1166 844174 : do iqgrid=1,mqgrid
1167 : ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
1168 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
1169 844174 : & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2)
1170 : end do
1171 280 : yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
1172 : ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
1173 280 : & (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0)
1174 4 : else if(iproj==3)then
1175 15622 : do iqgrid=1,mqgrid
1176 : ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
1177 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
1178 : & (two_pi*qgrid(iqgrid))*&
1179 15622 : & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4)
1180 : end do
1181 4 : yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
1182 : ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
1183 : & (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - &
1184 4 : & (two_pi*qmax*rrp)**6)
1185 : end if
1186 : call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
1187 886 : & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
1188 : end do
1189 :
1190 : ! Linear combination using the eigenvectors
1191 5212785 : ffspl(:,:,2,:)=0.0d0
1192 886 : do mu=1,nproj
1193 2047 : do nu=1,nproj
1194 3519053 : do iqgrid=1,mqgrid
1195 : ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) &
1196 10553082 : & +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu)
1197 : end do
1198 : end do
1199 : end do
1200 :
1201 301 : ABI_FREE(uu)
1202 301 : ABI_FREE(zz)
1203 : end if ! End condition on nproj(/=0)
1204 :
1205 : !-----------------------------------------------------------------------
1206 : !Now treat d channel.
1207 :
1208 326 : nproj=0
1209 326 : ap(:,:)=0.0d0
1210 : !If there is at least one projector
1211 326 : if ( abs(h11d) >= 1.0d-8 ) then
1212 276 : nproj=1 ; ldz=1 ; ap(1,1)=h11d
1213 : end if
1214 : !If there is a second projector
1215 326 : if ( abs(h22d) >= 1.0d-8 ) then
1216 10 : nproj=2 ; ldz=2 ; ap(1,3)=h22d
1217 10 : ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d
1218 : end if
1219 : !If there is a third projector. Warning : only two projectors are allowed.
1220 326 : if ( abs(h33d) >= 1.0d-8 ) then
1221 : write(msg, '(a,a,a)' )&
1222 0 : & ' only two d-projectors are allowed ',ch10,&
1223 0 : & ' Action: check your pseudopotential file.'
1224 0 : ABI_ERROR(msg)
1225 : ! nproj=3 ; ldz=3 ; ap(1,6)=h33d
1226 : ! ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d
1227 : ! ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d
1228 : end if
1229 :
1230 326 : if(nproj/=0)then
1231 :
1232 1104 : ABI_MALLOC(uu,(nproj,nproj))
1233 1104 : ABI_MALLOC(zz,(2,nproj,nproj))
1234 :
1235 276 : if (nproj > 1) then
1236 10 : call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1237 70 : uu(:,:)=zz(1,:,:)
1238 : else
1239 266 : ww(1)=h11d
1240 266 : uu(1,1)=1.0d0
1241 : end if
1242 :
1243 : ! Initialization of ekb, and spline fitting
1244 562 : do iproj=1,nproj
1245 286 : ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2
1246 286 : if(iproj==1)then
1247 828552 : do iqgrid=1,mqgrid
1248 : ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
1249 828552 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2
1250 : end do
1251 276 : yp1j(1)=0.0d0
1252 : ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
1253 276 : & exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2)
1254 10 : else if(iproj==2)then
1255 30020 : do iqgrid=1,mqgrid
1256 : ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
1257 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * &
1258 30020 : & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2)
1259 : end do
1260 10 : yp1j(2)=0.0d0
1261 : ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* &
1262 10 : & qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0)
1263 : end if
1264 : call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
1265 562 : & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
1266 : end do
1267 :
1268 : ! Linear combination using the eigenvectors
1269 4918371 : ffspl(:,:,3,:)=0.0d0
1270 562 : do mu=1,nproj
1271 868 : do nu=1,nproj
1272 918898 : do iqgrid=1,mqgrid
1273 : ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) &
1274 2755224 : & +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu)
1275 : end do
1276 : end do
1277 : end do
1278 :
1279 276 : ABI_FREE(uu)
1280 276 : ABI_FREE(zz)
1281 : end if ! End condition on nproj(/=0)
1282 :
1283 : !-----------------------------------------------------------------------
1284 : !Treat now f channel (max one projector ! - so do not use ppspl)
1285 :
1286 : !l=3 first projector
1287 326 : if (abs(h11f)>1.d-12) then
1288 2 : ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2
1289 6004 : do iqgrid=1,mqgrid
1290 : ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
1291 6004 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2)
1292 : end do
1293 : ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
1294 2 : yp1j(1)=0d0
1295 : ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*&
1296 2 : & (3.0d0-(two_pi*qmax*rrf)**2)
1297 : ! Fit spline to get second derivatives by spline fit
1298 : call spline(qgrid,ffspl(:,1,4,1),mqgrid,&
1299 2 : & yp1j(1),ypnj(1),ffspl(:,2,4,1))
1300 : end if
1301 :
1302 : !-----------------------------------------------------------------------
1303 :
1304 326 : ABI_FREE(ppspl)
1305 326 : ABI_FREE(work)
1306 :
1307 326 : end subroutine psp3nl
1308 : !!***
1309 :
1310 : !!****f* m_psp_hgh/psp10in
1311 : !! NAME
1312 : !! psp10in
1313 : !!
1314 : !! FUNCTION
1315 : !! Initialize pspcod=10 pseudopotentials (formalism is the same as in HGH psps
1316 : !! PRB58,3641(1998) [[cite:Hartwigsen1998]], but the full h and k matrices are read, allowing for using
1317 : !! also subsequent developments such as Theor. Chem. Acc. 114, 145 (2005) [[cite:Dolg2005]]:
1318 : !! continue to read the file, then compute the corresponding
1319 : !! local and non-local potentials.
1320 : !!
1321 : !! INPUTS
1322 : !! dtset <type(dataset_type)>=all input variables in this dataset
1323 : !! pspso=spin-orbit characteristics, govern the content of ffspl and ekb
1324 : !! if =0 : this input requires NO spin-orbit characteristics of the psp
1325 : !! if =2 : this input requires HGH characteristics of the psp
1326 : !! if =3 : this input requires HFN characteristics of the psp
1327 : !! ipsp=id in the array of the pseudo-potential.
1328 : !! zion=nominal valence of atom as specified in psp file
1329 : !!
1330 : !! OUTPUT
1331 : !! ekb(lnmax)=Kleinman-Bylander energy,
1332 : !! {{\ \begin{equation}
1333 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
1334 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
1335 : !! \end{equation} }}
1336 : !! for each (l,n)
1337 : !! if any, spin-orbit components begin at l=mpsang+1
1338 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$ (hartree)
1339 : !! ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
1340 : !! second derivative from spline fit for each angular momentum and
1341 : !! each projector; if any, spin-orbit components begin at l=mpsang+1
1342 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
1343 : !! or i=lmn (if useylm=1)
1344 : !! nproj(mpssoang)=number of projection functions for each angular momentum
1345 : !! vlspl(mqgrid_ff,2)=q^2 Vloc(q) and second derivatives from spline fit
1346 : !!
1347 : !! SIDE EFFECTS
1348 : !! Input/output
1349 : !! lmax : at input =value of lmax mentioned at the second line of the psp file
1350 : !! at output= 1
1351 : !! psps <type(pseudopotential_type)>=at output, values depending on the read
1352 : !! pseudo are set.
1353 : !! | lmnmax(IN)=if useylm=1, max number of (l,m,n) comp. over all type of psps
1354 : !! | =if useylm=0, max number of (l,n) comp. over all type of psps
1355 : !! | lnmax(IN)=max. number of (l,n) components over all type of psps
1356 : !! | angular momentum of nonlocal pseudopotential
1357 : !! | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
1358 : !! | mpssoang(IN)= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
1359 : !! | mqgrid_ff(IN)=dimension of q (or G) grid for arrays.
1360 : !! | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
1361 : !! | useylm(IN)=governs the way the nonlocal operator is to be applied:
1362 : !! | 1=using Ylm, 0=using Legendre polynomials
1363 : !!
1364 : !! SOURCE
1365 :
1366 12 : subroutine psp10in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion)
1367 :
1368 : !Arguments ------------------------------------
1369 : !scalars
1370 : integer,intent(in) :: ipsp,pspso
1371 : integer,intent(inout) :: lmax
1372 : real(dp),intent(in) :: zion
1373 : real(dp),intent(out) :: epsatm
1374 : type(dataset_type),intent(in) :: dtset
1375 : type(pseudopotential_type),intent(inout) :: psps
1376 : !arrays
1377 : integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang)
1378 : real(dp),intent(out) :: ekb(psps%lnmax) !vz_i
1379 : real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i
1380 : real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2)
1381 :
1382 : !Local variables-------------------------------
1383 : !scalars
1384 : integer :: ii,iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nnonloc,nprl,nso
1385 : real(dp) :: rloc,yp1,ypn
1386 : character(len=500) :: msg,errmsg
1387 : !arrays
1388 12 : integer,allocatable :: dummy_nproj(:)
1389 : real(dp) :: cc(4)
1390 12 : real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:)
1391 12 : real(dp),allocatable :: ffspl_sr(:,:,:,:),hij(:,:,:),kij(:,:,:),rr(:)
1392 12 : real(dp),allocatable :: work_space(:),work_spl(:)
1393 :
1394 : ! ***************************************************************************
1395 :
1396 : !Set various terms to 0 in case not defined below
1397 : !HGH values
1398 12 : rloc=zero ; cc(:)=zero
1399 50 : nproj(1:psps%mpssoang)=0
1400 :
1401 : !Read and write different lines of the pseudopotential file
1402 :
1403 12 : read (tmp_unit,*,err=10,iomsg=errmsg) rloc,nn,(cc(jj),jj=1,nn)
1404 12 : write(msg, '(a,f12.7)' ) ' rloc=',rloc
1405 12 : call wrtout(ab_out,msg,'COLL')
1406 12 : call wrtout(std_out, msg,'COLL')
1407 12 : write(msg, '(a,i1,a,4f12.7)' )' cc(1:',nn,')=',(cc(jj),jj=1,nn)
1408 12 : call wrtout(ab_out,msg,'COLL')
1409 12 : call wrtout(std_out, msg,'COLL')
1410 :
1411 : !Read the number of the non-local projectors
1412 12 : read (tmp_unit,*,err=10,iomsg=errmsg) nnonloc
1413 12 : if (nnonloc/=lmax+1) then
1414 0 : write(msg, '(a,a,a,a,i5,a,i5,a,a,a,a,i5)' ) ch10,&
1415 0 : & ' psp10in : COMMENT -',ch10,&
1416 0 : & ' input lmax=',lmax,' does not agree with input nnonloc=',nnonloc,ch10,&
1417 0 : & ' which has to be lmax+1.',ch10,&
1418 0 : & ' Setting lmax to ',nnonloc-1
1419 0 : call wrtout(ab_out,msg,'COLL')
1420 0 : call wrtout(std_out, msg,'COLL')
1421 0 : lmax=1
1422 : end if
1423 36 : ABI_MALLOC(rr,(0:lmax))
1424 60 : ABI_MALLOC(hij,(0:lmax,3,3))
1425 24 : ABI_MALLOC(kij,(0:lmax,3,3))
1426 799 : rr(:)=zero; hij(:,:,:)=zero; kij(:,:,:)=zero
1427 :
1428 : !Read and echo the coefficients of non-local projectors
1429 37 : prjloop: do ll=0,lmax
1430 25 : read (tmp_unit,*,err=10,iomsg=errmsg) rr(ll),nprl,(hij(ll,1,jj),jj=1,nprl)
1431 29 : do ii=2,nprl
1432 29 : read (tmp_unit,*,err=10,iomsg=errmsg) (hij(ll,ii,jj),jj=ii,nprl)
1433 : end do
1434 25 : nproj(ll+1)=nprl
1435 : write(msg, '(a,i3,a,f12.7,2a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )&
1436 25 : & ' for angular momentum l =',ll,' r(l) =',rr(ll),ch10,&
1437 25 : & ' h11, h12, h13 =', (hij(ll,1,jj),jj=1,3),ch10,&
1438 25 : & ' h22, h23 =', (hij(ll,2,jj),jj=2,3),ch10,&
1439 50 : & ' h33 =', (hij(ll,3,jj),jj=3,3)
1440 25 : call wrtout(ab_out,msg,'COLL')
1441 25 : call wrtout(std_out, msg,'COLL')
1442 25 : if (ll==0) cycle
1443 27 : do ii=1,nprl
1444 27 : read (tmp_unit,*,err=10,iomsg=errmsg) (kij(ll,ii,jj),jj=ii,nprl)
1445 : end do
1446 : write(msg, '(a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )&
1447 13 : & ' k11, k12, k13 =', (kij(ll,1,jj),jj=1,3),ch10,&
1448 13 : & ' k22, k23 =', (kij(ll,2,jj),jj=2,3),ch10,&
1449 26 : & ' k33 =', (kij(ll,3,jj),jj=3,3)
1450 13 : call wrtout(ab_out,msg,'COLL')
1451 37 : call wrtout(std_out, msg,'COLL')
1452 : end do prjloop
1453 :
1454 12 : if(pspso/=0) then
1455 :
1456 : ! MJV 10/2008: is this correct? For the normal HGH psp there are cases
1457 : ! where there are more SO projectors than SR ones! e.g. Pb with 12 electrons.
1458 0 : do ll=1,lmax
1459 0 : nproj(psps%mpsang+ll)=nproj(ll+1)
1460 : end do
1461 :
1462 : end if
1463 :
1464 : !Store the coefficients.
1465 12 : psps%gth_params%set(ipsp) = .true.
1466 96 : psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc(1), cc(2), cc(3), cc(4), zero, zero /)
1467 60 : do ii=1,4
1468 48 : ll=ii-1
1469 60 : if (ll>lmax) then
1470 184 : psps%gth_params%psppar(ii,:,ipsp) = (/ zero, zero, zero, zero, zero, zero, zero /)
1471 : else
1472 : psps%gth_params%psppar(ii,:,ipsp) =&
1473 200 : & (/ rr(ll), hij(ll,1,1), hij(ll,2,2), hij(ll,3,3), hij(ll,1,2), hij(ll,1,3), hij(ll,2,3) /)
1474 : end if
1475 : end do
1476 :
1477 : !Additionnal wavelet parameters
1478 12 : if (dtset%usewvl == 1) then
1479 0 : call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
1480 : end if
1481 :
1482 : !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn
1483 12 : nso=1;if(pspso/=0) nso=2
1484 313 : index=0;iln=0;indlmn(:,:)=0
1485 24 : do nn=1,nso
1486 49 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
1487 37 : if (nproj(ipsang)>0) then
1488 25 : ll=ipsang-(nn-1)*lmax-1
1489 54 : do kk=1,nproj(ipsang)
1490 29 : iln=iln+1
1491 83 : do mm=1,2*ll*psps%useylm+1
1492 29 : index=index+1
1493 29 : indlmn(1,index)=ll
1494 29 : indlmn(2,index)=mm-ll*psps%useylm-1
1495 29 : indlmn(3,index)=kk
1496 29 : indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
1497 29 : indlmn(5,index)=iln
1498 58 : indlmn(6,index)=nn
1499 : end do
1500 : end do
1501 : end if
1502 : end do
1503 : end do
1504 :
1505 36 : ABI_MALLOC(dvlspl,(psps%mqgrid_ff))
1506 : !First, the local potential -- compute on q grid and fit spline
1507 : call psp2lo(cc(1),cc(2),cc(3),cc(4),dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,&
1508 12 : & vlspl(:,1),rloc,.true.,yp1,ypn,zion)
1509 12 : ABI_FREE(dvlspl)
1510 :
1511 : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
1512 36 : ABI_MALLOC(work_space,(psps%mqgrid_ff))
1513 24 : ABI_MALLOC(work_spl,(psps%mqgrid_ff))
1514 12 : call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl)
1515 36024 : vlspl(:,2)=work_spl(:)
1516 12 : ABI_FREE(work_space)
1517 12 : ABI_FREE(work_spl)
1518 :
1519 : !Second, compute KB energies and form factors and fit splines
1520 55 : ekb(:)=zero
1521 :
1522 : !Check if any nonlocal projectors are being used
1523 50 : mproj=maxval(nproj)
1524 :
1525 12 : if (mproj>0) then
1526 :
1527 48 : ABI_MALLOC(ekb_sr,(psps%mpsang,mproj))
1528 60 : ABI_MALLOC(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj))
1529 36 : ABI_MALLOC(ekb_so,(psps%mpsang,mproj))
1530 48 : ABI_MALLOC(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj))
1531 :
1532 : call psp10nl(ekb_sr,ffspl_sr,hij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,&
1533 12 : & nproj,psps%qgrid_ff,rr)
1534 12 : if(pspso/=0) then
1535 0 : ABI_MALLOC(dummy_nproj,(psps%mpsang))
1536 0 : dummy_nproj(1)=0
1537 0 : do ll=1,lmax
1538 0 : dummy_nproj(ll+1)=nproj(psps%mpsang+ll)
1539 : end do
1540 : call psp10nl(ekb_so,ffspl_so,kij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,&
1541 0 : & dummy_nproj,psps%qgrid_ff,rr)
1542 0 : ABI_FREE(dummy_nproj)
1543 : end if
1544 :
1545 : ! Convert ekb and ffspl
1546 12 : iln0=0
1547 55 : do jj=1,psps%lmnmax
1548 43 : iln=indlmn(5,jj)
1549 55 : if (iln>iln0) then
1550 29 : iln0=iln
1551 29 : if (indlmn(6,jj)<=1) then
1552 29 : ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj))
1553 174145 : ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj))
1554 : else
1555 0 : ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj))
1556 0 : ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj))
1557 : end if
1558 : end if
1559 : end do
1560 :
1561 12 : ABI_FREE(ekb_sr)
1562 12 : ABI_FREE(ffspl_sr)
1563 12 : ABI_FREE(ekb_so)
1564 12 : ABI_FREE(ffspl_so)
1565 : end if
1566 :
1567 12 : ABI_FREE(rr)
1568 12 : ABI_FREE(hij)
1569 12 : ABI_FREE(kij)
1570 :
1571 12 : return
1572 :
1573 : ! Handle IO error
1574 : 10 continue
1575 0 : ABI_ERROR(errmsg)
1576 :
1577 12 : end subroutine psp10in
1578 : !!***
1579 :
1580 : !!****f* m_psp_hgh/psp10nl
1581 : !! NAME
1582 : !! psp10nl
1583 : !!
1584 : !! FUNCTION
1585 : !! Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998).
1586 : !! Uses Gaussians for fully nonlocal form, analytic expressions.
1587 : !!
1588 : !! INPUTS
1589 : !! hij(0:lmax,3,3)=factor defining strength of (max 3) projectors for each
1590 : !! angular momentum channel l among 0, 1, ..., lmax
1591 : !! lmax=maximum angular momentum
1592 : !! mproj=maximum number of projectors in any channel
1593 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
1594 : !! mqgrid=number of grid points for qgrid
1595 : !! nproj(1:lmax+1)=number of projectors in any channel
1596 : !! qgrid(mqgrid)=array of |G| values
1597 : !! rr(0:lmax)=core radius for each 0<l<lmax channel (bohr)
1598 : !!
1599 : !! OUTPUT
1600 : !! ekb(mpsang,mproj)=Kleinman-Bylander energies
1601 : !! ffspl(mqgrid,2,mpssang,mproj)=Kleinman-Bylander form factor f_l(q) and
1602 : !! second derivative from spline fit for each angular momentum and
1603 : !! each projectors
1604 : !!
1605 : !! SOURCE
1606 :
1607 12 : subroutine psp10nl(ekb,ffspl,hij,lmax,mproj,mpsang,mqgrid,nproj,qgrid,rr)
1608 :
1609 : !Arguments ------------------------------------
1610 : !scalars
1611 : integer,intent(in) :: lmax,mproj,mpsang,mqgrid
1612 : !arrays
1613 : integer,intent(in) :: nproj(mpsang)
1614 : real(dp),intent(in) :: hij(0:lmax,3,3),qgrid(mqgrid),rr(0:lmax)
1615 : real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
1616 :
1617 : !Local variables-------------------------------
1618 : !scalars
1619 : integer :: info,ipack,iproj,iqgrid,jproj,ll,numproj
1620 : real(dp) :: qmax,rrl
1621 : character(len=500) :: msg
1622 : character :: jobz,uplo
1623 : !arrays
1624 : real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
1625 12 : real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
1626 :
1627 : ! *************************************************************************
1628 :
1629 60 : ABI_MALLOC(ppspl,(mqgrid,2,mpsang,mproj))
1630 36 : ABI_MALLOC(work,(mqgrid))
1631 :
1632 12 : qmax=qgrid(mqgrid)
1633 12 : jobz='v'
1634 12 : uplo='u'
1635 60 : ekb(:,:)=zero
1636 :
1637 37 : lloop: do ll=0,lmax
1638 25 : ap(:,:)=zero
1639 25 : numproj=nproj(ll+1)
1640 :
1641 : ! Fill up the matrix in packed storage
1642 54 : prjloop: do jproj=1,numproj
1643 88 : priloop: do iproj=1,jproj
1644 34 : ipack=iproj+(jproj-1)*jproj/2
1645 34 : if(mod((jproj-1)*jproj,2)/=0) then
1646 0 : ABI_ERROR("odd")
1647 : end if
1648 63 : ap(1,ipack)=hij(ll,iproj,jproj)
1649 : end do priloop
1650 : end do prjloop
1651 :
1652 37 : if(numproj/=0)then
1653 :
1654 100 : ABI_MALLOC(uu,(numproj,numproj))
1655 100 : ABI_MALLOC(zz,(2,numproj,numproj))
1656 :
1657 25 : if (numproj > 1) then
1658 3 : call ZHPEV(jobz,uplo,numproj,ap,ww,zz,numproj,work1,rwork1,info)
1659 27 : uu(:,:)=zz(1,:,:)
1660 : else
1661 22 : ww(1)=hij(ll,1,1)
1662 22 : uu(1,1)=one
1663 : end if
1664 :
1665 : ! Initialization of ekb, and spline fitting
1666 :
1667 25 : if (ll==0) then ! s channel
1668 :
1669 12 : rrl=rr(0)
1670 27 : do iproj=1,numproj
1671 15 : ekb(1,iproj)=ww(iproj)*32.d0*(rrl**3)*(pi**2.5d0)/(4.d0*pi)**2
1672 15 : if(iproj==1)then
1673 36024 : do iqgrid=1,mqgrid
1674 36024 : ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2)
1675 : end do
1676 12 : yp1j(1)=zero
1677 12 : ypnj(1)=-(two_pi*rrl)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrl)**2)
1678 3 : else if(iproj==2)then
1679 6004 : do iqgrid=1,mqgrid
1680 : ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) &
1681 : & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) &
1682 6004 : & *( 3.d0-(two_pi*qgrid(iqgrid)*rrl)**2 )
1683 : end do
1684 2 : yp1j(2)=zero
1685 : ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrl)**2*qmax &
1686 2 : & *exp(-0.5d0*(two_pi*qmax*rrl)**2) * (-5.d0+(two_pi*qmax*rrl)**2)
1687 1 : else if(iproj==3)then
1688 3002 : do iqgrid=1,mqgrid
1689 : ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
1690 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1691 : & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrl)**2 + &
1692 3002 : & (two_pi*qgrid(iqgrid)*rrl)**4)
1693 : end do
1694 1 : yp1j(3)=zero
1695 : ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2) * &
1696 1 : & (two_pi*rrl)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrl)**2-(two_pi*qmax*rrl)**4)
1697 : end if
1698 : call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
1699 27 : & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
1700 : end do
1701 :
1702 13 : else if (ll==1) then ! p channel
1703 :
1704 12 : rrl=rr(1)
1705 25 : do iproj=1,numproj
1706 13 : ekb(2,iproj)=ww(iproj)*64.d0*(rrl**5)*(pi**2.5d0)/(4.d0*pi)**2
1707 13 : if(iproj==1)then
1708 36024 : do iqgrid=1,mqgrid
1709 : ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
1710 36024 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))
1711 : end do
1712 12 : yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
1713 : ypnj(1)=-two_pi*((two_pi*qmax*rrl)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)*&
1714 12 : & (1.0d0/sqrt(3.0d0))
1715 1 : else if(iproj==2)then
1716 3002 : do iqgrid=1,mqgrid
1717 : ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
1718 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1719 3002 : & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrl)**2)
1720 : end do
1721 1 : yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
1722 : ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1723 1 : & (-8*(two_pi*qmax*rrl)**2 + (two_pi*qmax*rrl)**4 + 5.0d0)
1724 0 : else if(iproj==3)then
1725 0 : do iqgrid=1,mqgrid
1726 : ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
1727 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1728 : & (two_pi*qgrid(iqgrid))*&
1729 0 : & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrl)**2+(two_pi*qgrid(iqgrid)*rrl)**4)
1730 : end do
1731 0 : yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
1732 : ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1733 : & (35.0d0-77.0d0*(two_pi*qmax*rrl)**2+19.0d0*(two_pi*qmax*rrl)**4 - &
1734 0 : & (two_pi*qmax*rrl)**6)
1735 : end if
1736 : call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
1737 25 : & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
1738 : end do
1739 :
1740 1 : else if (ll==2) then ! d channel
1741 :
1742 : ! If there is a third projector. Warning : only two projectors are allowed.
1743 1 : if ( numproj>2 ) then
1744 : write(msg, '(3a)' )&
1745 0 : & ' only two d-projectors are allowed ',ch10,&
1746 0 : & ' Action: check your pseudopotential file.'
1747 0 : ABI_ERROR(msg)
1748 : end if
1749 :
1750 1 : rrl=rr(2)
1751 2 : do iproj=1,numproj
1752 1 : ekb(3,iproj)=ww(iproj)*128.d0*(rrl**7)*(pi**2.5d0)/(4.d0*pi)**2
1753 1 : if(iproj==1)then
1754 3002 : do iqgrid=1,mqgrid
1755 : ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
1756 3002 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))**2
1757 : end do
1758 1 : yp1j(1)=zero
1759 : ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
1760 1 : & exp(-0.5d0*(two_pi*qmax*rrl)**2)*qmax*(2d0-(two_pi*qmax*rrl)**2)
1761 0 : else if(iproj==2)then
1762 0 : do iqgrid=1,mqgrid
1763 : ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
1764 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1765 0 : & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrl)**2)
1766 : end do
1767 0 : yp1j(2)=zero
1768 : ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1769 0 : & qmax*(two_pi**2)*( (two_pi*qmax*rrl)**4 - 11.0d0*(two_pi*qmax*rrl)**2 + 14.0d0)
1770 : end if
1771 : call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
1772 2 : & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
1773 : end do
1774 :
1775 0 : else if (ll==3) then ! f channel
1776 :
1777 : ! If there is a second projector. Warning : only one projector is allowed.
1778 0 : if ( numproj>1 ) then
1779 : write(msg, '(a,a,a)' )&
1780 0 : & 'only one f-projector is allowed ',ch10,&
1781 0 : & 'Action: check your pseudopotential file.'
1782 0 : ABI_ERROR(msg)
1783 : end if
1784 :
1785 0 : rrl=rr(3)
1786 0 : ekb(4,1)=ww(1)*(256.0d0/105.0d0)*(rrl**9)*(pi**2.5d0)/(4.d0*pi)**2
1787 0 : do iqgrid=1,mqgrid
1788 : ppspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
1789 0 : & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2)
1790 : end do
1791 : ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
1792 0 : yp1j(1)=zero
1793 : ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrl)**2)*&
1794 0 : & (3.0d0-(two_pi*qmax*rrl)**2)
1795 : ! Fit spline to get second derivatives by spline fit
1796 : call spline(qgrid,ppspl(:,1,4,1),mqgrid,&
1797 0 : & yp1j(1),ypnj(1),ppspl(:,2,4,1))
1798 :
1799 : else
1800 0 : ABI_ERROR("lmax>3?")
1801 : end if
1802 :
1803 : ! Linear combination using the eigenvectors
1804 198190 : ffspl(:,:,ll+1,:)=zero
1805 54 : do jproj=1,numproj
1806 93 : do iproj=1,numproj
1807 117107 : do iqgrid=1,mqgrid
1808 : ffspl(iqgrid,1:2,ll+1,jproj)=ffspl(iqgrid,1:2,ll+1,jproj) &
1809 351156 : & +uu(iproj,jproj)*ppspl(iqgrid,1:2,ll+1,iproj)
1810 : end do
1811 : end do
1812 : end do
1813 :
1814 25 : ABI_FREE(uu)
1815 25 : ABI_FREE(zz)
1816 :
1817 : ! End condition on numproj(/=0)
1818 : end if
1819 :
1820 : end do lloop
1821 :
1822 12 : ABI_FREE(ppspl)
1823 12 : ABI_FREE(work)
1824 :
1825 12 : end subroutine psp10nl
1826 : !!***
1827 :
1828 : end module m_psp_hgh
1829 : !!***
|
