Branch data Line data Source code
1 : : !!****m* ABINIT/m_psp5
2 : : !! NAME
3 : : !! m_psp5
4 : : !!
5 : : !! FUNCTION
6 : : !! Initialize pspcod=5 ("Phoney pseudopotentials" with Hamman grid):
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, FrD, FJ, MT)
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_psp5
23 : :
24 : : use defs_basis
25 : : use m_splines
26 : : use m_errors
27 : : use m_abicore
28 : :
29 : : use m_psptk, only : psp1cc, psp5lo, psp5nl
30 : :
31 : : implicit none
32 : :
33 : : private
34 : : !!***
35 : :
36 : : public :: psp5in
37 : : !!***
38 : :
39 : : contains
40 : : !!***
41 : :
42 : : !!****f* ABINIT/psp5in
43 : : !! NAME
44 : : !! psp5in
45 : : !!
46 : : !! FUNCTION
47 : : !! Initialize pspcod=5 ("Phoney pseudopotentials" with Hamman grid):
48 : : !! continue to read the corresponding file, then compute the
49 : : !! local and non-local potentials.
50 : : !!
51 : : !! INPUTS
52 : : !! lloc=angular momentum choice of local pseudopotential
53 : : !! lmax=value of lmax mentioned at the second line of the psp file
54 : : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
55 : : !! =if useylm=0, max number of (l,n) comp. over all type of psps
56 : : !! lnmax=max. number of (l,n) components over all type of psps
57 : : !! mmax=maximum number of points in real space grid in the psp file
58 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
59 : : !! mpssoang= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
60 : : !! mqgrid=dimension of q (or G) grid for arrays.
61 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
62 : : !! pspso= spin orbit signal
63 : : !! qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
64 : : !! useylm=governs the way the nonlocal operator is to be applied:
65 : : !! 1=using Ylm, 0=using Legendre polynomials
66 : : !! zion=nominal valence of atom as specified in psp file
67 : : !! znucl=nuclear number of atom as specified in psp file
68 : : !!
69 : : !! OUTPUT
70 : : !! ekb(lnmax)=Kleinman-Bylander energy,
71 : : !! {{\ \begin{equation}
72 : : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
73 : : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
74 : : !! \end{equation} }}
75 : : !! for each (l,n)
76 : : !! if any, spin-orbit components begin at l=mpsang+1
77 : : !! ekb1(mpssoang)= Kleinman-Bylander energy from the psp file, for iproj=1
78 : : !! ekb2(mpssoang)= Kleinman-Bylander energy from the psp file, for iproj=2
79 : : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree)
80 : : !! epspsp(mpssoang)=values of epsatm for different angular momenta, from the psp file
81 : : !! e990(mpssoang)=ecut at which 0.99 of the kinetic energy is recovered
82 : : !! e999(mpssoang)=ecut at which 0.999 of the kinetic energy is recovered
83 : : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
84 : : !! second derivative from spline fit for each angular momentum and
85 : : !! each projector; if any, spin-orbit components begin at l=mpsang+1
86 : : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
87 : : !! or i=lmn (if useylm=1)
88 : : !! nproj(mpssoang)=number of projection functions for each angular momentum
89 : : !! qchrg is the total (integrated) core charge
90 : : !! rcpsp(mpssoang)=cut-off radius for each angular momentum
91 : : !! rms(mpssoang)=root mean square of the KB psp
92 : : !! vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit
93 : : !! xcccrc=XC core correction cutoff radius (bohr) from psp file
94 : : !! xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
95 : : !!
96 : : !! SOURCE
97 : :
98 : 145 : subroutine psp5in(ekb,ekb1,ekb2,epsatm,epspsp,e990,e999,ffspl,indlmn,&
99 : : & lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid,&
100 : 145 : & nproj,n1xccc,pspso,qchrg,qgrid,rcpsp,rms,&
101 : 145 : & useylm,vlspl,xcccrc,xccc1d,zion,znucl)
102 : :
103 : : !Arguments ------------------------------------
104 : : !scalars
105 : : integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid
106 : : integer,intent(in) :: n1xccc,pspso,useylm
107 : : real(dp),intent(in) :: zion,znucl
108 : : real(dp),intent(out) :: epsatm,qchrg,xcccrc
109 : : !arrays
110 : : integer,intent(out) :: indlmn(6,lmnmax) !vz_i
111 : : integer,intent(inout) :: nproj(mpssoang) !vz_i
112 : : real(dp),intent(in) :: qgrid(mqgrid)
113 : : real(dp),intent(out) :: e990(mpssoang),e999(mpssoang),ekb(lnmax)
114 : : real(dp),intent(out) :: ekb1(mpssoang),ekb2(mpssoang),epspsp(mpssoang)
115 : : real(dp),intent(out) :: rcpsp(mpssoang),rms(mpssoang) !vz_i
116 : : real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
117 : : real(dp),intent(out) :: vlspl(mqgrid,2) !vz_i
118 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
119 : :
120 : : !Local variables-------------------------------
121 : : !scalars
122 : : integer :: i1,i2,ii,iln,index,ipsang,kk,lhigh,ll,mm,mproj,nn,nso,pspso0
123 : : real(dp) :: al,fchrg,r1,rchrg,yp1,ypn
124 : : logical :: test
125 : : character(len=500) :: message,errmsg
126 : : !arrays
127 : 145 : real(dp),allocatable :: ekb_so(:),ekb_sr(:),ekb_tmp(:,:),ffspl_so(:,:,:)
128 : 145 : real(dp),allocatable :: ffspl_sr(:,:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:)
129 : 145 : real(dp),allocatable :: vpspll(:,:),vpspll_so(:,:),wfll(:,:),wfll_so(:,:)
130 : 145 : real(dp),allocatable :: work_space(:),work_spl(:)
131 : :
132 : : ! ***************************************************************************
133 : :
134 : : !File format of formatted Phoney psp input (the 3 first lines
135 : : !have already been read in calling -pspatm- routine) :
136 : :
137 : : !(1) title (character) line
138 : : !(2) znucl,zion,pspdat
139 : : !(3) pspcod,pspxc,lmax,lloc,mmax,r2well
140 : : !(4) r1,al,pspso
141 : : !For each angular momentum :
142 : : !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll)
143 : : !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll)
144 : : !(6) rchrg,fchrg,qchrg
145 : : !(7) ll
146 : : !(8) (vpsp(j,ll),j=0,nmax)
147 : : !Then for iproj=1 to 2
148 : : !for ll=0,lmax
149 : : !(10) ll
150 : : !(11) ((upsp(j,ll,iproj),j=0,nmax)
151 : :
152 : : !Read fourth line of the file ; parameter pspso is optional (but
153 : : !this is not treated correctly by all machines - problems with SGI )
154 : 145 : pspso0=1
155 [ + - ]: 145 : read (tmp_unit,fmt=*,err=50,end=50) r1,al,pspso0
156 : : 50 continue
157 : :
158 [ - + ]: 145 : if(pspso0/=1 .and. pspso0/=2)then
159 : : write(message, '(3a,i0,2a)' )&
160 : 0 : & 'Problem reading the fourth line of pseudopotential file.',ch10,&
161 : 0 : & 'The parameter pspso should be 1 or 2, but it is pspso= ',pspso0,ch10,&
162 : 0 : & 'Action: check your pseudopotential input file.'
163 : 0 : ABI_ERROR(message)
164 : : end if
165 : :
166 : 145 : write(message, '(2es16.6,t47,a)' ) r1,al,'r1 and al (Hamman grid)'
167 : 145 : call wrtout(ab_out,message,'COLL')
168 : 145 : call wrtout(std_out, message,'COLL')
169 : :
170 [ + + ]: 145 : if (pspso0/=1) then
171 : 1 : write(message,'(a)') ' Pseudopotential is in spin-orbit format '
172 : 1 : call wrtout(ab_out,message,'COLL')
173 : 1 : call wrtout(std_out, message,'COLL')
174 : : end if
175 : :
176 [ + + - + ]: 145 : if (pspso/=0.and.pspso0==1) then
177 : : write(message, '(a,a,a,a,a)' )&
178 : 0 : & 'The treatment of spin-orbit interaction is required (pspso/=0)',ch10,&
179 : 0 : & 'but pseudopotential file format cannot contain spin-orbit information !',ch10,&
180 : 0 : & 'Action: check your pseudopotential input file.'
181 : 0 : ABI_ERROR(message)
182 : : end if
183 : :
184 [ + + ]: 145 : nso=1;if (pspso0/=1) nso=2
185 [ + + ]: 559 : do ipsang=1,(nso*lmax)+1
186 [ + - ]: 414 : read (tmp_unit,*, err=10, iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang)
187 : : write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) &
188 : 414 : & ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),'l,e99.0,e99.9,nproj,rcpsp'
189 : 414 : call wrtout(ab_out,message,'COLL')
190 : 414 : call wrtout(std_out, message,'COLL')
191 [ + - ]: 414 : read (tmp_unit,*, err=10, iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang)
192 : : write(message, '(4f13.8,t55,a)' ) &
193 : 414 : & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),' rms, ekb1, ekb2, epsatm'
194 : 414 : call wrtout(ab_out,message,'COLL')
195 : 559 : call wrtout(std_out, message,'COLL')
196 : : end do
197 : :
198 : : !If pspso/=0 and nproj/=2, forces nproj to be 2
199 : : !(for compatibility with old psp-file format)
200 [ + + + - ]: 145 : if (pspso/=0.and.lmax>0) then
201 : : test=.false.
202 [ + + ]: 6 : do ipsang=1,(nso*lmax)+1
203 [ + + - + ]: 6 : if (ipsang>1.and.nproj(ipsang)/=2) then
204 : 0 : test=.true.;nproj(ipsang)=2
205 : : end if
206 : : end do
207 [ - + ]: 1 : if (test) then
208 : : write(message, '(a,a,a,a,a)' )&
209 : 0 : & 'Pseudopotential file is spin-orbit (pspso=2)',ch10,&
210 : 0 : & 'and number of projector for l/=0 is not 2 !',ch10,&
211 : 0 : & 'It has been forced to 2.'
212 : 0 : call wrtout(std_out,message,'COLL')
213 : 0 : ABI_WARNING(message)
214 : : end if
215 : : end if
216 : :
217 : : !mproj=maxval(nproj(1:lmax+1))
218 : : !mjv 10/2008: I believe this is correct. Perhaps unnecessary if the normal
219 : : !projectors are always more numerous, but should be conservative anyway with
220 : : !maxval.
221 [ + + ]: 575 : mproj=maxval(nproj)
222 [ + + + + ]: 2896 : index=0;iln=0;indlmn(:,:)=0
223 [ + + ]: 291 : do nn=1,nso
224 [ + + ]: 705 : do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
225 [ + + ]: 560 : if (nproj(ipsang)>0) then
226 : 271 : ll=ipsang-(nn-1)*lmax-1
227 [ + + ]: 578 : do kk=1,nproj(ipsang)
228 : 307 : iln=iln+1
229 [ + + ]: 893 : do mm=1,2*ll*useylm+1
230 : 315 : index=index+1
231 : 315 : indlmn(1,index)=ll
232 : 315 : indlmn(2,index)=mm-ll*useylm-1
233 : 315 : indlmn(3,index)=kk
234 : 315 : indlmn(4,index)=ll*ll+(1-useylm)*ll+mm
235 : 315 : indlmn(5,index)=iln
236 : 622 : indlmn(6,index)=nn
237 : : end do
238 : : end do
239 : : end if
240 : : end do
241 : : end do
242 : :
243 [ + - ]: 145 : read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
244 : 145 : write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
245 : 145 : call wrtout(ab_out,message,'COLL')
246 : 145 : call wrtout(std_out, message,'COLL')
247 : :
248 : : !Generate core charge function and derivatives, if needed
249 : 145 : xcccrc=zero
250 [ + + ]: 145 : if(n1xccc>0)then
251 : : ! Use the revised expression of 5 Nov 1992, also used for format=1.
252 : 1 : call psp1cc(fchrg,n1xccc,xccc1d)
253 : 1 : xcccrc=3*rchrg
254 : : end if
255 : :
256 : : !--------------------------------------------------------------------
257 : : !Will now proceed at the reading of pots and wfs, as well as their treatment
258 : :
259 : : !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials
260 : : !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg.
261 : : !rad(:)=radial grid r(i)
262 : : !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions
263 [ - + + - : 435 : ABI_MALLOC(vloc,(mmax))
+ - + - -
+ - + ]
264 [ - + + - : 580 : ABI_MALLOC(vpspll,(mmax,mpsang))
- + + - +
- + - - +
- + ]
265 : :
266 : : !(1) Read atomic pseudopotential for each l, filling up array vpspll
267 : : !Note: put each l into vpspll(:,l+1)
268 : :
269 [ + + ]: 145 : if (pspso==0) then
270 : :
271 : : ! --NON SPIN-ORBIT
272 [ + + ]: 553 : do ipsang=1,lmax+1
273 [ + - ]: 409 : read (tmp_unit,*, err=10, iomsg=errmsg) ll
274 [ + - ]: 553 : read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax)
275 : : ! write(std_out,*) 'END OF READING PSP',ll,'OK'
276 : : end do
277 : : else
278 : :
279 : : ! --SPIN-ORBIT
280 [ - + + - : 3 : ABI_MALLOC(vpspll_so,(mmax,mpsang))
- + + - -
+ - + ]
281 [ + - ]: 1 : read (tmp_unit,*, err=10, iomsg=errmsg) ll
282 [ + - ]: 1 : read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll(ii,1),ii=1,mmax)
283 [ + + ]: 1047 : vpspll_so(:,1)=0.0d0
284 [ + + ]: 3 : do ipsang=2,lmax+1
285 [ + - ]: 2 : read (tmp_unit,*, err=10, iomsg=errmsg) ll
286 [ + - ]: 2 : read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax)
287 [ + - ]: 2 : read (tmp_unit,*, err=10, iomsg=errmsg) ll
288 [ + - ]: 3 : read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll_so(ii,ipsang),ii=1,mmax)
289 : : end do
290 : : end if
291 : :
292 : : !Copy appropriate nonlocal psp for use as local one
293 [ + + ]: 145 : if (pspso==0) then
294 [ + + ]: 62216 : vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 )
295 : : else
296 [ + - ]: 1 : if(lloc<=0) then
297 [ + + ]: 1047 : vloc( 1:mmax ) = vpspll( 1:mmax , -lloc+1 )
298 : : else
299 [ # # ]: 0 : vloc( 1:mmax ) = vpspll_so( 1:mmax , lloc+1 )
300 : : end if
301 : : end if
302 : : !DEBUG
303 : : !write(std_out,*) 'VLOC=',vloc(1),vloc(2),vloc(3)
304 : : !write(std_out,*) 'VLOC=',vloc(4),vloc(5),vloc(6)
305 : : !ENDDEBUG
306 : :
307 : :
308 : : !(2) Create radial grid, and associated quantities
309 : :
310 : : !Now compute Hamman Grid
311 [ - + + - : 290 : ABI_MALLOC(rad,(mmax))
- + - + ]
312 [ + + ]: 63263 : do ii=1,mmax
313 : 63263 : rad (ii)=r1*exp(dble(ii-1)*al)
314 : : end do
315 : : !DEBUG
316 : : !write(std_out,*) 'HAMMAN RADIAL GRID r1 and al',r1,al
317 : : !write(std_out,*) 'rad(1)=',rad(1)
318 : : !write(std_out,*) 'rad(10)=',rad(10)
319 : : !write(std_out,*) 'rad(100)=',rad(100)
320 : : !ENDDEBUG
321 : :
322 : :
323 : : !(3)Carry out calculations for local (lloc) pseudopotential.
324 : : !Obtain Fourier transform (1-d sine transform)
325 : : !to get q^2 V(q).
326 : :
327 : : call psp5lo(al,epsatm,mmax,mqgrid,qgrid,&
328 : 145 : & vlspl(:,1),rad,vloc,yp1,ypn,zion)
329 : :
330 : :
331 : : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
332 [ - + + - : 435 : ABI_MALLOC(work_space,(mqgrid))
+ - + - -
+ - + ]
333 [ - + + - : 290 : ABI_MALLOC(work_spl,(mqgrid))
- + - + ]
334 : 145 : call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
335 [ + + ]: 435656 : vlspl(:,2)=work_spl(:)
336 : :
337 [ - + ]: 145 : ABI_FREE(work_space)
338 [ - + ]: 145 : ABI_FREE(work_spl)
339 : :
340 : : !(4)Take care of non-local part
341 : :
342 : : !DEBUG
343 : : !write(std_out,*)' psp5in : before nonlocal corrections '
344 : : !write(std_out,*)' psp5in : lloc, lmax = ',lloc,lmax
345 : : !ENDDEBUG
346 : :
347 : : !Zero out all Kleinman-Bylander energies to initialize
348 [ + + ]: 530 : ekb(:)=0.0d0
349 : :
350 : : !Allow for option of no nonlocal corrections (lloc=lmax=0)
351 [ + + + + ]: 145 : if (lloc==0.and.lmax==0) then
352 : :
353 : 11 : write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
354 : 11 : call wrtout(ab_out,message,'COLL')
355 : 11 : call wrtout(std_out, message,'COLL')
356 : :
357 : : else
358 : :
359 : : ! Proceed to make Kleinman-Bylander form factors for
360 : : ! each l up to lmax
361 : :
362 : : ! Read wavefunctions for each l up to lmax
363 [ - + + - : 402 : ABI_MALLOC( wfll,(mmax,mpsang))
- + + - -
+ - + ]
364 : : ! -----------------------------------------------------------------
365 : :
366 [ + + ]: 134 : if (pspso==0) then
367 : :
368 : : ! --NON SPIN-ORBIT
369 [ + + ]: 531 : do ipsang=1,lmax+1
370 [ + + ]: 531 : if (nproj(ipsang)/=0) then
371 [ + - ]: 266 : read (tmp_unit,*, err=10, iomsg=errmsg) ll
372 [ - + ]: 266 : if (ipsang/=ll+1) then
373 : : write(message, '(a,a,a,a,a,a,2i6,a,a)' )&
374 : 0 : & 'Pseudopotential input file does not have',ch10,&
375 : 0 : & 'angular momenta in order expected for first projection',&
376 : 0 : & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,&
377 : 0 : & 'Action: check your pseudopotential input file.'
378 : 0 : ABI_ERROR(message)
379 : : end if
380 [ + - ]: 266 : read (tmp_unit,*, err=10, iomsg=errmsg) wfll(:,ipsang)
381 : : else
382 [ + + ]: 53704 : wfll(:,ipsang)=0.0d0
383 : : end if
384 : : end do
385 : : else
386 : :
387 : : ! --SPIN-ORBIT
388 [ - + + - : 3 : ABI_MALLOC(wfll_so,(mmax,mpsang))
- + + - -
+ - + ]
389 [ + - ]: 1 : if (nproj(1)/=0) then
390 [ + - ]: 1 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
391 [ + - ]: 1 : read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,1)
392 : : else
393 [ # # ]: 0 : wfll(:,1)=0.0d0
394 : : end if
395 [ + + ]: 1047 : wfll_so(:,1)=0.0d0
396 [ + + ]: 3 : do ipsang=2,lmax+1
397 [ + - ]: 3 : if (nproj(ipsang)/=0) then
398 [ + - ]: 2 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
399 [ + - ]: 2 : read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang)
400 [ + - ]: 2 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
401 [ + - ]: 2 : read (tmp_unit,*,err=10,iomsg=errmsg) wfll_so(:,ipsang)
402 : : else
403 [ # # ]: 0 : wfll(:,ipsang)=0.0d0
404 [ # # ]: 0 : wfll_so(:,ipsang)=0.0d0
405 : : end if
406 : : end do
407 : : end if
408 : :
409 : : ! ----------------------------------------------------------------------
410 : : ! Compute KB form factors and fit splines
411 [ - + + - : 536 : ABI_MALLOC(ekb_tmp,(mpssoang,max(nso,mproj)))
- + + - +
- + - - +
- + ]
412 [ - + + - : 670 : ABI_MALLOC(ffspl_tmp,(mqgrid,2,mpssoang,max(nso,mproj)))
- + + - -
+ + - + -
+ - - + -
+ ]
413 [ + + + + ]: 747 : ekb_tmp(:,:)=0.d0
414 : :
415 [ - + + - : 402 : ABI_MALLOC(ekb_sr,(mpsang))
+ - + - -
+ - + ]
416 [ - + + - : 536 : ABI_MALLOC(ffspl_sr,(mqgrid,2,mpsang))
- + + - +
- + - - +
- + ]
417 : : call psp5nl(al,ekb_sr(:),ffspl_sr(:,:,:),lmax,mmax,mpsang,mqgrid,&
418 : 134 : & qgrid,rad,vloc,vpspll,wfll)
419 [ + + ]: 537 : ekb_tmp(1:mpsang,1)=ekb_sr(1:mpsang)
420 [ + + + + : 2422345 : ffspl_tmp(:,:,1:mpsang,1)=ffspl_sr(:,:,1:mpsang)
+ + ]
421 : :
422 [ + + ]: 134 : if (pspso/=0) then
423 [ - + + - : 2 : ABI_MALLOC(ekb_so,(mpsang))
- + - + ]
424 [ - + + - : 3 : ABI_MALLOC(ffspl_so,(mqgrid,2,mpsang))
- + + - -
+ - + ]
425 : : call psp5nl(al,ekb_so,ffspl_so,lmax,mmax,mpsang,mqgrid,&
426 : 1 : & qgrid,rad,vloc,vpspll_so,wfll_so)
427 [ + + ]: 3 : ekb_tmp(mpsang+1:mpssoang,1)=ekb_so(2:mpsang)
428 [ + + ]: 3 : do ipsang=2,lmax+1
429 [ - + ]: 3 : if((ekb_sr(ipsang)*ekb_so(ipsang))<0.0) then
430 : 0 : ABI_ERROR('BIG PROBLEM WITH THE SPIN ORBIT IN PSP5NL')
431 : : end if
432 : : end do
433 : :
434 [ - + ]: 1 : if(lloc<0) ekb_sr(-lloc+1)=ekb_so(-lloc+1)
435 [ - + ]: 1 : if(lloc<0) ekb_tmp(-lloc+1,1)=ekb_tmp(-lloc+1+lmax,1)
436 [ - + ]: 1 : if(lloc>0) ekb_so(lloc+1)=ekb_sr(lloc+1)
437 [ - + ]: 1 : if(lloc>0) ekb_tmp(lmax+lloc+1,1)=ekb_tmp(lloc+1,1)
438 [ + + ]: 6 : do ipsang=1,mpssoang
439 [ + + ]: 5 : if(ekb_tmp(ipsang,1)>0) ekb_tmp(ipsang,1)= 1.d0
440 [ + + ]: 6 : if(ekb_tmp(ipsang,1)<0) ekb_tmp(ipsang,1)=-1.d0
441 : : end do
442 : :
443 : : ! v_ion is calculated in ffspl_tmp(:,:,1:mpsang,1) and v_so in
444 : : ! ffspl_tmp(:,:,mpsang+1:mpssoang,1) taking into account sqrt(ekb)
445 [ + + ]: 3002 : do i1=1,mqgrid
446 [ + + ]: 9004 : do i2=1,2
447 : 6002 : ffspl_tmp(i1,i2,1,1)=ffspl_sr(i1,i2,1)*sqrt(abs(ekb_sr(1)))
448 [ + + ]: 21007 : do ipsang=2,mpsang
449 : : ffspl_tmp(i1,i2,ipsang,1)=((ffspl_sr(i1,i2,ipsang)*&
450 : : & sqrt(abs(ekb_sr(ipsang)))*(ipsang-1))+&
451 : : & (ffspl_so(i1,i2,ipsang)*&
452 : : & sqrt(abs(ekb_so(ipsang)))*(ipsang)))&
453 : 12004 : & /(2.d0*ipsang-1)
454 : : ffspl_tmp(i1,i2,mpsang+ipsang-1,1)=(-ffspl_sr(i1,i2,ipsang)*&
455 : : & sqrt(abs(ekb_sr(ipsang)))+&
456 : : & ffspl_so(i1,i2,ipsang)*&
457 : : & sqrt(abs(ekb_so(ipsang))))*2.d0&
458 : 18006 : & /(2.d0*ipsang-1)
459 : : end do
460 : : end do
461 : : end do
462 [ - + ]: 1 : ABI_FREE(ekb_so)
463 [ - + ]: 1 : ABI_FREE(ffspl_so)
464 [ - + ]: 1 : ABI_FREE(vpspll_so)
465 [ - + ]: 1 : ABI_FREE(wfll_so)
466 : :
467 : : ! The non local contribution is written as quadratic form of the vector
468 : : ! V=(v_ion,v_so)
469 : : ! t_V (Q1+Q2 L.S) V
470 : : ! with Q1= (1 0 ) et Q2=(0 1 )
471 : : ! (0 l(l+1)/4) (1 -1/2)
472 : : ! The LS independent part is already diagonal. V is therefore built
473 : : ! putting v_so in the second projector of ffspl for the non spin-orbit
474 : : ! part and taking the eigenvalues of Q1 as new ekb (apart the sign)
475 [ + + ]: 3 : do ipsang=2,mpsang
476 [ + + ]: 6004 : do i1=1,mqgrid
477 [ + + ]: 18008 : do i2=1,2
478 : 18006 : ffspl_tmp(i1,i2,ipsang,2)= ffspl_tmp(i1,i2,mpsang+ipsang-1,1)
479 : : end do
480 : : end do
481 : 3 : ekb_tmp(ipsang,2)=ekb_tmp(mpsang+ipsang-1,1)*ipsang*(ipsang-1)*0.25d0
482 : : end do
483 : :
484 : : ! For the spin orbit part, after diagonalisation of Q2, the eigenvectors
485 : : ! are: ((1-sqrt(17))/4 , 1) and ((1+sqrt(17))/4 ,1)
486 : : ! The passage matrix is therefore P=((1-sqrt(17))/4 (1+sqrt(17))/4)
487 : : ! ( 1 1 )
488 : : ! t_P*Q2*P=( -sqrt(17)/2 0 )
489 : : ! ( 0 sqrt(17)/2)
490 : : ! The diagonal values are the new ekb and the new ffspl are
491 : : ! P^-1 (v_ion)
492 : : ! (v_so )
493 [ + + ]: 3 : do ipsang=2,mpsang
494 [ + + ]: 6004 : do i1=1,mqgrid
495 [ + + ]: 18008 : do i2=1,2
496 : : ffspl_tmp(i1,i2,mpsang+ipsang-1,1)=-2.d0/sqrt(17.d0)*&
497 : : & (ffspl_tmp(i1,i2,ipsang,1)-&
498 : : & ((sqrt(17.d0)+1)*0.25d0)*&
499 : 12004 : ffspl_tmp(i1,i2,ipsang,2))
500 : : ffspl_tmp(i1,i2,mpsang+ipsang-1,2)=2.d0/sqrt(17.d0)*&
501 : : & (ffspl_tmp(i1,i2,ipsang,1)+&
502 : : & ((sqrt(17.d0)-1)*0.25d0)*&
503 : 18006 : & ffspl_tmp(i1,i2,ipsang,2))
504 : : end do
505 : : end do
506 : 2 : ekb_tmp(mpsang+ipsang-1,1)=-(sqrt(17.d0)*0.5d0)*ekb_tmp(ipsang,1)
507 : 3 : ekb_tmp(mpsang+ipsang-1,2)= (sqrt(17.d0)*0.5d0)*ekb_tmp(ipsang,1)
508 : : end do
509 : :
510 : : end if
511 : :
512 [ - + ]: 134 : ABI_FREE(ekb_sr)
513 [ - + ]: 134 : ABI_FREE(ffspl_sr)
514 : :
515 : : ! FJ WARNING : No spin orbit if nproj>1
516 [ + + ]: 134 : if (pspso==0) then
517 : :
518 : : ! Read second wavefunction for second projection operator
519 : : ! (only read cases where nproj(ll)=2)
520 : : ! --also find highest l for which nproj(l)=2
521 : 133 : lhigh=-1
522 [ + + ]: 531 : do ipsang=1,min(lmax+1,mpsang)
523 [ + + ]: 531 : if (nproj(ipsang)==2) then
524 : 32 : lhigh=ipsang-1
525 [ + - ]: 32 : read (tmp_unit,*, err=10, iomsg=errmsg) ll
526 [ - + ]: 32 : if (ipsang/=ll+1) then
527 : : write(message, '(a,a,a,a,a,a,2i6,a,a)' )&
528 : 0 : & 'Pseudopotential input file does not have',ch10,&
529 : 0 : & 'angular momenta in order expected for second projection',&
530 : 0 : & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,&
531 : 0 : & 'Action: check your pseudopotential input file.'
532 : 0 : ABI_ERROR(message)
533 : : end if
534 [ + - ]: 32 : read (tmp_unit,*, err=10, iomsg=errmsg) wfll(:,ipsang)
535 : : ! DEBUG
536 : : ! write(std_out,*) 'WF second',ipsang-1,wfll(1,ipsang),wfll(2,ipsang),wfll(3,ipsang)
537 : : ! ENDDEBUG
538 : : else
539 [ + + ]: 143882 : wfll(:,ipsang)=0.0d0
540 : : end if
541 : :
542 : : end do
543 : :
544 : : ! Compute KB form factors and fit splines for second wf if any
545 [ + + ]: 133 : if (lhigh>-1) then
546 : : call psp5nl(al,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lmax,&
547 : 17 : & mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll)
548 : : end if
549 : :
550 : : end if
551 : :
552 : : ! Convert ekb and ffspl
553 : : iln=0
554 [ + + ]: 509 : do ii=1,lmnmax
555 : 375 : kk=indlmn(5,ii)
556 [ + + ]: 509 : if (kk>iln) then
557 : 307 : iln=kk
558 : 307 : ll=indlmn(1,ii);nn=indlmn(3,ii)
559 [ + + ]: 307 : if (indlmn(6,ii)==1) then
560 : 303 : ekb(kk)=ekb_tmp(1+ll,nn)
561 [ + + + + ]: 1820979 : ffspl(:,:,kk)=ffspl_tmp(:,:,1+ll,nn)
562 : : else
563 : 4 : ekb(kk)=ekb_tmp(mpsang+ll,nn)
564 [ + + + + ]: 24020 : ffspl(:,:,kk)=ffspl_tmp(:,:,mpsang+ll,nn)
565 : : end if
566 : : end if
567 : : end do
568 : :
569 [ - + ]: 134 : ABI_FREE(ekb_tmp)
570 [ - + ]: 134 : ABI_FREE(ffspl_tmp)
571 [ - + ]: 134 : ABI_FREE(wfll)
572 : :
573 : : ! end of if concerning lloc
574 : : end if
575 : :
576 [ - + ]: 145 : ABI_FREE(vpspll)
577 [ - + ]: 145 : ABI_FREE(rad)
578 [ - + ]: 145 : ABI_FREE(vloc)
579 : :
580 : 145 : return
581 : :
582 : : ! Handle IO error
583 : : 10 continue
584 : 0 : ABI_ERROR(errmsg)
585 : :
586 [ - + - + : 145 : end subroutine psp5in
- + - + -
+ - + - +
- + ]
587 : : !!***
588 : :
589 : : end module m_psp5
590 : : !!***
|
