Branch data 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 : 696 : subroutine psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,&
107 : 696 : & e990,e999,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
108 : 696 : & mmax,mpsang,mqgrid,nproj,n1xccc,pspcod,&
109 : 696 : & 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 : 696 : real(dp),allocatable :: drad(:),ekb_tmp(:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:)
134 : 696 : real(dp),allocatable :: vpspll(:,:),wfll(:,:),wksincos(:,:,:),work_space(:)
135 : 696 : 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 [ - + ]: 696 : 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 [ + + ]: 2360 : do ipsang=1,lmax+1
166 : :
167 [ + - ]: 1664 : 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 : 1664 : & ipsang-1,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),&
170 : 3328 : & 'l,e99.0,e99.9,nproj,rcpsp'
171 : 1664 : call wrtout(ab_out,message,'COLL')
172 : 1664 : call wrtout(std_out, message,'COLL')
173 : :
174 [ + - ]: 1664 : read (tmp_unit,*,err=10,iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang)
175 : : write(message, '(4f13.8,t55,a)' ) &
176 : 1664 : & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),&
177 : 3328 : & ' rms, ekb1, ekb2, epsatm'
178 : 1664 : call wrtout(ab_out,message,'COLL')
179 : 2360 : 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 [ + + + + ]: 9943 : index=0;iln=0;indlmn(:,:)=0
185 [ + + ]: 2360 : do ipsang=1,lmax+1
186 [ + + ]: 2360 : if(nproj(ipsang)>0)then
187 : 970 : ll=ipsang-1
188 [ + + ]: 2020 : do kk=1,nproj(ipsang)
189 : 1050 : iln=iln+1
190 [ + + ]: 3086 : do mm=1,2*ll*useylm+1
191 : 1066 : index=index+1
192 : 1066 : indlmn(1,index)=ll
193 : 1066 : indlmn(2,index)=mm-ll*useylm-1
194 : 1066 : indlmn(3,index)=kk
195 : 1066 : indlmn(4,index)=ll*ll+(1-useylm)*ll+mm
196 : 1066 : indlmn(5,index)=iln
197 : 2116 : indlmn(6,index)=1
198 : : end do
199 : : end do
200 : : end if
201 : : end do
202 : :
203 [ + - ]: 696 : read (tmp_unit,*,err=10,iomsg=errmsg) rchrg,fchrg,qchrg
204 : 696 : write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
205 : 696 : call wrtout(ab_out,message,'COLL')
206 : 696 : call wrtout(std_out, message,'COLL')
207 : :
208 : : ! Generate core charge function and derivatives, if needed
209 [ + + ]: 696 : if(fchrg>1.0d-15)then
210 [ + + ]: 591 : if(pspcod==1)then
211 : 588 : 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 : 588 : 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 [ - + + - : 2088 : ABI_MALLOC(vloc,(mmax))
+ - + - -
+ - + ]
235 [ - + + - : 2784 : ABI_MALLOC(vpspll,(mmax,mpsang))
- + + - +
- + - - +
- + ]
236 [ + + + + : 4700 : 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 [ + + ]: 2360 : do ipsang=1,lmax+1
241 [ + - ]: 1664 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
242 [ + - ]: 2360 : 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 [ + + ]: 1393392 : vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 )
247 : :
248 : : !(2) Create radial grid, and associated quantities
249 [ - + + - : 1392 : ABI_MALLOC(rad,(mmax))
- + - + ]
250 [ - + + - : 1392 : ABI_MALLOC(drad,(mmax))
- + - + ]
251 [ - + + - : 2784 : ABI_MALLOC(wksincos,(mmax,2,2))
+ - + - +
- - + -
+ ]
252 : :
253 : : !Teter grid--need both r and dr in this case
254 [ + + ]: 1393392 : do ii=0,mmax-1
255 : 1392696 : xx=dble(ii)/dble(mmax-1)
256 : 1392696 : rad (ii+1)=100.d0*(xx+.01d0)**5-1.d-8
257 : 1393392 : 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 : 696 : dq2pi=2.0d0*pi*dq
263 [ + + ]: 1393392 : do ii=1,mmax
264 : 1392696 : arg=dq2pi*rad(ii)
265 : 1392696 : drad(ii)=1.0d0/drad(ii)
266 : 1392696 : wksincos(ii,1,1)=sin(arg)
267 : 1393392 : 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 [ - + + - : 2088 : ABI_MALLOC(work_space,(mqgrid))
+ - + - -
+ - + ]
273 [ - + + - : 1392 : ABI_MALLOC(work_spl1,(mqgrid))
- + - + ]
274 [ - + + - : 1392 : ABI_MALLOC(work_spl2,(mqgrid))
- + - + ]
275 : : call psp1lo(drad,epsatm,mmax,mqgrid,qgrid,&
276 : 696 : & work_spl1,rad,vloc,wksincos,yp1,ypn,zion)
277 : :
278 : : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
279 : 696 : call spline (qgrid,work_spl1,mqgrid,yp1,ypn,work_spl2)
280 [ + + ]: 2089957 : vlspl(:,1)=work_spl1(:)
281 [ + + ]: 2089957 : vlspl(:,2)=work_spl2(:)
282 : :
283 [ - + ]: 696 : ABI_FREE(work_space)
284 [ - + ]: 696 : ABI_FREE(work_spl1)
285 [ - + ]: 696 : ABI_FREE(work_spl2)
286 : :
287 : : !(4)Take care of non-local part
288 : :
289 : : !Zero out all Kleinman-Bylander energies to initialize
290 [ + + ]: 1987 : 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 [ + + + + ]: 696 : 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 [ - + + - : 1797 : ABI_MALLOC(wfll,(mmax,mpsang))
- + + - -
+ - + ]
306 [ + + ]: 2166 : do ipsang=1,lmax+1
307 : : ! For pspcod==4, wfs for the local angular momentum are not written
308 [ + + + + ]: 2166 : if (nproj(ipsang)/=0 .or. pspcod==1) then
309 [ + - ]: 1558 : read (tmp_unit,*,err=10,iomsg=errmsg) ll
310 [ - + ]: 1558 : 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 [ + - ]: 1558 : 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 : 599 : nlmax=lmax
330 [ + + ]: 599 : if (lloc==lmax) nlmax=lmax-1
331 : : ! write(std_out,*)' psp1in : lmax,lloc=',lmax,lloc
332 [ - + + - : 1797 : ABI_MALLOC(ekb_tmp,(mpsang,2))
+ - + - -
+ - + ]
333 [ - + + - : 2995 : ABI_MALLOC(ffspl_tmp,(mqgrid,2,nlmax+1,2))
+ + + + +
- + - + -
+ + - + -
+ ]
334 : :
335 : : call psp1nl(drad,ekb_tmp(:,1),ffspl_tmp(:,:,:,1),lloc,&
336 : 599 : & 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 : 599 : lhigh=-1
341 [ + + ]: 2166 : do ipsang=1,min(lmax+1,mpsang)
342 [ + + ]: 2166 : 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 [ + + ]: 2976974 : 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 [ + + ]: 599 : 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 [ + + ]: 1817 : do ii=1,lmnmax
371 : 1218 : kk=indlmn(5,ii)
372 [ + + ]: 1817 : if (kk>iln) then
373 : 1050 : iln=kk
374 : 1050 : 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 [ + + + + ]: 6307510 : ffspl(:,:,kk)=ffspl_tmp(:,:,1+indlmn(1,ii),indlmn(3,ii))
377 : : end if
378 : : end do
379 : :
380 [ - + ]: 599 : ABI_FREE(ekb_tmp)
381 [ - + ]: 599 : ABI_FREE(ffspl_tmp)
382 [ - + ]: 599 : ABI_FREE(wfll)
383 : : end if
384 : :
385 [ - + ]: 696 : ABI_FREE(vpspll)
386 [ - + ]: 696 : ABI_FREE(rad)
387 [ - + ]: 696 : ABI_FREE(drad)
388 [ - + ]: 696 : ABI_FREE(vloc)
389 [ - + ]: 696 : ABI_FREE(wksincos)
390 : :
391 : 696 : return
392 : :
393 : : ! Handle IO error
394 : : 10 continue
395 : 0 : ABI_ERROR(errmsg)
396 : :
397 [ - + - + : 1392 : 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 : 696 : subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
431 : 696 : & 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 : 696 : wk(1)=0.0d0
458 : :
459 [ + + ]: 1392696 : 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 : 1392000 : 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 [ + + ]: 1392696 : if (abs(test)<1.0d-20) then
468 : 404102 : wk(ir)=0.0d0
469 : : else
470 : 987898 : 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 : 696 : call der_int(wk,wk2,rad,drad,mmax-1,result)
477 : :
478 : 696 : epsatm=4.d0*pi*(result)
479 : : !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi
480 : 696 : q2vq(1)=-zion/pi
481 : :
482 : : !Prepare loop over q values
483 : 696 : irmax=mmax+1
484 [ + + ]: 1392696 : do ir=mmax,2,-1
485 : 1392000 : test=vloc(ir)+zion/rad(ir)
486 : 1392000 : wk1(ir)=test*rad(ir)
487 : : ! Will ignore tail within decade of machine precision
488 [ + + + + ]: 1392696 : if ((scale+abs(test))==scale .and. irmax==ir+1) then
489 : 580401 : irmax=ir
490 : : end if
491 : : end do
492 : : !Increase irmax a bit : this is copied from psp1nl
493 : 696 : irmax=irmax+4
494 [ + + ]: 696 : if(irmax>mmax)irmax=mmax
495 : :
496 : : !Loop over q values
497 [ + + ]: 2089261 : do iq=2,mqgrid
498 : 2088565 : tpiq=two_pi*qgrid(iq)
499 : 2088565 : 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 : 2088565 : wk(1)=0.0d0
504 [ + + ]: 2447659533 : do ir=2,irmax
505 : 2447659533 : wk(ir)=wksincos(ir,1,2)*wk1(ir)
506 : : end do
507 : : ! do integral from 0 to r(max)
508 [ + + ]: 2088565 : if(irmax>mmax-1)irmax=mmax-1
509 : :
510 : 2088565 : call der_int(wk,wk2,rad,drad,irmax,result)
511 : : ! store q^2 v(q)
512 : 2089261 : 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 : 696 : 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 : 696 : wk(1)=0.0d0
522 : 696 : tpiq=two_pi*qgrid(mqgrid)
523 [ + + ]: 1392696 : do ir=2,mmax
524 : 1392000 : test=vloc(ir)+zion/rad(ir)
525 : : ! Ignore contributions within decade of machine precision
526 [ + + ]: 1392696 : if ((scale+abs(test))==scale) then
527 : 581786 : wk(ir)=0.0d0
528 : : else
529 : : wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * &
530 : 810214 : & (rad(ir)*vloc(ir)+zion)
531 : : end if
532 : : end do
533 : 696 : call der_int(wk,wk2,rad,drad,mmax-1,result)
534 : :
535 : 696 : ypn=2.0d0*result
536 : :
537 : 696 : 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 : 640 : subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,&
590 : 640 : & 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 : 640 : real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:)
610 : 640 : real(dp),allocatable :: work_spl(:)
611 : :
612 : : ! *************************************************************************
613 : :
614 : : !Zero out Kleinman-Bylander energies ekb
615 [ + + ]: 2400 : 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 [ + + ]: 640 : if (lmax/=-1) then
623 : :
624 : : ! Check that lmax is within allowed range
625 [ - + ]: 638 : 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 [ + - - + : 1914 : ABI_MALLOC(work1,(mmax+1))
+ - + - +
- - + -
+ ]
639 [ - + + - : 1276 : ABI_MALLOC(work2,(mmax+1))
- + - + ]
640 [ - + + - : 1914 : ABI_MALLOC(work_spl,(mqgrid))
+ - + - -
+ - + ]
641 [ - + + - : 1914 : ABI_MALLOC(work5,(mmax))
+ - + - -
+ - + ]
642 [ - + + - : 1276 : ABI_MALLOC(besjx,(mmax))
- + - + ]
643 : :
644 [ + + ]: 1771 : do lp1=1,lmax+1
645 : :
646 : : ! Only do the work if nonlocal correction is nonzero
647 [ + + ]: 1771 : if (lp1 /= lloc+1) then
648 : :
649 : : ! integrand for 0 to r(mmax)
650 [ + + ]: 2102100 : do ir=1,mmax
651 : 2101050 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
652 : 2102100 : 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 : 1050 : call der_int(work1,work2,rad,dr,mmax-1,result)
659 : 1050 : 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 [ + + ]: 2102100 : do ir=1,mmax
674 : 2101050 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
675 : 2102100 : work1(ir)=dvwf**2
676 : : end do
677 : 1050 : call der_int(work1,work2,rad,dr,mmax-1,result)
678 : :
679 : 1050 : 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 [ + - ]: 1050 : if (dvms(lp1)/=0.0d0) then
693 : 1050 : ekb(lp1)=dvms(lp1)/eta(lp1)
694 : 1050 : 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 [ + + ]: 1771 : do lp1=1,lmax+1
705 : :
706 : : ! Compute form factor if ekb(lp1) not 0
707 [ + + ]: 1771 : if (ekb(lp1)/=0.0d0) then
708 : :
709 : : ! do q=0 separately, non-zero if l=0
710 [ + + ]: 1050 : if(lp1==1)then
711 [ + + ]: 1111110 : do ir=1,mmax
712 : 1110555 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
713 : 1111110 : work1(ir)=rad(ir)*dvwf
714 : : end do
715 : 555 : call der_int(work1,work2,rad,dr,mmax-1,result)
716 : 555 : ffspl(1,1,lp1)=result
717 : : else
718 : : ! For l non-zero, f(q=0) vanishes !
719 : 495 : ffspl(1,1,lp1)=0.0d0
720 : : end if
721 : :
722 : : ! Prepare loop over q values
723 : 1050 : irmax=mmax+1
724 [ + + ]: 2101050 : do ir=mmax,2,-1
725 : 2100000 : test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir)
726 : 2100000 : work5(ir)=test
727 : 2100000 : work1(ir)=0.0d0
728 : : ! Will ignore tail within decade of machine precision
729 [ + + + + ]: 2101050 : if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then
730 : 1147846 : irmax=ir
731 : : end if
732 : : end do
733 : : ! Increase irmax a bit
734 : 1050 : irmax=irmax+4
735 : : ! Ask irmax to be lower than mmax
736 [ - + ]: 1050 : if(irmax>mmax-1)irmax=mmax-1
737 : :
738 [ + - - + : 3150 : ABI_MALLOC(work3,(irmax-1))
+ - + - +
- - + -
+ ]
739 [ - + + - : 2100 : ABI_MALLOC(work4,(irmax-1))
- + - + ]
740 : :
741 : : ! Loop over q values
742 [ + + ]: 3152180 : do iq=2,mqgrid
743 : 3151130 : tpiq=two_pi*qgrid(iq)
744 : 3151130 : call sincos(iq,irmax,mmax,wksincos,rad,tpiq)
745 [ + + ]: 2876388040 : work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility)
746 [ + + ]: 2876388040 : work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility)
747 : :
748 : : ! Handle r=0 separately
749 : 3151130 : work1(1)=0.0d0
750 : 3151130 : call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax))
751 [ + + ]: 2876388040 : do ir=2,irmax
752 : 2876388040 : work1(ir)=besjx(ir)*work5(ir)
753 : : end do
754 : : ! do integral
755 : 3151130 : call der_int(work1,work2,rad,dr,irmax,result)
756 : 3152180 : ffspl(iq,1,lp1)=result
757 : : end do
758 : :
759 : : ! Compute yp1=derivative of f(q) at q=0
760 [ + + ]: 1050 : if(lp1/=2)then
761 : : ! For l/=1, yp1=0
762 : 638 : yp1=0.0d0
763 : : else
764 : : ! For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3
765 [ + + ]: 380609 : do ir=1,irmax
766 : 380197 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
767 : 380609 : work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0
768 : : end do
769 : 412 : call der_int(work1,work2,rad,dr,irmax,result)
770 : 412 : yp1=result
771 : : end if
772 : :
773 : : ! Compute ypn=derivative of f(q) at q=qgrid(mqgrid)
774 : 1050 : tpiq=two_pi*qgrid(mqgrid)
775 : : ! Treat ir=1, r=0, separately
776 : 1050 : work1(1)=0.0d0
777 : : ! Here, must distinguish l==0 from others
778 [ + + ]: 1050 : 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 : 555 : call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax))
783 [ + + ]: 497558 : do ir=2,irmax
784 : 497003 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
785 : 497558 : 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 [ + + ]: 460896 : work1(2:irmax)=besjx(2:irmax)
795 : 495 : call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax))
796 [ + + ]: 460896 : do ir=2,irmax
797 : 460401 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
798 : : work1(ir)=(two_pi*rad(ir)**2)*dvwf*&
799 : 460896 : & ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) )
800 : : end do
801 : : end if
802 : : ! work1 is ready for integration
803 : 1050 : call der_int(work1,work2,rad,dr,irmax,result)
804 : 1050 : ypn=result
805 : :
806 : : ! Fit spline to get second derivatives by spline fit
807 : : call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,&
808 : 1050 : & ffspl(:,2,lp1))
809 : :
810 [ - + ]: 1050 : ABI_FREE(work3)
811 [ - + ]: 2100 : 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 [ - + ]: 638 : ABI_FREE(work1)
821 [ - + ]: 638 : ABI_FREE(work2)
822 [ - + ]: 638 : ABI_FREE(work_spl)
823 [ - + ]: 638 : ABI_FREE(work5)
824 [ - + ]: 638 : ABI_FREE(besjx)
825 : : end if ! End of lmax/=-1 condition
826 : :
827 : 640 : 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 : 5245204 : 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 [ - + ]: 5245204 : 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 : 5245204 : & -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 : 5245204 : & -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 [ + + ]: 5322989778 : do jj=2,nlast-2
889 : 5322989778 : 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 [ - + ]: 5245204 : 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 : 5245204 : & +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 : 5245204 : & +25.d0/12.d0*ff(nlast)
908 : :
909 : : !Apply correct normalization over full range
910 [ + + ]: 5343970594 : do jj=0,nlast
911 : 5343970594 : df(jj)=df(jj)*dr(jj)
912 : : end do
913 : :
914 : 5245204 : smf=0.d0
915 [ + + ]: 5338725390 : do jj=0,nlast-1
916 : 5333480186 : hh=rr(jj+1)-rr(jj)
917 : 5338725390 : smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1)))
918 : : end do
919 : 5245204 : smf=smf/12.d0
920 : :
921 : 5245204 : 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 : 5239695 : 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 [ + + ]: 5239695 : if(iq==2)then
972 : :
973 : : ! Here set up the sin and cos at iq=2
974 [ + + ]: 1774167 : do ir=2,irmax
975 : 1772421 : pspwk(ir,1,2)=pspwk(ir,1,1)
976 : 1774167 : 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 : 5237949 : testmipspro=.false.
992 : 5237949 : nstep=40
993 [ + + ]: 5237949 : if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then
994 : :
995 : : ! Every nstep steps, uses the hard way
996 [ + + ]: 133100821 : do ir=2,irmax
997 : 132969829 : pspwk(ir,1,2)=sin(tpiq*rad(ir))
998 : 133100821 : pspwk(ir,2,2)=cos(tpiq*rad(ir))
999 : : end do
1000 : :
1001 : : else
1002 : :
1003 : : ! Here the fastest way, iteratively
1004 [ + + ]: 5189172585 : do ir=2,irmax
1005 : 5184065628 : prevsin=pspwk(ir,1,2)
1006 : 5184065628 : prevcos=pspwk(ir,2,2)
1007 : 5184065628 : pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1)
1008 : 5189172585 : 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 : 5239695 : 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 : : !!***
|
