Line data Source code
1 : !!****m* ABINIT/m_psp1
2 : !! NAME
3 : !! m_psp1
4 : !!
5 : !! FUNCTION
6 : !! Initialize pspcod=1 or 4 pseudopotential (Teter format)
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, FrD, 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_psp1
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_splines
28 :
29 : use m_special_funcs, only : besjm
30 : use m_psptk, only : psp1cc
31 :
32 : implicit none
33 :
34 : private
35 : !!***
36 :
37 : public :: psp1in ! Initialize pspcod=1 or 4 pseudopotential (Teter format)
38 : !!***
39 :
40 : contains
41 : !!***
42 :
43 : !!****f* m_psp1/psp1in
44 : !! NAME
45 : !! psp1in
46 : !!
47 : !! FUNCTION
48 : !! Initialize pspcod=1 or 4 pseudopotential (Teter format):
49 : !! continue to read the corresponding file, then compute the
50 : !! local and non-local potentials.
51 : !!
52 : !! INPUTS
53 : !! dq= spacing of the q-grid
54 : !! lloc=angular momentum choice of local pseudopotential
55 : !! lmax=value of lmax mentioned at the second line of the psp file
56 : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
57 : !! =if useylm=0, max number of (l,n) comp. over all type of psps
58 : !! lnmax=max. number of (l,n) components over all type of psps
59 : !! mmax=maximum number of points in real space grid in the psp file
60 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
61 : !! mqgrid=dimension of q (or G) grid for arrays.
62 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
63 : !! pspcod=pseudopotential type
64 : !! qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
65 : !! useylm=governs the way the nonlocal operator is to be applied:
66 : !! 1=using Ylm, 0=using Legendre polynomials
67 : !! zion=nominal valence of atom as specified in psp file
68 : !! znucl=atomic number of atom as specified in psp file
69 : !!
70 : !! OUTPUT
71 : !! ekb(lnmax)=Kleinman-Bylander energy,
72 : !! {{\ \begin{equation}
73 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
74 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
75 : !! \end{equation} }}
76 : !! for each (l,n)
77 : !! ekb1(mpsang)= Kleinman-Bylander energy from the psp file, for iproj=1
78 : !! ekb2(mpsang)= 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(mpsang)=values of epsatm for different angular momenta, from the psp file
81 : !! e990(mpsang)=ecut at which 0.99 of the kinetic energy is recovered
82 : !! e999(mpsang)=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
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(mpsang)=number of projection functions for each angular momentum
89 : !! qchrg is the total (integrated) core charge
90 : !! rcpsp(mpsang)=cut-off radius for each angular momentum
91 : !! rms(mpsang)=root mean square of the KB psp
92 : !! vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit
93 : !! xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
94 : !! xcccrc=XC core correction cutoff radius (bohr)
95 : !!
96 : !! NOTES
97 : !! there are only minor differences in the two formats
98 : !! 1) With pspcod=1, even for the LOCAL angular momentum, there is
99 : !! a block for the wfs (can be set to zero, though)
100 : !! 2) The core charge density differs: for pspcod=1, it is a
101 : !! revised expression for core density of 5 Nov 1992, while
102 : !! for pspcod=4, it is an older expression, of 7 May 1992 .
103 : !!
104 : !! SOURCE
105 :
106 694 : subroutine psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,&
107 694 : & e990,e999,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
108 694 : & mmax,mpsang,mqgrid,nproj,n1xccc,pspcod,&
109 694 : & qchrg,qgrid,rcpsp,rms,useylm,vlspl,xcccrc,xccc1d,&
110 : & zion,znucl)
111 :
112 : !Arguments ------------------------------------
113 : !scalars
114 : integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc,pspcod
115 : integer,intent(in) :: useylm
116 : real(dp),intent(in) :: dq,zion,znucl
117 : real(dp),intent(out) :: epsatm,qchrg,xcccrc
118 : !arrays
119 : integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang)
120 : real(dp),intent(in) :: qgrid(mqgrid)
121 : real(dp),intent(out) :: e990(mpsang),e999(mpsang),ekb(lnmax),ekb1(mpsang)
122 : real(dp),intent(out) :: ekb2(mpsang),epspsp(mpsang)
123 : real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax)
124 : real(dp),intent(out) :: rcpsp(mpsang),rms(mpsang),vlspl(mqgrid,2)
125 : real(dp),intent(inout) :: xccc1d(n1xccc,6)
126 :
127 : !Local variables-------------------------------
128 : !scalars
129 : integer :: ii,iln,index,ipsang,kk,lhigh,ll,mm,nlmax
130 : real(dp) :: arg,dq2pi,fchrg,rchrg,xx,yp1,ypn
131 : character(len=500) :: message,errmsg
132 : !arrays
133 694 : real(dp),allocatable :: drad(:),ekb_tmp(:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:)
134 694 : real(dp),allocatable :: vpspll(:,:),wfll(:,:),wksincos(:,:,:),work_space(:)
135 694 : real(dp),allocatable :: work_spl1(:),work_spl2(:)
136 : ! ***************************************************************************
137 :
138 : !Note: Teter s grid is hard-coded at mmax=2001
139 : !mmax was read from the pseudopotential file in the calling routine
140 694 : if (mmax/=2001) then
141 : write(message, '(a,i12,a,a,a,a)' )&
142 0 : & 'Using Teter grid (pspcod=1 and 4) but mmax=',mmax,ch10,&
143 0 : & 'mmax must be 2001 for Teter grid.',ch10,&
144 0 : & 'Action: check your pseudopotential input file.'
145 0 : ABI_ERROR(message)
146 : end if
147 :
148 : !File format of formatted Teter psp input (the 3 first lines
149 : !have already been read in calling -pspatm- routine) :
150 :
151 : !(1) title (character) line
152 : !(2) znucl,zion,pspdat
153 : !(3) pspcod,pspxc,lmax,lloc,mmax,r2well
154 : !For each angular momentum :
155 : !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll)
156 : !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll)
157 : !(6) rchrg,fchrg,qchrg
158 : !(7) ll
159 : !(8) (vpsp(j,ll),j=0,nmax)
160 : !Then for iproj=1 to 2
161 : !for ll=0,lmax
162 : !(10) ll
163 : !(11) ((upsp(j,ll,iproj),j=0,nmax)
164 :
165 2353 : do ipsang=1,lmax+1
166 :
167 1659 : read (tmp_unit,*,err=10,iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang)
168 : write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) &
169 1659 : & ipsang-1,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),&
170 3318 : & 'l,e99.0,e99.9,nproj,rcpsp'
171 1659 : call wrtout(ab_out,message,'COLL')
172 1659 : call wrtout(std_out, message,'COLL')
173 :
174 1659 : read (tmp_unit,*,err=10,iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang)
175 : write(message, '(4f13.8,t55,a)' ) &
176 1659 : & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),&
177 3318 : & ' rms, ekb1, ekb2, epsatm'
178 1659 : call wrtout(ab_out,message,'COLL')
179 2353 : call wrtout(std_out, message,'COLL')
180 :
181 : end do
182 :
183 : !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn
184 9913 : index=0;iln=0;indlmn(:,:)=0
185 2353 : do ipsang=1,lmax+1
186 2353 : if(nproj(ipsang)>0)then
187 967 : ll=ipsang-1
188 2014 : do kk=1,nproj(ipsang)
189 1047 : iln=iln+1
190 3077 : do mm=1,2*ll*useylm+1
191 1063 : index=index+1
192 1063 : indlmn(1,index)=ll
193 1063 : indlmn(2,index)=mm-ll*useylm-1
194 1063 : indlmn(3,index)=kk
195 1063 : indlmn(4,index)=ll*ll+(1-useylm)*ll+mm
196 1063 : indlmn(5,index)=iln
197 2110 : indlmn(6,index)=1
198 : end do
199 : end do
200 : end if
201 : end do
202 :
203 694 : read (tmp_unit,*,err=10,iomsg=errmsg) rchrg,fchrg,qchrg
204 694 : write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
205 694 : call wrtout(ab_out,message,'COLL')
206 694 : call wrtout(std_out, message,'COLL')
207 :
208 : ! Generate core charge function and derivatives, if needed
209 694 : if(fchrg>1.0d-15)then
210 589 : if(pspcod==1)then
211 586 : call psp1cc(fchrg,n1xccc,xccc1d)
212 : ! The core charge function for pspcod=1 becomes zero beyond 3*rchrg only.
213 : ! Thus xcccrc must be set equal to 3*rchrg .
214 586 : xcccrc=3*rchrg
215 3 : else if(pspcod==4)then
216 3 : call psp4cc(fchrg,n1xccc,xccc1d)
217 : ! For pspcod=4, the core charge cut off exactly beyond rchrg
218 3 : xcccrc=rchrg
219 : end if
220 : else
221 105 : xcccrc=0.0d0
222 135789 : xccc1d(:,:)=0.0d0
223 : end if
224 :
225 : !--------------------------------------------------------------------
226 : !Will now proceed at the reading of pots and wfs, as well as their treatment
227 :
228 : !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials
229 : !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg.
230 : !rad(:)=radial grid r(i)
231 : !drad(:)= inverse of d(r(i))/d(i) for radial grid
232 : !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions
233 :
234 2082 : ABI_MALLOC(vloc,(mmax))
235 2776 : ABI_MALLOC(vpspll,(mmax,mpsang))
236 4698 : if(lmax==-1) vpspll(:,:)=zero
237 :
238 : !(1) Read atomic pseudopotential for each l, filling up array vpspll
239 : !Note: put each l into vpspll(:,l+1)
240 2353 : do ipsang=1,lmax+1
241 1659 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
242 2353 : read (tmp_unit,*,err=10,iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax)
243 : end do
244 :
245 : !Copy appropriate nonlocal psp for use as local one
246 1389388 : vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 )
247 :
248 : !(2) Create radial grid, and associated quantities
249 1388 : ABI_MALLOC(rad,(mmax))
250 1388 : ABI_MALLOC(drad,(mmax))
251 2776 : ABI_MALLOC(wksincos,(mmax,2,2))
252 :
253 : !Teter grid--need both r and dr in this case
254 1389388 : do ii=0,mmax-1
255 1388694 : xx=dble(ii)/dble(mmax-1)
256 1388694 : rad (ii+1)=100.d0*(xx+.01d0)**5-1.d-8
257 1389388 : drad(ii+1)=500.d0*(xx+.01d0)**4/dble(mmax-1)
258 : end do
259 :
260 : !here compute sin(r(:)*dq) and cos(r(:)*dq)
261 : !NOTE: also invert dr !!
262 694 : dq2pi=2.0d0*pi*dq
263 1389388 : do ii=1,mmax
264 1388694 : arg=dq2pi*rad(ii)
265 1388694 : drad(ii)=1.0d0/drad(ii)
266 1388694 : wksincos(ii,1,1)=sin(arg)
267 1389388 : wksincos(ii,2,1)=cos(arg)
268 : end do
269 :
270 : !(3)Carry out calculations for local (lloc) pseudopotential.
271 : !Obtain Fourier transform (1-d sine transform) to get q^2 V(q).
272 2082 : ABI_MALLOC(work_space,(mqgrid))
273 1388 : ABI_MALLOC(work_spl1,(mqgrid))
274 1388 : ABI_MALLOC(work_spl2,(mqgrid))
275 : call psp1lo(drad,epsatm,mmax,mqgrid,qgrid,&
276 694 : & work_spl1,rad,vloc,wksincos,yp1,ypn,zion)
277 :
278 : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
279 694 : call spline (qgrid,work_spl1,mqgrid,yp1,ypn,work_spl2)
280 2083953 : vlspl(:,1)=work_spl1(:)
281 2083953 : vlspl(:,2)=work_spl2(:)
282 :
283 694 : ABI_FREE(work_space)
284 694 : ABI_FREE(work_spl1)
285 694 : ABI_FREE(work_spl2)
286 :
287 : !(4)Take care of non-local part
288 :
289 : !Zero out all Kleinman-Bylander energies to initialize
290 1981 : ekb(:)=0.0d0
291 : !write(std_out,*)' psp1in : before nonlocal corrections '
292 : !write(std_out,*)' psp1in : lloc, lmax = ',lloc,lmax
293 :
294 : !Allow for option of no nonlocal corrections (lloc=lmax=0)
295 694 : if (lloc==0.and.lmax==0) then
296 97 : write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
297 97 : call wrtout(ab_out,message,'COLL')
298 97 : call wrtout(std_out, message,'COLL')
299 :
300 : else
301 :
302 : ! Proceed to make Kleinman-Bylander form factors for each l up to lmax
303 :
304 : ! Read wavefunctions for each l up to lmax
305 1791 : ABI_MALLOC(wfll,(mmax,mpsang))
306 2159 : do ipsang=1,lmax+1
307 : ! For pspcod==4, wfs for the local angular momentum are not written
308 2159 : if (nproj(ipsang)/=0 .or. pspcod==1) then
309 1553 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
310 1553 : if (ipsang/=ll+1) then
311 : write(message, '(a,a,a,a,a,a,2i6,a,a)' )&
312 0 : & 'Pseudopotential input file does not have',ch10,&
313 0 : & 'angular momenta in order expected for first projection',&
314 0 : & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,&
315 0 : & 'Action: check your pseudopotential input file.'
316 0 : ABI_ERROR(message)
317 : end if
318 1553 : read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang)
319 :
320 : else
321 18018 : wfll(:,ipsang)=0.0d0
322 : end if
323 :
324 : end do
325 : ! ----------------------------------------------------------------------
326 : ! Compute KB form factors and fit splines
327 :
328 : ! nlmax is highest l for which a nonlocal correction is being computed
329 597 : nlmax=lmax
330 597 : if (lloc==lmax) nlmax=lmax-1
331 : ! write(std_out,*)' psp1in : lmax,lloc=',lmax,lloc
332 1791 : ABI_MALLOC(ekb_tmp,(mpsang,2))
333 2985 : ABI_MALLOC(ffspl_tmp,(mqgrid,2,nlmax+1,2))
334 :
335 : call psp1nl(drad,ekb_tmp(:,1),ffspl_tmp(:,:,:,1),lloc,&
336 597 : & nlmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos)
337 :
338 : ! Read second wavefunction for second projection operator
339 : ! (only read cases where nproj(ll)=2) --also find highest l for which nproj(l)=2
340 597 : lhigh=-1
341 2159 : do ipsang=1,min(lmax+1,mpsang)
342 2159 : if (nproj(ipsang)==2) then
343 80 : lhigh=ipsang-1
344 80 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
345 80 : if (ipsang/=ll+1) then
346 : write(message, '(a,a,a,a,a,a,2i6,a,a)' )&
347 0 : & 'Pseudopotential input file does not have',ch10,&
348 0 : & 'angular momenta in order expected for second projection',&
349 0 : & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,&
350 0 : & 'Action: check your pseudopotential input file.'
351 0 : ABI_ERROR(message)
352 : end if
353 80 : read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang)
354 :
355 : else
356 2966964 : wfll(:,ipsang)=0.0d0
357 :
358 : end if
359 : end do
360 :
361 : ! Compute KB form factors and fit splines for second wf if any
362 :
363 597 : if (lhigh>-1) then
364 : call psp1nl(drad,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lloc,&
365 41 : & lhigh,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos)
366 : end if
367 :
368 : ! Convert ekb and ffspl
369 : iln=0
370 1811 : do ii=1,lmnmax
371 1214 : kk=indlmn(5,ii)
372 1811 : if (kk>iln) then
373 1047 : iln=kk
374 1047 : ekb(kk)=ekb_tmp(1+indlmn(1,ii),indlmn(3,ii))
375 : ! write(std_out,*)' psp1in : lmnmax,ii,indlmn(1,ii)=',lmnmax,ii,indlmn(1,ii)
376 6289495 : ffspl(:,:,kk)=ffspl_tmp(:,:,1+indlmn(1,ii),indlmn(3,ii))
377 : end if
378 : end do
379 :
380 597 : ABI_FREE(ekb_tmp)
381 597 : ABI_FREE(ffspl_tmp)
382 597 : ABI_FREE(wfll)
383 : end if
384 :
385 694 : ABI_FREE(vpspll)
386 694 : ABI_FREE(rad)
387 694 : ABI_FREE(drad)
388 694 : ABI_FREE(vloc)
389 694 : ABI_FREE(wksincos)
390 :
391 694 : return
392 :
393 : ! Handle IO error
394 : 10 continue
395 0 : ABI_ERROR(errmsg)
396 :
397 1388 : end subroutine psp1in
398 : !!***
399 :
400 : !!****f* m_psp1/psp1lo
401 : !! NAME
402 : !! psp1lo
403 : !!
404 : !! FUNCTION
405 : !! Compute sine transform to transform from v(r) to q^2 v(q)
406 : !! using subroutines related to Teter atomic structure grid.
407 : !!
408 : !! INPUTS
409 : !! drad(mmax)=inverse of r grid spacing at each point
410 : !! mmax=number of radial r grid points (Teter grid)
411 : !! mqgrid=number of grid points in q from 0 to qmax.
412 : !! qgrid(mqgrid)=q grid values (bohr**-1).
413 : !! rad(mmax)=r grid values (bohr).
414 : !! vloc(mmax)=v(r) on radial grid.
415 : !! wksincos(mmax,2,2)=contains sine and cosine of 2*pi*r(:)*dq and 2*pi*r(:)*q
416 : !! at input : wksincos(:,1,1)=sine of 2*pi*r(:)*dq
417 : !! wksincos(:,2,1)=cosine of 2*pi*r(:)*dq
418 : !! wksincos(:,:,2) is not initialized, will be used inside the routine
419 : !! zion=nominal valence charge of atom.
420 : !!
421 : !! OUTPUT
422 : !! epsatm= $4\pi \int[r^2 (v(r)+Zv/r) dr]$
423 : !! q2vq(mqgrid)=$q^2 v(q)$
424 : !! =$\displaystyle -Zv/\pi+q^2 4\pi\int(\frac{\sin(2\pi q r)}{2 \pi q r})(r^2 v(r)+r Zv)dr$.
425 : !! yp1,ypn=derivative of q^2 v(q) wrt q at q=0 and q=qmax
426 : !! (needed for spline fitter).
427 : !!
428 : !! SOURCE
429 :
430 694 : subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
431 694 : & vloc,wksincos,yp1,ypn,zion)
432 :
433 : !Arguments ------------------------------------
434 : !scalars
435 : integer,intent(in) :: mmax,mqgrid
436 : real(dp),intent(in) :: zion
437 : real(dp),intent(out) :: epsatm,yp1,ypn
438 : !arrays
439 : real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax)
440 : real(dp),intent(inout) :: wksincos(mmax,2,2)
441 : real(dp),intent(out) :: q2vq(mqgrid)
442 :
443 : !Local variables-------------------------------
444 : !scalars
445 : integer,parameter :: mma0=2001
446 : integer :: iq,ir,irmax
447 : real(dp),parameter :: scale=10.0d0
448 : real(dp) :: result,test,tpiq
449 : !arrays
450 : real(dp) :: wk(mma0),wk1(mma0),wk2(mma0)
451 :
452 : ! *************************************************************************
453 :
454 : !Do q=0 separately (compute epsatm)
455 : !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr]
456 : !Treat r=0 by itself
457 694 : wk(1)=0.0d0
458 :
459 1388694 : do ir=2,mmax
460 : ! (at large r do not want prefactor of r^2 and should see
461 : ! V(r)+Zv/r go to 0 at large r)
462 1388000 : test=vloc(ir)+zion/rad(ir)
463 : ! write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test
464 : ! In this routine, NO cut-off radius is imposed : the input
465 : ! vloc MUST be in real(dp) to obtain numerically
466 : ! accurate values. The error can be on the order of 0.001 Ha !
467 1388694 : if (abs(test)<1.0d-20) then
468 402838 : wk(ir)=0.0d0
469 : else
470 985162 : wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion)
471 : end if
472 : end do
473 : !Do integral from 0 to r(max) (disregard contrib beyond r(max)
474 : !(need numerical derivatives to do integral)
475 : !Use mmax-1 to convert to Teter s dimensioning starting at 0
476 694 : call der_int(wk,wk2,rad,drad,mmax-1,result)
477 :
478 694 : epsatm=4.d0*pi*(result)
479 : !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi
480 694 : q2vq(1)=-zion/pi
481 :
482 : !Prepare loop over q values
483 694 : irmax=mmax+1
484 1388694 : do ir=mmax,2,-1
485 1388000 : test=vloc(ir)+zion/rad(ir)
486 1388000 : wk1(ir)=test*rad(ir)
487 : ! Will ignore tail within decade of machine precision
488 1388694 : if ((scale+abs(test))==scale .and. irmax==ir+1) then
489 578595 : irmax=ir
490 : end if
491 : end do
492 : !Increase irmax a bit : this is copied from psp1nl
493 694 : irmax=irmax+4
494 694 : if(irmax>mmax)irmax=mmax
495 :
496 : !Loop over q values
497 2083259 : do iq=2,mqgrid
498 2082565 : tpiq=two_pi*qgrid(iq)
499 2082565 : call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
500 : ! set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral
501 : !$\displaystyle -Zv/\pi + q^2 4\pi \int[\frac{\sin(2\pi q r)}{2\pi q r}(r^2 v(r)+r Zv)dr]$.
502 : ! Handle r=0 separately
503 2082565 : wk(1)=0.0d0
504 2441041533 : do ir=2,irmax
505 2441041533 : wk(ir)=wksincos(ir,1,2)*wk1(ir)
506 : end do
507 : ! do integral from 0 to r(max)
508 2082565 : if(irmax>mmax-1)irmax=mmax-1
509 :
510 2082565 : call der_int(wk,wk2,rad,drad,irmax,result)
511 : ! store q^2 v(q)
512 2083259 : q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result
513 : end do
514 :
515 : !Compute derivatives of q^2 v(q) at ends of interval
516 694 : yp1=0.0d0
517 : !ypn=$\displaystyle 2\int_0^\infty (\sin (2\pi qmax r)+(2\pi qmax r)\cos (2\pi qmax r)(r V(r)+Z)dr]$
518 : !integral from r(mmax) to infinity is overkill; ignore
519 : !set up integrand
520 : !Handle r=0 separately
521 694 : wk(1)=0.0d0
522 694 : tpiq=two_pi*qgrid(mqgrid)
523 1388694 : do ir=2,mmax
524 1388000 : test=vloc(ir)+zion/rad(ir)
525 : ! Ignore contributions within decade of machine precision
526 1388694 : if ((scale+abs(test))==scale) then
527 579979 : wk(ir)=0.0d0
528 : else
529 : wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * &
530 808021 : & (rad(ir)*vloc(ir)+zion)
531 : end if
532 : end do
533 694 : call der_int(wk,wk2,rad,drad,mmax-1,result)
534 :
535 694 : ypn=2.0d0*result
536 :
537 694 : end subroutine psp1lo
538 : !!***
539 :
540 : !!****f* m_psp1/psp1nl
541 : !! NAME
542 : !! psp1nl
543 : !!
544 : !! FUNCTION
545 : !! Make Kleinman-Bylander form factors f_l(q) for each l from 0 to lmax.
546 : !! Vloc is assumed local potential.
547 : !!
548 : !! INPUTS
549 : !! dr(mmax)=inverse of grid spacing for radial grid
550 : !! lloc=angular momentum of local channel (avoid doing integrals for this l)
551 : !! lmax=maximum ang momentum for which nonlocal form factor is desired.
552 : !! mmax=number of radial grid points for atomic grid
553 : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
554 : !! mqgrid=number of grid points for q grid
555 : !! qgrid(mqgrid)=values at which form factors are returned
556 : !! rad(mmax)=radial grid values
557 : !! vloc(mmax)=local pseudopotential on radial grid
558 : !! vpspll(mmax,lmax+1)=nonlocal pseudopotentials for each l on radial grid
559 : !! wfll(mmax,lmax+1)=reference state wavefunctions on radial grid
560 : !! wksincos(mmax,2,2)=contains sine and cosine of 2*pi*r(:)*dq and 2*pi*r(:)*q
561 : !! at input : wksincos(:,1,1)=sine of 2*pi*r(:)*dq
562 : !! wksincos(:,2,1)=cosine of 2*pi*r(:)*dq
563 : !! wksincos(:,:,2) is not initialized, will be used inside the routine
564 : !!
565 : !! OUTPUT
566 : !! ekb(mpsang)=Kleinman-Bylander energy,
567 : !! {{\ \begin{equation}
568 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
569 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
570 : !! \end{equation} }}
571 : !! for each l
572 : !! ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and
573 : !! second derivative from spline fit for each angular momentum
574 : !!
575 : !! NOTES
576 : !! u_l(r) is reference state wavefunction (input as wfll);
577 : !! j_l(q) is a spherical Bessel function;
578 : !! dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
579 : !! f_l(q) =$ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
580 : !! where dvms=$\displaystyle \int_0^{rmax}[(u_l(r) dV_l(r))^2 dr]$ is the mean
581 : !! square value of the nonlocal correction for angular momentum l.
582 : !! E_KB = $\displaystyle \frac{dvms}{\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]}$.
583 : !! This is the eigenvalue of the Kleinman-Bylander operator and sets
584 : !! the energy scale of the nonlocal psp corrections.
585 : !! Bessel functions replaced by besj, which accomodates args near 0.
586 : !!
587 : !! SOURCE
588 :
589 638 : subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,&
590 638 : & qgrid,rad,vloc,vpspll,wfll,wksincos)
591 :
592 : !Arguments ------------------------------------
593 : !scalars
594 : integer,intent(in) :: lloc,lmax,mmax,mpsang,mqgrid
595 : !arrays
596 : real(dp),intent(in) :: dr(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax)
597 : real(dp),intent(in) :: vpspll(mmax,mpsang),wfll(mmax,mpsang)
598 : real(dp),intent(inout) :: wksincos(mmax,2,2)
599 : real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang)
600 :
601 : !Local variables-------------------------------
602 : !scalars
603 : integer,parameter :: dpsang=5
604 : integer :: iq,ir,irmax,lp1
605 : real(dp) :: dvwf,result,test,tpiq,yp1,ypn
606 : character(len=500) :: message
607 : !arrays
608 : real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang)
609 638 : real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:)
610 638 : real(dp),allocatable :: work_spl(:)
611 :
612 : ! *************************************************************************
613 :
614 : !Zero out Kleinman-Bylander energies ekb
615 2392 : ekb(:)=0.0d0
616 : !Zero out eta and other parameters too (so 0 s show up in output later)
617 : eta(:)=0.0d0
618 : dvms(:)=0.0d0
619 : ckb(:)=0.0d0
620 :
621 : !Allow for no nonlocal correction (lmax=-1)
622 638 : if (lmax/=-1) then
623 :
624 : ! Check that lmax is within allowed range
625 636 : if (lmax<0.or.lmax>3) then
626 : write(message, '(a,i12,a,a,a,a,a,a,a)' )&
627 0 : & 'lmax=',lmax,' is not an allowed value.',ch10,&
628 0 : & 'Allowed values are -1 for no nonlocal correction or else',ch10,&
629 0 : & '0, 1, 2, or 3 for maximum l nonlocal correction.',ch10,&
630 0 : & 'Action: check the input atomic psp data file for lmax.'
631 0 : ABI_ERROR(message)
632 : end if
633 :
634 : ! Compute normalizing integrals eta=<dV> and mean square
635 : ! nonlocal psp correction dvms=<dV^2>
636 : ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction
637 :
638 1908 : ABI_MALLOC(work1,(mmax+1))
639 1272 : ABI_MALLOC(work2,(mmax+1))
640 1908 : ABI_MALLOC(work_spl,(mqgrid))
641 1908 : ABI_MALLOC(work5,(mmax))
642 1272 : ABI_MALLOC(besjx,(mmax))
643 :
644 1766 : do lp1=1,lmax+1
645 :
646 : ! Only do the work if nonlocal correction is nonzero
647 1766 : if (lp1 /= lloc+1) then
648 :
649 : ! integrand for 0 to r(mmax)
650 2096094 : do ir=1,mmax
651 2095047 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
652 2096094 : work1(ir)=wfll(ir,lp1)*dvwf
653 : end do
654 :
655 : ! do integral
656 : ! first need derivative of function; note use of
657 : ! shifted indices to accomodate Mike Teter s choice of 0:mmax-1
658 1047 : call der_int(work1,work2,rad,dr,mmax-1,result)
659 1047 : eta(lp1)=result
660 :
661 : ! DEBUG
662 : ! write(std_out,*)' psp1nl : write eta(lp1)'
663 : ! write(std_out,*)result
664 : ! do ir=1,mmax,61
665 : ! write(std_out,*)vpspll(ir,lp1),vloc(ir),wfll(ir,lp1)
666 : ! end do
667 : ! write(std_out,*)
668 : ! do ir=1,mmax,61
669 : ! write(std_out,*)work1(ir),rad(ir),dr(ir)
670 : ! end do
671 : ! ENDDEBUG
672 :
673 2096094 : do ir=1,mmax
674 2095047 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
675 2096094 : work1(ir)=dvwf**2
676 : end do
677 1047 : call der_int(work1,work2,rad,dr,mmax-1,result)
678 :
679 1047 : dvms(lp1)=result
680 :
681 : ! If dvms is not 0 for any given angular momentum l,
682 : ! compute Xavier Gonze s definition of the Kleinman-Bylander
683 : ! energy E_KB = dvms/eta. In this case also renormalize
684 : ! the projection operator to u_KB(r)=$u_l(r) dV(r)/\sqrt{dvms}$.
685 : ! This means dvwf gets multiplied by the normalization factor
686 : ! "renorm"=$1/\sqrt{dvms}$ as seen below.
687 : ! With dvwf=dV(r)*wf(r) for wf(r)=``radial'' wf, the integrand
688 : ! for each angular momentum l is
689 : ! Bessel_l(2 $\pi$ q r) * wf(r) * dV(r) * r;
690 : ! NOTE presence of extra r in integrand.
691 :
692 1047 : if (dvms(lp1)/=0.0d0) then
693 1047 : ekb(lp1)=dvms(lp1)/eta(lp1)
694 1047 : renorm(lp1)=1.0d0/sqrt(dvms(lp1))
695 : ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze)
696 : ckb(lp1)=eta(lp1)/sqrt(dvms(lp1))
697 : else
698 0 : ekb(lp1)=0.0d0
699 : end if
700 : end if
701 : end do
702 :
703 : ! Loop on angular momenta
704 1766 : do lp1=1,lmax+1
705 :
706 : ! Compute form factor if ekb(lp1) not 0
707 1766 : if (ekb(lp1)/=0.0d0) then
708 :
709 : ! do q=0 separately, non-zero if l=0
710 1047 : if(lp1==1)then
711 1107106 : do ir=1,mmax
712 1106553 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
713 1107106 : work1(ir)=rad(ir)*dvwf
714 : end do
715 553 : call der_int(work1,work2,rad,dr,mmax-1,result)
716 553 : ffspl(1,1,lp1)=result
717 : else
718 : ! For l non-zero, f(q=0) vanishes !
719 494 : ffspl(1,1,lp1)=0.0d0
720 : end if
721 :
722 : ! Prepare loop over q values
723 1047 : irmax=mmax+1
724 2095047 : do ir=mmax,2,-1
725 2094000 : test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir)
726 2094000 : work5(ir)=test
727 2094000 : work1(ir)=0.0d0
728 : ! Will ignore tail within decade of machine precision
729 2095047 : if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then
730 1144491 : irmax=ir
731 : end if
732 : end do
733 : ! Increase irmax a bit
734 1047 : irmax=irmax+4
735 : ! Ask irmax to be lower than mmax
736 1047 : if(irmax>mmax-1)irmax=mmax-1
737 :
738 3141 : ABI_MALLOC(work3,(irmax-1))
739 2094 : ABI_MALLOC(work4,(irmax-1))
740 :
741 : ! Loop over q values
742 3143177 : do iq=2,mqgrid
743 3142130 : tpiq=two_pi*qgrid(iq)
744 3142130 : call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
745 2868399040 : work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility)
746 2868399040 : work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility)
747 :
748 : ! Handle r=0 separately
749 3142130 : work1(1)=0.0d0
750 3142130 : call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax))
751 2868399040 : do ir=2,irmax
752 2868399040 : work1(ir)=besjx(ir)*work5(ir)
753 : end do
754 : ! do integral
755 3142130 : call der_int(work1,work2,rad,dr,irmax,result)
756 3143177 : ffspl(iq,1,lp1)=result
757 : end do
758 :
759 : ! Compute yp1=derivative of f(q) at q=0
760 1047 : if(lp1/=2)then
761 : ! For l/=1, yp1=0
762 636 : yp1=0.0d0
763 : else
764 : ! For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3
765 379700 : do ir=1,irmax
766 379289 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
767 379700 : work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0
768 : end do
769 411 : call der_int(work1,work2,rad,dr,irmax,result)
770 411 : yp1=result
771 : end if
772 :
773 : ! Compute ypn=derivative of f(q) at q=qgrid(mqgrid)
774 1047 : tpiq=two_pi*qgrid(mqgrid)
775 : ! Treat ir=1, r=0, separately
776 1047 : work1(1)=0.0d0
777 : ! Here, must distinguish l==0 from others
778 1047 : if(lp1==1)then
779 : ! l==0 : ypn=$\int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$
780 : ! The sine and cosine of the last point were computed in the previous loop
781 : ! So, there is no need to call sincos. Note that the rank of besj is 1.
782 553 : call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax))
783 495803 : do ir=2,irmax
784 495250 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
785 495803 : work1(ir)=-besjx(ir)*two_pi*rad(ir)*rad(ir)*dvwf
786 : end do
787 : else
788 : ! l==1 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$
789 : ! l==2 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$
790 : ! l==3 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$
791 : ! The sine and cosine of the last point were computed in the previous loop
792 : ! Store first previously computed value with besj of order l, then use
793 : ! besj of order l-1 (=lp1-2)
794 459988 : work1(2:irmax)=besjx(2:irmax)
795 494 : call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax))
796 459988 : do ir=2,irmax
797 459494 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
798 : work1(ir)=(two_pi*rad(ir)**2)*dvwf*&
799 459988 : & ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) )
800 : end do
801 : end if
802 : ! work1 is ready for integration
803 1047 : call der_int(work1,work2,rad,dr,irmax,result)
804 1047 : ypn=result
805 :
806 : ! Fit spline to get second derivatives by spline fit
807 : call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,&
808 1047 : & ffspl(:,2,lp1))
809 :
810 1047 : ABI_FREE(work3)
811 2094 : ABI_FREE(work4)
812 :
813 : else
814 : ! KB energy is zero, put nonlocal correction at l=0 to 0
815 498415 : ffspl(:,:,lp1)=0.0d0
816 : end if
817 :
818 : end do ! End loop on angular momenta
819 :
820 636 : ABI_FREE(work1)
821 636 : ABI_FREE(work2)
822 636 : ABI_FREE(work_spl)
823 636 : ABI_FREE(work5)
824 636 : ABI_FREE(besjx)
825 : end if ! End of lmax/=-1 condition
826 :
827 638 : end subroutine psp1nl
828 : !!***
829 :
830 : !!****f* m_psp1/der_int
831 : !! NAME
832 : !! der_int
833 : !!
834 : !! FUNCTION
835 : !! Given input function f(i) on Teter radial grid, and grid spacing
836 : !! dr(i), compute function derivative df/dr on points from 0 to n.
837 : !! Integrate function f(i) on grid r(i) from r(0) to r(nlast).
838 : !! Note that array dimensions start at 0.
839 : !!
840 : !! INPUTS
841 : !! f(0 to nlast)=function values on grid
842 : !! r(0 to nlast)=radial grid points
843 : !! dr(0 to nlast)=INVERSE of spacing on grid
844 : !! nlast=radial grid point for upper limit
845 : !!
846 : !! OUTPUT
847 : !! df(0 to n)=derivative $ \frac{df}{dr}$ on grid
848 : !! smf= $ \int_{r(0)}^{r(nlast)} f(r) dr $.
849 : !!
850 : !! SOURCE
851 :
852 5230188 : subroutine der_int(ff,df,rr,dr,nlast,smf)
853 :
854 : !Arguments ------------------------------------
855 : !nmax sets standard number of grid points ! SHOULD BE REMOVED
856 : !scalars
857 : integer,parameter :: nmax=2000
858 : integer,intent(in) :: nlast
859 : real(dp),intent(out) :: smf
860 : !no_abirules
861 : !Note that dimension here starts at 0
862 : real(dp), intent(in) :: dr(0:nmax),ff(0:nmax),rr(0:nmax)
863 : real(dp), intent(out) :: df(0:nmax)
864 :
865 : !Local variables-------------------------------
866 : !scalars
867 : integer :: jj
868 : real(dp),parameter :: div12=1.d0/12.d0
869 : real(dp) :: hh
870 : character(len=500) :: message
871 :
872 : ! *************************************************************************
873 :
874 : !Check that nlast lie within 0 to nmax
875 5230188 : if (nlast<0.or.nlast>nmax) then
876 : write(message, '(a,i12,a,i12)' )&
877 0 : & ' nlast=',nlast,' lies outside range [0,nmax] with dimension nmax=',nmax
878 0 : ABI_BUG(message)
879 : end if
880 :
881 : !Compute derivatives at lower end, near r=0
882 : df(0)=-25.d0/12.d0*ff(0)+4.d0*ff(1)-3.d0*ff(2)+4.d0/3.d0*ff(3)&
883 5230188 : & -1.d0/4.d0*ff(4)
884 : df(1)=-1.d0/4.d0*ff(0)-5.d0/6.d0*ff(1)+3.d0/2.d0*ff(2)&
885 5230188 : & -1.d0/2.d0*ff(3)+1.d0/12.d0*ff(4)
886 :
887 : !Run over range from just past r=0 to near r(n), using central differences
888 5308385239 : do jj=2,nlast-2
889 5308385239 : df(jj)=(ff(jj-2)-8.d0*(ff(jj-1)-ff(jj+1))-ff(jj+2))*div12
890 : end do
891 :
892 : !Compute derivative at upper end of range
893 5230188 : if (nlast < 4) then
894 0 : message = ' der_int: ff does not have enough elements. nlast is too low'
895 0 : ABI_ERROR(message)
896 : end if
897 :
898 : df(nlast-1)=-1.d0/12.d0*ff(nlast-4)&
899 : & +1.d0/2.d0*ff(nlast-3)&
900 : & -3.d0/2.d0*ff(nlast-2)&
901 : & +5.d0/6.d0*ff(nlast-1)&
902 5230188 : & +1.d0/4.d0*ff(nlast)
903 : df(nlast)=1.d0/4.d0*ff(nlast-4)&
904 : & -4.d0/3.d0*ff(nlast-3)&
905 : & +3.d0*ff(nlast-2)&
906 : & -4.d0*ff(nlast-1)&
907 5230188 : & +25.d0/12.d0*ff(nlast)
908 :
909 : !Apply correct normalization over full range
910 5329305991 : do jj=0,nlast
911 5329305991 : df(jj)=df(jj)*dr(jj)
912 : end do
913 :
914 5230188 : smf=0.d0
915 5324075803 : do jj=0,nlast-1
916 5318845615 : hh=rr(jj+1)-rr(jj)
917 5324075803 : smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1)))
918 : end do
919 5230188 : smf=smf/12.d0
920 :
921 5230188 : end subroutine der_int
922 : !!***
923 :
924 : !!****f* m_psp1/sincos
925 : !! NAME
926 : !! sincos
927 : !!
928 : !! FUNCTION
929 : !! Update the sine and cosine values, needed inside the
930 : !! pseudopotential routines psp1lo and psp1nl.
931 : !!
932 : !! INPUTS
933 : !! iq = number of current wavevector q
934 : !! irmax = number of values of r on the radial grid to be computed
935 : !! mmax = dimension of pspwk and rad
936 : !! pspwk(:,1,1) and pspwk(:,2,1) : sine and cosine of 2$\pi$ dq * rad
937 : !! pspwk(:,1,2) and pspwk(:,2,2) : sine and cosine of 2$\pi$ previous q * rad
938 : !! rad(mmax) radial grid
939 : !! tpiq = 2 $\pi$ * current wavevector q
940 : !!
941 : !! OUTPUT
942 : !! pspwk(*,1,2) and pspwk(*,2,2) : sine and cosine of 2$\pi$ current q * rad
943 : !!
944 : !! NOTES
945 : !! The speed was a special concern, so iterative computation
946 : !! based on addition formula is possible. Interestingly,
947 : !! this algorithm places strong constraints on accuracy,
948 : !! so this routine is machine-dependent.
949 : !!
950 : !! SOURCE
951 :
952 5224695 : subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq)
953 :
954 : !Arguments ------------------------------------
955 : !scalars
956 : integer,intent(in) :: iq,irmax,mmax
957 : real(dp),intent(in) :: tpiq
958 : !arrays
959 : real(dp),intent(in) :: rad(mmax)
960 : real(dp),intent(inout) :: pspwk(mmax,2,2)
961 :
962 : !Local variables-------------------------------
963 : !scalars
964 : integer :: ir,nstep
965 : real(dp) :: prevcos,prevsin
966 : logical :: testmipspro
967 :
968 :
969 : ! *************************************************************************
970 :
971 5224695 : if(iq==2)then
972 :
973 : ! Here set up the sin and cos at iq=2
974 1769298 : do ir=2,irmax
975 1767557 : pspwk(ir,1,2)=pspwk(ir,1,1)
976 1769298 : pspwk(ir,2,2)=pspwk(ir,2,1)
977 : end do
978 :
979 : else
980 : !
981 : ! The sensitivity of the algorithm to changes of nstep
982 : ! has been tested : for all the machines except SGI - R10000 ,
983 : ! either using only the hard way, or
984 : ! using up to nstep=40 causes changes at the level
985 : ! of 1.0d-16 in the total energy. Larger values of
986 : ! nstep might be possible, but the associated residual
987 : ! is already very small ! The accelerated computation of
988 : ! sine and cosine is essential for a good speed on IBM, but,
989 : ! fortunately, on the SGI - R10000 the normal computation is fast enough.
990 :
991 5222954 : testmipspro=.false.
992 5222954 : nstep=40
993 5222954 : if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then
994 :
995 : ! Every nstep steps, uses the hard way
996 132735646 : do ir=2,irmax
997 132605029 : pspwk(ir,1,2)=sin(tpiq*rad(ir))
998 132735646 : pspwk(ir,2,2)=cos(tpiq*rad(ir))
999 : end do
1000 :
1001 : else
1002 :
1003 : ! Here the fastest way, iteratively
1004 5174935629 : do ir=2,irmax
1005 5169843292 : prevsin=pspwk(ir,1,2)
1006 5169843292 : prevcos=pspwk(ir,2,2)
1007 5169843292 : pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1)
1008 5174935629 : pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1)
1009 : end do
1010 :
1011 : end if
1012 :
1013 : end if ! iq==2
1014 :
1015 5224695 : end subroutine sincos
1016 : !!***
1017 :
1018 : !!****f* m_psp1/psp4cc
1019 : !! NAME
1020 : !! psp4cc
1021 : !!
1022 : !! FUNCTION
1023 : !! Compute the core charge density, for use in the XC core
1024 : !! correction, following the function definition valid
1025 : !! for the format 4 of pseudopotentials.
1026 : !! This is a even polynomial of 24th order for core density,
1027 : !! that is cut off exactly beyond rchrg.
1028 : !! It has been produced on 7 May 1992 by M. Teter.
1029 : !!
1030 : !! INPUTS
1031 : !! fchrg=magnitude of the core charge correction
1032 : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
1033 : !!
1034 : !! OUTPUT
1035 : !! xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
1036 : !!
1037 : !! NOTES
1038 : !! The argument of xccc1d is assumed to be normalized, and to vary
1039 : !! from xx=0 to 1 (from r=0 to r=xcccrc)
1040 : !!
1041 : !! WARNINGS
1042 : !! the fifth derivative is not yet delivered.
1043 : !!
1044 : !! SOURCE
1045 :
1046 3 : subroutine psp4cc(fchrg,n1xccc,xccc1d)
1047 :
1048 : !Arguments ------------------------------------
1049 : !scalars
1050 : integer,intent(in) :: n1xccc
1051 : real(dp),intent(in) :: fchrg
1052 : !arrays
1053 : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
1054 :
1055 : !Local variables-------------------------------
1056 : !scalars
1057 : integer :: i1xccc,ider
1058 : real(dp),parameter :: a10=-0.1156854803757563d5,a12=+0.2371534625455588d5
1059 : real(dp),parameter :: a14=-0.3138755797827918d5,a16=+0.2582842713241039d5
1060 : real(dp),parameter :: a18=-0.1200356429115204d5,a20=+0.2405099057118771d4
1061 : real(dp),parameter :: a2=-0.8480751097855989d1,a4=+0.9684600878284791d2
1062 : real(dp),parameter :: a6=-0.7490894651588015d3,a8=+0.3670890998130434d4
1063 : real(dp) :: der1,dern,factor
1064 : character(len=500) :: message
1065 : !arrays
1066 3 : real(dp),allocatable :: ff(:),ff2(:),work(:),xx(:)
1067 : real(dp) :: x
1068 :
1069 : ! *************************************************************************
1070 :
1071 9 : ABI_MALLOC(ff,(n1xccc))
1072 6 : ABI_MALLOC(ff2,(n1xccc))
1073 6 : ABI_MALLOC(work,(n1xccc))
1074 6 : ABI_MALLOC(xx,(n1xccc))
1075 :
1076 :
1077 3 : if(n1xccc > 1)then
1078 3 : factor=1.0d0/dble(n1xccc-1)
1079 7506 : do i1xccc=1,n1xccc
1080 7506 : xx(i1xccc)=(i1xccc-1)*factor
1081 : end do
1082 : else
1083 0 : write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc
1084 0 : ABI_BUG(message)
1085 : end if
1086 :
1087 : !Initialization, to avoid some problem with some compilers
1088 39 : xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
1089 :
1090 : !Take care of each derivative separately
1091 12 : do ider=0,2
1092 :
1093 9 : if(ider==0)then
1094 : ! Generate spline fitting for the function gg
1095 7506 : do i1xccc=1,n1xccc
1096 : ! ff(i1xccc)=fchrg*gg(xx(i1xccc))
1097 7506 : ff(i1xccc)=fchrg*gg_psp4(xx(i1xccc))
1098 : end do
1099 : ! Complete with derivatives at end points
1100 3 : der1=0.0d0
1101 : ! dern=fchrg*gp(1.0d0)
1102 3 : dern=fchrg*gp_psp4(1.0d0)
1103 6 : else if(ider==1)then
1104 : ! Generate spline fitting for the function gp
1105 7506 : do i1xccc=1,n1xccc
1106 : ! ff(i1xccc)=fchrg*gp(xx(i1xccc))
1107 7506 : ff(i1xccc)=fchrg*gp_psp4(xx(i1xccc))
1108 : end do
1109 : ! Complete with derivatives at end points, already estimated
1110 3 : der1=xccc1d(1,ider+2)
1111 3 : dern=xccc1d(n1xccc,ider+2)
1112 : else if(ider==2)then
1113 : ! Generate spline fitting for the function gpp
1114 : ! (note : the function gpp has already been estimated, for the spline
1115 : ! fitting of the function gg, but it is replaced here by the more
1116 : ! accurate analytic derivative)
1117 7506 : do i1xccc=1,n1xccc
1118 7503 : x=xx(i1xccc)
1119 7506 : ff(i1xccc)=fchrg*(gpp_1_psp4(x)+gpp_2_psp4(x)+gpp_3_psp4(x))
1120 : ! ff(i1xccc)=fchrg*gpp(xx(i1xccc))
1121 : end do
1122 : ! Complete with derivatives of end points
1123 3 : der1=xccc1d(1,ider+2)
1124 3 : dern=xccc1d(n1xccc,ider+2)
1125 : end if
1126 :
1127 : ! Produce second derivative numerically, for use with splines
1128 9 : call spline(xx,ff,n1xccc,der1,dern,ff2)
1129 22518 : xccc1d(:,ider+1)=ff(:)
1130 22521 : xccc1d(:,ider+3)=ff2(:)
1131 : end do
1132 :
1133 7506 : xccc1d(:,6)=zero
1134 :
1135 3 : ABI_FREE(ff)
1136 3 : ABI_FREE(ff2)
1137 3 : ABI_FREE(work)
1138 3 : ABI_FREE(xx)
1139 :
1140 : !DEBUG
1141 : !write(std_out,*)' psp1cc : output of core charge density and derivatives '
1142 : !write(std_out,*)' xx gg gp '
1143 : !do i1xccc=1,n1xccc
1144 : !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
1145 : !end do
1146 : !write(std_out,*)' xx gpp gg2 '
1147 : !do i1xccc=1,n1xccc
1148 : !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
1149 : !end do
1150 : !write(std_out,*)' xx gp2 gpp2 '
1151 : !do i1xccc=1,n1xccc
1152 : !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
1153 : !end do
1154 : !write(std_out,*)' psp1cc : debug done, stop '
1155 : !stop
1156 : !ENDDEBUG
1157 :
1158 : contains
1159 :
1160 7503 : function gg_psp4(x)
1161 : !Expression of 7 May 1992
1162 : real(dp) :: gg_psp4
1163 : real(dp),intent(in) :: x
1164 : gg_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + &
1165 : & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ &
1166 7503 : & x**2*(a18+x**2*(a20))))))))))) *(1.0d0-x**2)**2
1167 7503 : end function gg_psp4
1168 :
1169 7506 : function gp_psp4(x)
1170 : !gp(x) is the derivative of gg(x) wrt x
1171 : real(dp) :: gp_psp4
1172 : real(dp),intent(in) :: x
1173 : gp_psp4=2.d0*x*((a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( &
1174 : & 4.d0*a8+x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( &
1175 : & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(10.d0*a20))))))))))*&
1176 : & (1.d0-x**2)**2 &
1177 : & -2.0d0*(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + &
1178 : & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ &
1179 7506 : & x**2*(a18+x**2*a20)))))))))) *(1.0d0-x**2) )
1180 7506 : end function gp_psp4
1181 :
1182 7503 : function gpp_1_psp4(x)
1183 : !gpp(x) is the second derivative of gg(x) wrt x
1184 : real(dp) :: gpp_1_psp4
1185 : real(dp),intent(in) :: x
1186 : gpp_1_psp4= ( 2.d0*a4+ x**2*(3.d0*2.d0*a6 +x**2*( &
1187 : & 4.d0*3.d0*a8+ x**2*(5.d0*4.d0*a10+x**2*( &
1188 : & 6.d0*5.d0*a12+x**2*(7.d0*6.d0*a14+x**2*( &
1189 : & 8.d0*7.d0*a16+x**2*(9.d0*8.d0*a18+x**2*( &
1190 : & 10.d0*9.d0*a20) &
1191 7503 : & ))))))))*(2.d0*x*(1.d0-x**2))**2
1192 7503 : end function gpp_1_psp4
1193 :
1194 7503 : function gpp_2_psp4(x)
1195 :
1196 : real(dp) :: gpp_2_psp4
1197 : real(dp),intent(in) :: x
1198 : gpp_2_psp4=(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( &
1199 : & 4.d0*a8 +x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( &
1200 : & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*( &
1201 : & 10.d0*a20) &
1202 7503 : & )))))))))*(1.d0-x**2)*2*(1.d0-9.d0*x**2)
1203 7503 : end function gpp_2_psp4
1204 :
1205 7503 : function gpp_3_psp4(x)
1206 :
1207 : real(dp) :: gpp_3_psp4
1208 : real(dp),intent(in) :: x
1209 : gpp_3_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + &
1210 : & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ &
1211 : & x**2*(a18+x**2*a20 &
1212 7503 : & ))))))))))*(1.0d0-3.d0*x**2)*(-4.d0)
1213 7503 : end function gpp_3_psp4
1214 :
1215 : end subroutine psp4cc
1216 : !!***
1217 :
1218 : end module m_psp1
1219 : !!***
|
