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