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