Branch data Line data Source code
1 : : !!****m* ABINIT/m_psptk
2 : : !! NAME
3 : : !! m_psptk
4 : : !!
5 : : !! FUNCTION
6 : : !! This module collects low-level procedures used by the other psp modules
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (XG, DCA, MM, DRH, FrD, GZ, AF)
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_psptk
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_splines
28 : :
29 : : use m_numeric_tools, only : ctrap
30 : : use m_special_funcs, only : sbf8
31 : :
32 : : implicit none
33 : :
34 : : private
35 : : !!***
36 : :
37 : : public :: psp1cc
38 : : public :: psp5lo
39 : : public :: psp5nl
40 : : public :: psp8lo
41 : : public :: psp8nl
42 : : public :: cc_derivatives
43 : : !!***
44 : :
45 : : contains
46 : : !!***
47 : :
48 : : !!****f* ABINIT/psp1cc
49 : : !! NAME
50 : : !! psp1cc
51 : : !!
52 : : !! FUNCTION
53 : : !! Compute the core charge density, for use in the XC core
54 : : !! correction, following the function definition valid
55 : : !! for the format 1 and 5 of pseudopotentials.
56 : : !! WARNING : the fifth derivate is actually set to zero
57 : : !!
58 : : !! INPUTS
59 : : !! fchrg=magnitude of the core charge correction
60 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
61 : : !
62 : : !! OUTPUT
63 : : !! xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
64 : : !!
65 : : !! NOTES
66 : : !! This is a revised expression for core density (5 Nov 1992) :
67 : : !! density(r)=fchrg*gg(xx)
68 : : !! with
69 : : !! $ gg(xx)=(\frac{\sin(2\pi xx)}{(2\pi xx)(1-4 xx^2)(1-xx^2)})^2 $
70 : : !! and
71 : : !! $ xx=\frac{r}{rchrg}=\frac{r}{xcccrc/3.0d0}=3*\frac{r}{xcccrc}=3*yy $
72 : : !!
73 : : !! Code for gg(xx), gp(xx), and gpp(xx) has been tested by numerical
74 : : !! derivatives--looks ok. gpp(x) should still be rewritten.
75 : : !! The argument of xccc1d is assumed to be normalized, and to vary
76 : : !! from yy=0 to 1 (from r=0 to r=xcccrc, or from xx=0 to 3)
77 : : !! Thus :
78 : : !!{{\ \begin{equation}
79 : : !! xccc1d(yy)=fchrg*[\frac{\sin(2*\pi*(3yy))}
80 : : !! {(6*\pi*(3yy))(1-4*(3yy)^2)(1-(3yy)^2)}]^2
81 : : !!\end{equation} }}
82 : : !!
83 : : !! WARNINGS
84 : : !! Warning: the fifth derivative is not yet delivered.
85 : : !!
86 : : !! SOURCE
87 : :
88 : 589 : subroutine psp1cc(fchrg,n1xccc,xccc1d)
89 : :
90 : : !Arguments ------------------------------------
91 : : !scalars
92 : : integer,intent(in) :: n1xccc
93 : : real(dp),intent(in) :: fchrg
94 : : !arrays
95 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
96 : :
97 : : !Local variables-------------------------------
98 : : !scalars
99 : : integer :: i1xccc,ider
100 : : real(dp) :: der1,dern,factor,gg1cc_xx,gp1cc_xx,gpp1cc_xx,xx
101 : : character(len=500) :: message
102 : : !arrays
103 : 589 : real(dp),allocatable :: ff(:),ff2(:),work(:),yy(:)
104 : : ! *************************************************************************
105 : :
106 [ - + + - : 1767 : ABI_MALLOC(ff,(n1xccc))
+ - + - -
+ - + ]
107 [ - + + - : 1178 : ABI_MALLOC(ff2,(n1xccc))
- + - + ]
108 [ - + + - : 1178 : ABI_MALLOC(work,(n1xccc))
- + - + ]
109 [ - + + - : 1178 : ABI_MALLOC(yy,(n1xccc))
- + - + ]
110 : :
111 [ + - ]: 589 : if(n1xccc > 1)then
112 : 589 : factor=one/dble(n1xccc-1)
113 [ + + ]: 1473678 : do i1xccc=1,n1xccc
114 : 1473678 : yy(i1xccc)=(i1xccc-1)*factor
115 : : end do
116 : : else
117 : 0 : write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc
118 : 0 : ABI_BUG(message)
119 : : end if
120 : :
121 : : !Initialization, to avoid some problem with some compilers
122 [ + + + + ]: 7657 : xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
123 : :
124 : : !Take care of each derivative separately
125 [ + + ]: 2356 : do ider=0,2
126 : :
127 [ + + ]: 1767 : if(ider==0)then
128 : : ! Generate spline fitting for the function gg
129 [ + + ]: 1473678 : do i1xccc=1,n1xccc
130 : 1473089 : xx=three*yy(i1xccc)
131 : 1473089 : call gg1cc(gg1cc_xx,xx)
132 : 1473678 : ff(i1xccc)=fchrg*gg1cc_xx
133 : : end do
134 : : ! Complete with derivatives at end points
135 : 589 : der1=zero
136 : 589 : call gp1cc(gp1cc_xx,three)
137 : 589 : dern=three*fchrg*gp1cc_xx
138 [ + + ]: 1178 : else if(ider==1)then
139 : : ! Generate spline fitting for the function gp
140 [ + + ]: 1473678 : do i1xccc=1,n1xccc
141 : 1473089 : xx=three*yy(i1xccc)
142 : 1473089 : call gp1cc(gp1cc_xx,xx)
143 : 1473678 : ff(i1xccc)=three*fchrg*gp1cc_xx
144 : : end do
145 : : ! Complete with derivatives at end points, already estimated
146 : 589 : der1=xccc1d(1,ider+2)
147 : 589 : dern=xccc1d(n1xccc,ider+2)
148 : : else if(ider==2)then
149 : : ! Generate spline fitting for the function gpp
150 : : ! (note : the function gpp has already been estimated, for the spline
151 : : ! fitting of the function gg, but it is replaced here by the more
152 : : ! accurate analytic derivative)
153 [ + + ]: 1473678 : do i1xccc=1,n1xccc
154 : 1473089 : xx=three*yy(i1xccc)
155 : 1473089 : call gpp1cc(gpp1cc_xx,xx)
156 : 1473678 : ff(i1xccc)=9.0_dp*fchrg*gpp1cc_xx
157 : : end do
158 : : ! Complete with derivatives of end points
159 : 589 : der1=xccc1d(1,ider+2)
160 : 589 : dern=xccc1d(n1xccc,ider+2)
161 : : end if
162 : :
163 : : ! Produce second derivative numerically, for use with splines
164 : 1767 : call spline(yy,ff,n1xccc,der1,dern,ff2)
165 [ + + ]: 4421034 : xccc1d(:,ider+1)=ff(:)
166 [ + + ]: 4421623 : xccc1d(:,ider+3)=ff2(:)
167 : : end do
168 : :
169 [ + + ]: 1473678 : xccc1d(:,6)=zero
170 : :
171 : : !DEBUG
172 : : !write(std_out,*)' psp1cc : output of core charge density and derivatives '
173 : : !write(std_out,*)' yy gg gp '
174 : : !do i1xccc=1,n1xccc
175 : : !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
176 : : !end do
177 : : !write(std_out,*)' yy gpp gg2 '
178 : : !do i1xccc=1,n1xccc
179 : : !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
180 : : !end do
181 : : !write(std_out,*)' yy gp2 gpp2 '
182 : : !do i1xccc=1,n1xccc
183 : : !write(std_out,'(3es14.6)' ) yy(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
184 : : !end do
185 : : !write(std_out,*)' psp1cc : debug done, stop '
186 : : !stop
187 : : !ENDDEBUG
188 : :
189 [ - + ]: 589 : ABI_FREE(ff)
190 [ - + ]: 589 : ABI_FREE(ff2)
191 [ - + ]: 589 : ABI_FREE(work)
192 [ - + ]: 589 : ABI_FREE(yy)
193 : :
194 : 589 : end subroutine psp1cc
195 : : !!***
196 : :
197 : : !!****f* ABINIT/gg1cc
198 : : !! NAME
199 : : !! gg1cc
200 : : !!
201 : : !! FUNCTION
202 : : !! gg1cc_xx=$(\frac{\sin(2\pi xx)}{(2\pi xx)(1-4xx^2)(1-xx^2)})^2$
203 : : !!
204 : : !! INPUTS
205 : : !! xx= abscisse to which gg1cc_xx is calculated
206 : : !!
207 : : !! OUTPUT
208 : : !! gg1cc_xx= gg1cc_x(xx)
209 : : !!
210 : : !! SOURCE
211 : :
212 : 1473089 : subroutine gg1cc(gg1cc_xx,xx)
213 : :
214 : : !Arguments ------------------------------------
215 : : !scalars
216 : : real(dp),intent(in) :: xx
217 : : real(dp),intent(out) :: gg1cc_xx
218 : :
219 : : !Local variables -------------------------------------------
220 : : !The c s are coefficients for Taylor expansion of the analytic form near xx=0, 1/2, and 1.
221 : : !scalars
222 : : real(dp) :: c21=4.d0/9.d0,c22=-40.d0/27.d0,c23=20.d0/3.d0-16.d0*pi**2/27.d0
223 : : real(dp) :: c24=-4160.d0/243.d0+160.d0*pi**2/81.d0,c31=1.d0/36.d0
224 : : real(dp) :: c32=-25.d0/108.d0,c33=485.d0/432.d0-pi**2/27.d0
225 : : real(dp) :: c34=-4055.d0/972.d0+25.d0*pi**2/81.d0
226 : : ! *************************************************************************
227 : :
228 : : !Cut off beyond 3/gcut=xcccrc
229 [ - + ]: 1473089 : if (xx>3.0d0) then
230 : 0 : gg1cc_xx=0.0d0
231 : : ! Take care of difficult limits near x=0, 1/2, and 1
232 [ + + ]: 1473089 : else if (abs(xx)<=1.d-09) then
233 : 589 : gg1cc_xx=1.d0
234 [ - + ]: 1472500 : else if (abs(xx-0.5d0)<=1.d-04) then
235 : : ! (this limit and next are more troublesome for numerical cancellation)
236 : 0 : gg1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*c24))
237 [ - + ]: 1472500 : else if (abs(xx-1.d0)<=1.d-04) then
238 : 0 : gg1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*c34))
239 : : else
240 : : ! The following is the square of the Fourier transform of a
241 : : ! function built out of two spherical bessel functions in G
242 : : ! space and cut off absolutely beyond gcut
243 : : gg1cc_xx=(sin(2.0d0*pi*xx)/( (2.0d0*pi*xx) * &
244 : 1472500 : & (1.d0-4.0d0*xx**2)*(1.d0-xx**2) ) )**2
245 : : end if
246 : :
247 : 1473089 : end subroutine gg1cc
248 : : !!***
249 : :
250 : : !!****f* ABINIT/gp1cc
251 : : !! NAME
252 : : !! gp1cc
253 : : !!
254 : : !! FUNCTION
255 : : !! Derivative of gg(xx) wrt xx.
256 : : !!
257 : : !! INPUTS
258 : : !! xx=abscisse to which gp1cc_xx is calculated
259 : : !!
260 : : !! OUTPUT
261 : : !! gp1cc_xx=derivative of gg(xx) wrt xx.
262 : : !!
263 : : !! NOTES
264 : : !! $ phi(x) = \frac{\sin(2\pi x)}{(2\pi x)(1-4x^2)(1-x^2)}$
265 : : !! $ gg(x)= phi(x)^2$
266 : : !! $ gp(x)= 2 * phi(x) * phi''(x)$
267 : : !! $ phi''(x)=\frac{\cos(2\pi x)-(1-15x^2+20x^4) phi(x)}{x(1-4x^2)(1-x^2)}$
268 : : !!
269 : : !! SOURCE
270 : :
271 : 1473678 : subroutine gp1cc(gp1cc_xx,xx)
272 : :
273 : : !Arguments ------------------------------------
274 : : !scalars
275 : : real(dp),intent(in) :: xx
276 : : real(dp),intent(out) :: gp1cc_xx
277 : :
278 : : !Local variables -------------------------------------------
279 : : !scalars
280 : : real(dp),parameter :: c11=20.d0-8.d0*pi**2/3.d0
281 : : real(dp),parameter :: c12=268.d0-160.d0/3.d0*pi**2+128.d0/45.d0*pi**4
282 : : real(dp),parameter :: c21=-40.d0/27.d0,c22=40.d0/3.d0-32.d0*pi**2/27.d0
283 : : real(dp),parameter :: c23=-4160.d0/81.d0+160.d0*pi**2/27.d0
284 : : real(dp),parameter :: c24=157712.d0/729.d0-320.d0*pi**2/9.d0+512.d0*pi**4/405.d0
285 : : real(dp),parameter :: c25=-452200.d0/729.d0+83200.d0*pi**2/729.d0-1280.d0*pi**4/243.d0
286 : : real(dp),parameter :: c31=-25.d0/108.d0,c32=485.d0/216.d0-2.d0*pi**2/27.d0
287 : : real(dp),parameter :: c33=-4055.d0/324.d0+25.d0*pi**2/27.d0
288 : : real(dp),parameter :: c34=616697.d0/11664.d0-485.d0*pi**2/81.d0+32.d0*pi**4/405.d0
289 : : real(dp),parameter :: c35=-2933875.d0/15552.d0+20275.d0*pi**2/729.d0-200.d0*pi**4/243.d0
290 : : real(dp),parameter :: two_pim1=1.0d0/two_pi
291 : : real(dp) :: denom,phi,phip
292 : : ! *************************************************************************
293 : :
294 : : !Cut off beyond r=3*xcccrc is already done at the calling level
295 [ + + ]: 1473678 : if (xx>1.001d0) then
296 : : ! The part that follows will be repeated later, but written in this way,
297 : : ! only one "if" condition is tested in most of the cases (1.001 < x < 3.0)
298 : 981863 : denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
299 : 981863 : phi=denom*sin(two_pi*xx)*two_pim1
300 : 981863 : phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
301 : 981863 : gp1cc_xx=2.d0*phi*phip
302 : : ! Handle limits where denominator vanishes
303 [ + + ]: 491815 : else if (abs(xx)<1.d-03) then
304 : 589 : gp1cc_xx=xx*(c11+xx**2*c12)
305 [ + + ]: 491226 : else if (abs(xx-0.5d0)<=1.d-03) then
306 : 1178 : gp1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*(c24+(xx-0.5d0)*c25)))
307 [ + + ]: 490048 : else if (abs(xx-1.d0)<=1.d-03) then
308 : 1178 : gp1cc_xx=c31+(xx-1.0d0)*(c32+(xx-1.0d0)*(c33+(xx-1.0d0)*(c34+(xx-1.0d0)*c35)))
309 : : else
310 : : ! Here is the repeated part ...
311 : 488870 : denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
312 : 488870 : phi=denom*sin(two_pi*xx)*two_pim1
313 : 488870 : phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
314 : 488870 : gp1cc_xx=2.d0*phi*phip
315 : : end if
316 : :
317 : 1473678 : end subroutine gp1cc
318 : : !!***
319 : :
320 : : !!****f* ABINIT/gpp1cc
321 : : !! NAME
322 : : !! gpp1cc
323 : : !!
324 : : !! FUNCTION
325 : : !! Second derivative of gg wrt xx.
326 : : !!
327 : : !! INPUTS
328 : : !! xx= abscisse to which gpp1cc_xx is calculated
329 : : !!
330 : : !! OUTPUT
331 : : !! gpp1cc_xx=second derivative of gg wrt xx.
332 : : !!
333 : : !! SOURCE
334 : :
335 : 1473089 : subroutine gpp1cc(gpp1cc_xx,xx)
336 : :
337 : : !Arguments ------------------------------------
338 : : !scalars
339 : : real(dp),intent(in) :: xx
340 : : real(dp),intent(out) :: gpp1cc_xx
341 : :
342 : : !Local variables -------------------------------------------
343 : : !scalars
344 : : real(dp),parameter :: c1=20.d0-8.d0*pi**2/3.d00
345 : : real(dp),parameter :: c2=40.d0/3.d0-32.d0*pi**2/27.d0
346 : : real(dp),parameter :: c3=-8320.d0/81.d0+320.d0*pi**2/27.d0
347 : : real(dp),parameter :: c4=157712.d0/243.d0-320.d0*pi**2/3.d0+512.d0*pi**4/135.d0
348 : : real(dp),parameter :: c5=-18088.d2/729.d0+3328.d2*pi**2/729.d0-5120.d0*pi**4/243.d0
349 : : real(dp),parameter :: c6=485.d0/216.d0-2.d0*pi**2/27.d0
350 : : real(dp),parameter :: c7=-4055.d0/162.d0+50.d0*pi**2/27.d0
351 : : real(dp),parameter :: c8=616697.d0/3888.d0-485.d0*pi**2/27.d0+32.d0*pi**4/135.d0
352 : : real(dp),parameter :: c9=-2933875.d0/3888.d0+81100.d0*pi**2/729.d0-800.d0*pi**4/243.d0
353 : : real(dp) :: t1,t10,t100,t11,t12,t120,t121,t122,t127,t138,t14,t140,t15,t152
354 : : real(dp) :: t157,t16,t160,t17,t174,t175,t18,t19,t2,t20,t21,t23,t24,t3,t31,t33
355 : : real(dp) :: t34,t4,t41,t42,t44,t45,t46,t5,t54,t55,t56,t57,t6,t62,t64,t65,t7
356 : : real(dp) :: t72,t78,t79,t8,t85,t9,t93
357 : : ! *************************************************************************
358 : :
359 [ - + ]: 1473089 : if (xx>3.0d0) then
360 : : ! Cut off beyond 3/gcut=3*xcccrc
361 : 0 : gpp1cc_xx=0.0d0
362 : : ! Take care of difficult limits near xx=0, 1/2, and 1
363 [ + + ]: 1473089 : else if (abs(xx)<=1.d-09) then
364 : 589 : gpp1cc_xx=c1
365 [ - + ]: 1472500 : else if (abs(xx-0.5d0)<=1.d-04) then
366 : : ! (this limit and next are more troublesome for numerical cancellation)
367 : 0 : gpp1cc_xx=c2+(xx-0.5d0)*(c3+(xx-0.5d0)*(c4+(xx-0.5d0)*c5))
368 [ - + ]: 1472500 : else if (abs(xx-1.d0)<=1.d-04) then
369 : 0 : gpp1cc_xx=c6+(xx-1.0d0)*(c7+(xx-1.0d0)*(c8+(xx-1.0d0)*c9))
370 : : else
371 : :
372 : : ! Should fix up this Maple fortran later
373 : 1472500 : t1 = xx**2
374 : 1472500 : t2 = 1/t1
375 : 1472500 : t3 = 1/Pi
376 : 1472500 : t4 = 2*xx
377 : 1472500 : t5 = t4-1
378 : 1472500 : t6 = t5**2
379 : 1472500 : t7 = 1/t6
380 : 1472500 : t8 = t4+1
381 : 1472500 : t9 = t8**2
382 : 1472500 : t10 = 1/t9
383 : 1472500 : t11 = xx-1
384 : 1472500 : t12 = t11**2
385 : 1472500 : t14 = 1/t12/t11
386 : 1472500 : t15 = xx+1
387 : 1472500 : t16 = t15**2
388 : 1472500 : t17 = 1/t16
389 : 1472500 : t18 = Pi*xx
390 : 1472500 : t19 = sin(t18)
391 : 1472500 : t20 = cos(t18)
392 : 1472500 : t21 = t20**2
393 : 1472500 : t23 = t19*t21*t20
394 : 1472500 : t24 = t17*t23
395 : 1472500 : t31 = t19**2
396 : 1472500 : t33 = t31*t19*t20
397 : 1472500 : t34 = t17*t33
398 : 1472500 : t41 = Pi**2
399 : 1472500 : t42 = 1/t41
400 : 1472500 : t44 = 1/t16/t15
401 : 1472500 : t45 = t31*t21
402 : 1472500 : t46 = t44*t45
403 : 1472500 : t54 = 1/t1/xx
404 : 1472500 : t55 = 1/t12
405 : 1472500 : t56 = t55*t46
406 : 1472500 : t57 = t10*t56
407 : 1472500 : t62 = t9**2
408 : 1472500 : t64 = t17*t45
409 : 1472500 : t65 = t55*t64
410 : 1472500 : t72 = 1/t9/t8
411 : 1472500 : t78 = t14*t64
412 : 1472500 : t79 = t10*t78
413 : 1472500 : t85 = t12**2
414 : 1472500 : t93 = t21**2
415 : 1472500 : t100 = t31**2
416 : 1472500 : t120 = 1/t6/t5
417 : 1472500 : t121 = t55*t34
418 : 1472500 : t122 = t10*t121
419 : 1472500 : t127 = t16**2
420 : 1472500 : t138 = t6**2
421 : 1472500 : t140 = t10*t65
422 : 1472500 : t152 = t72*t65
423 : 1472500 : t157 = t7*t140
424 : 1472500 : t160 = t1**2
425 : 1472500 : t174 = t55*t24
426 : 1472500 : t175 = t10*t174
427 : : gpp1cc_xx = 8*t2*t3*t7*t10*t14*t34+8*t2*t42*t7*t10*t14*t46&
428 : : & -8*t2*t3*t7*t10*t14*t24+8*t2*t3*t7*t10*t55*t44*t33+&
429 : : & 6*t2*t42*t7*t10*t55/t127*t45+24*t2*t42/t138*t140+&
430 : : & 16*t54*t42*t120*t140+16*t2*t3*t120*t122+16*t2&
431 : : & *t42*t7*t72*t78-8*t2*t3*t7*t10*t55*t44*t23-8*t54*t3*t7*t175&
432 : : & +2*t2*t7*t10*t55*t17*t100+2*t2*t7*t10*t55*t17*t93+&
433 : : & 8*t54*t42*t7*t79+16*t2*t42*t7*t72*t56+6*t2*t42*t7*t10/t85&
434 : : & *t64+24*t2*t42*t7/t62*t65+8*t54*t42*t7*t57-&
435 : : & 16*t2*t3*t7*t72*t174+8*t54*t3*t7*t122-16*t2*t3*t120*t175&
436 : : & +16*t2*t42*t120*t79+16*t2*t42*t120*t57+16*t54*t42*t7*t152+&
437 : : & 32*t2*t42*t120*t152+16*t2*t3*t7*t72*t121-12*t2*t157+&
438 : 1472500 : & 6/t160*t42*t157
439 : : end if
440 : :
441 : 1473089 : end subroutine gpp1cc
442 : : !!***
443 : :
444 : : !!****f* ABINIT/psp5lo
445 : : !! NAME
446 : : !! psp5lo
447 : : !!
448 : : !! FUNCTION
449 : : !! Compute sine transform to transform from V(r) to q^2 V(q).
450 : : !! Computes integrals on logarithmic grid using related uniform
451 : : !! grid in exponent and corrected trapezoidal integration.
452 : : !!
453 : : !! INPUTS
454 : : !! al=spacing in exponent for radial atomic grid.
455 : : !! mmax=number of radial r grid points (logarithmic atomic grid).
456 : : !! mqgrid=number of grid points in q from 0 to qmax.
457 : : !! qgrid(mqgrid)=q grid values (bohr**-1).
458 : : !! rad(mmax)=r grid values (bohr).
459 : : !! vloc(mmax)=V(r) on radial grid.
460 : : !! zion=nominal valence charge of atom.
461 : : !!
462 : : !! OUTPUT
463 : : !! epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
464 : : !!{{\\ \begin{equation}
465 : : !! q2vq(mqgrid)
466 : : !! =q^2 V(q)
467 : : !! = -\frac{Zv}{\pi}
468 : : !! + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
469 : : !!\end{equation} }}
470 : : !! yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax
471 : : !! (needed for spline fitter).
472 : : !!
473 : : !! SOURCE
474 : :
475 : 324 : subroutine psp5lo(al,epsatm,mmax,mqgrid,qgrid,q2vq,rad,&
476 : 324 : & vloc,yp1,ypn,zion)
477 : :
478 : : !Arguments----------------------------------------------------------
479 : : !scalars
480 : : integer,intent(in) :: mmax,mqgrid
481 : : real(dp),intent(in) :: al,zion
482 : : real(dp),intent(out) :: epsatm,yp1,ypn
483 : : !arrays
484 : : real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
485 : : real(dp),intent(out) :: q2vq(mqgrid)
486 : :
487 : : !Local variables-------------------------------
488 : : !scalars
489 : : integer :: iq,ir
490 : : real(dp),parameter :: scale=10.0d0
491 : : real(dp) :: arg,result,rmtoin,test,ztor1
492 : : !arrays
493 : 324 : real(dp),allocatable :: work(:)
494 : : ! *************************************************************************
495 : :
496 [ - + + - : 972 : ABI_MALLOC(work,(mmax))
+ - + - -
+ - + ]
497 : :
498 : : !Do q=0 separately (compute epsatm)
499 : : !Do integral from 0 to r1
500 : 324 : ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
501 : :
502 : : !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
503 : : !with extra factor of r to convert to uniform grid in exponent
504 [ + + ]: 159840 : do ir=1,mmax
505 : : ! First handle tail region
506 : 159516 : test=vloc(ir)+zion/rad(ir)
507 : : ! DEBUG
508 : : ! write(std_out,*)ir,rad(ir),test
509 : : ! ENDDEBUG
510 : : ! Ignore small contributions, or impose a cut-off in the case
511 : : ! the pseudopotential data are in single precision.
512 : : ! (it is indeed expected that vloc is very close to zero beyond 20,
513 : : ! so a value larger than 2.0d-8 is considered anomalous)
514 [ + + + + : 159840 : if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then
+ + ]
515 : 3064 : work(ir)=zero
516 : : else
517 : 156452 : work(ir)=(rad(ir)*rad(ir))*(rad(ir)*vloc(ir)+zion)
518 : : end if
519 : : end do
520 : : !write(std_out,*)' psp5lo : stop '; stop
521 : :
522 : : !Do integral from r(1) to r(max)
523 : 324 : call ctrap(mmax,work,al,result)
524 : : !Do integral from r(mmax) to infinity
525 : : !compute decay length lambda at r(mmax)
526 : : !$\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ &
527 : : !$(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$
528 : : !rmtoin=$(rad(mmax)*vloc(mmax)+zion)*(rad(mmax)+1.d0/\lambda)/\lambda$
529 : : !Due to inability to fit exponential decay to r*V(r)+Zv
530 : : !in tail, NO TAIL CORRECTION IS APPLIED
531 : : !(numerical trouble might be removed if atomic code is
532 : : !cleaned up in tail region)
533 : 324 : rmtoin=0.0d0
534 : :
535 : 324 : epsatm=4.d0*pi*(result+ztor1+rmtoin)
536 : :
537 : 324 : q2vq(1)=-zion/pi
538 : :
539 : : !Loop over q values
540 [ + + ]: 970438 : do iq=2,mqgrid
541 : 970114 : arg=2.d0*pi*qgrid(iq)
542 : : ! ztor1=$ -Zv/\pi+2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$
543 : : ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion)* &
544 : 970114 : & cos(arg*rad(1)) )/pi
545 : :
546 : : ! set up integrand
547 [ + + ]: 476319419 : do ir=1,mmax
548 : 475349305 : test=vloc(ir)+zion/rad(ir)
549 : : ! Ignore contributions within decade of machine precision
550 [ + + ]: 476319419 : if ((scale+abs(test)).eq.scale) then
551 : 20359150 : work(ir)=zero
552 : : else
553 : 454990155 : work(ir)=rad(ir)*sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion)
554 : : end if
555 : : end do
556 : : ! do integral from r(1) to r(mmax)
557 : 970114 : call ctrap(mmax,work,al,result)
558 : :
559 : : ! do integral from r(mmax) to infinity
560 : : ! rmtoin=(r(mmax)*vr(mmax)+zion)*(lambda*sin(arg*r(mmax))+
561 : : ! arg*cos(arg*r(mmax)))/(arg**2+lambda**2)
562 : : ! See comment above; no tail correction
563 : 970114 : rmtoin=0.0d0
564 : :
565 : : ! store q^2 v(q)
566 : 970438 : q2vq(iq)=ztor1+2.d0*qgrid(iq)*(result+rmtoin)
567 : :
568 : : end do
569 : :
570 : : !Compute derivatives of q^2 v(q) at ends of interval
571 : 324 : yp1=0.0d0
572 : : !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
573 : : !integral from 0 to r1
574 : 324 : arg=2.0d0*pi*qgrid(mqgrid)
575 : 324 : ztor1=zion*rad(1)*sin(arg*rad(1))
576 : : ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + &
577 : 324 : & (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1))
578 : : !integral from r(mmax) to infinity is overkill; ignore
579 : : !set up integrand
580 [ + + ]: 159840 : do ir=1,mmax
581 : 159516 : test=vloc(ir)+zion/rad(ir)
582 : : ! Ignore contributions within decade of machine precision
583 [ + + ]: 159840 : if ((scale+abs(test)).eq.scale) then
584 : 6785 : work(ir)=0.0d0
585 : : else
586 : : work(ir)=rad(ir)*(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * &
587 : 152731 : & (rad(ir)*vloc(ir)+zion)
588 : : end if
589 : : end do
590 : 324 : call ctrap(mmax,work,al,result)
591 : 324 : ypn=2.0d0 * (ztor1 + result)
592 : :
593 [ - + ]: 324 : ABI_FREE(work)
594 : :
595 : 648 : end subroutine psp5lo
596 : : !!***
597 : :
598 : : !!****f* ABINIT/psp5nl
599 : : !! NAME
600 : : !! psp5nl
601 : : !!
602 : : !! FUNCTION
603 : : !! Make Kleinman-Bylander form factors f_l(q) for each l from 0 to lmax.
604 : : !! Vloc is assumed local potential.
605 : : !!
606 : : !! INPUTS
607 : : !! al=grid spacing in exponent for radial grid
608 : : !! lmax=maximum ang momentum for which nonlocal form factor is desired.
609 : : !! Usually lmax=1, sometimes = 0 (e.g. for oxygen); lmax <= 2 allowed.
610 : : !! mmax=number of radial grid points for atomic grid
611 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
612 : : !! mqgrid=number of grid points for q grid
613 : : !! qgrid(mqgrid)=values at which form factors are returned
614 : : !! rad(mmax)=radial grid values
615 : : !! vloc(mmax)=local pseudopotential on radial grid
616 : : !! vpspll(mmax,3)=nonlocal pseudopotentials for each l on radial grid
617 : : !! wfll(mmax,3)=reference state wavefunctions on radial grid mmax and mqgrid
618 : : !!
619 : : !! OUTPUT
620 : : !! ekb(mpsang)=Kleinman-Bylander energy,
621 : : !! {{\\ \begin{equation}
622 : : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
623 : : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
624 : : !! \end{equation} }}
625 : : !! for each l
626 : : !! ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and
627 : : !! second derivative from spline fit for each angular momentum
628 : : !!
629 : : !! NOTES
630 : : !! u_l(r) is reference state wavefunction (input as wf);
631 : : !! j_l(q) is a spherical Bessel function;
632 : : !! dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
633 : : !! f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
634 : : !! where dvms = $\int_0^{rmax} [(u_l(r) dV_l(r))^2 dr]$ is the mean
635 : : !! square value of the nonlocal correction for angular momentum l.
636 : : !! Xavier Gonze s E_KB = $ dvms/\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]$.
637 : : !! This is the eigenvalue of the Kleinman-Bylander operator and sets
638 : : !! the energy scale of the nonlocal psp corrections.
639 : : !!
640 : : !! SOURCE
641 : :
642 : 322 : subroutine psp5nl(al,ekb,ffspl,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll)
643 : :
644 : : !Arguments ------------------------------------
645 : : !scalars
646 : : real(dp),intent(in) :: al
647 : : integer,intent(in) :: lmax,mmax,mpsang,mqgrid
648 : : !arrays
649 : : real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax),vpspll(mmax,mpsang)
650 : : real(dp),intent(in) :: wfll(mmax,mpsang)
651 : : real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang)
652 : :
653 : : !Local variables-------------------------------
654 : : !scalars
655 : : integer,parameter :: dpsang=5
656 : : integer :: iq,ir,lp1
657 : : real(dp) :: arg,bessel,dvwf,qr,result,yp1,ypn,ztor1
658 : : character(len=500) :: message
659 : : !arrays
660 : : real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang)
661 : 322 : real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:)
662 : : !*************************************************************************
663 : :
664 : : !l=0,1,2 and 3 spherical Bessel functions
665 : : !The accuracy of the bes1, bes2, bes3 functions for small arguments
666 : : !may be insufficient. In the present version
667 : : !of the routines, some care is taken with the value of the argument.
668 : : !If smaller than 1.d-3, a two terms
669 : : !Taylor series expansion is prefered.
670 : : ! bes0(arg)=sin(arg)/arg
671 : : ! bes1(arg)=(sin(arg)-arg*cos(arg))/arg**2
672 : : ! bes2(arg)=( (3.0d0-arg**2)*sin(arg)-&
673 : : !& 3.0d0*arg*cos(arg) ) /arg**3
674 : :
675 : : ! bes3(arg)=(15.d0*sin(arg)-15.d0*arg*cos(arg) &
676 : : !& -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4
677 : :
678 : : !Zero out Kleinman-Bylander energies ekb
679 [ + + ]: 1283 : ekb(:)=0.0d0
680 : :
681 [ - + + - : 966 : ABI_MALLOC(work1,(mmax))
+ - + - -
+ - + ]
682 [ - + + - : 644 : ABI_MALLOC(work2,(mmax))
- + - + ]
683 [ - + + - : 644 : ABI_MALLOC(work3,(mmax))
- + - + ]
684 [ - + + - : 644 : ABI_MALLOC(work4,(mmax))
- + - + ]
685 : :
686 : : !Allow for no nonlocal correction (lmax=-1)
687 [ + - ]: 322 : if (lmax/=-1) then
688 : :
689 : : ! Check that lmax is within allowed range
690 [ - + ]: 322 : if (lmax<0.or.lmax>3) then
691 : : write(message, '(a,i12,a,a,a,a,a,a,a)' )&
692 : 0 : & 'lmax=',lmax,' is not an allowed value.',ch10,&
693 : 0 : & 'Allowed values are -1 for no nonlocal correction or else',ch10,&
694 : 0 : & '0, 1,2 or 3 for maximum l nonlocal correction.',ch10,&
695 : 0 : & 'Action: check the input atomic psp data file for lmax.'
696 : 0 : ABI_ERROR(message)
697 : : end if
698 : :
699 : : ! Compute normalizing integrals eta=<dV> and mean square
700 : : ! nonlocal psp correction dvms=<dV^2>
701 : : ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction
702 [ + + ]: 1277 : do lp1=1,lmax+1
703 : :
704 : : ! integral from 0 to r1
705 : 955 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)
706 : 955 : ztor1=(wfll(1,lp1)*dvwf)*rad(1)/dble(2*(lp1-1)+3)
707 : : ! integrand for r1 to r(mmax) (incl extra factor of r)
708 [ + + ]: 469215 : do ir=1,mmax
709 : 468260 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
710 : 469215 : work1(ir)=rad(ir)*(wfll(ir,lp1)*dvwf)
711 : : end do
712 : : ! do integral by corrected trapezoidal integration
713 : 955 : call ctrap(mmax,work1,al,result)
714 : 955 : eta(lp1)=ztor1+result
715 : :
716 : 955 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)
717 : 955 : ztor1=dvwf**2*rad(1)/dble(2*(lp1-1)+3)
718 [ + + ]: 469215 : do ir=1,mmax
719 : 468260 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)
720 : 469215 : work1(ir)=rad(ir)*(dvwf**2)
721 : : end do
722 : 955 : call ctrap(mmax,work1,al,result)
723 : 955 : dvms(lp1)=ztor1+result
724 : :
725 : : ! DEBUG
726 : : ! Compute the norm of wfll
727 : : ! wf=wfll(1,lp1)
728 : : ! ztor1=wf**2*rad(1)/dble(2*(lp1-1)+3)
729 : : ! do ir=1,mmax
730 : : ! wf=wfll(ir,lp1)
731 : : ! work1(ir)=rad(ir)*(wf**2)
732 : : ! end do
733 : : ! call ctrap(mmax,work1,al,result)
734 : : ! norm=ztor1+result
735 : : ! write(std_out,*)' lp1, norm',lp1,norm
736 : : ! ENDDEBUG
737 : :
738 : : ! If dvms is not 0 for any given angular momentum l,
739 : : ! compute Xavier Gonze's definition of the Kleinman-Bylander
740 : : ! energy E_KB = dvms/eta. In this case also renormalize
741 : : ! the projection operator to u_KB(r)=$u_l(r)*dV(r)/\sqrt{dvms}$.
742 : : ! This means dvwf gets multiplied by the normalization factor
743 : : ! "renorm"=$1/\sqrt{dvms}$ as seen below.
744 [ + + ]: 2232 : if (dvms(lp1)/=0.0d0) then
745 : 633 : ekb(lp1)=dvms(lp1)/eta(lp1)
746 : 633 : renorm(lp1)=1.0d0/sqrt(dvms(lp1))
747 : : ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze)
748 : : ckb(lp1)=eta(lp1)/sqrt(dvms(lp1))
749 : : else
750 : 322 : ekb(lp1)=0.0d0
751 : : end if
752 : :
753 : : end do
754 : :
755 : : ! l=0 form factor if ekb(1) not 0 (lmax always at least 0)
756 [ + + ]: 322 : if (ekb(1)/=0.0d0) then
757 : :
758 : : ! do q=0 separately
759 : 280 : lp1=1
760 : : ! 0 to r1 integral
761 : 280 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
762 : 280 : ztor1=(rad(1)*dvwf)*rad(1)/3.0d0
763 : : ! integrand
764 [ + + ]: 128567 : do ir=1,mmax
765 : 128287 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
766 : 128567 : work1(ir)=rad(ir)*(rad(ir)*dvwf)
767 : : end do
768 : 280 : call ctrap(mmax,work1,al,result)
769 : 280 : ffspl(1,1,1)=ztor1+result
770 : :
771 : : ! do rest of q points
772 [ + + ]: 837845 : do iq=2,mqgrid
773 : 837565 : arg=two_pi*qgrid(iq)
774 : : ! 0 to r1 integral
775 : 837565 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
776 : 837565 : ztor1=(bes0_psp5(arg*rad(1))*rad(1)*dvwf)*rad(1)/3.0d0
777 : : ! integrand
778 [ + + ]: 382181450 : do ir=1,mmax
779 : 381343885 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
780 : 837565 : work1(ir)=rad(ir)*(rad(ir)*bes0_psp5(arg*rad(ir))*dvwf)
781 : : end do
782 : 837565 : call ctrap(mmax,work1,al,result)
783 : 837845 : ffspl(iq,1,1)=ztor1+result
784 : : end do
785 : :
786 : : ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
787 : : ! yp1=0 for l=0
788 : 280 : yp1=0.0d0
789 : : ! ypn=$ \int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$
790 : 280 : arg=two_pi*qgrid(mqgrid)
791 : 280 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
792 : 280 : qr=arg*rad(1)
793 [ + + ]: 280 : if(qr<1.d-3)then
794 : 92 : bessel=(10.d0-qr*qr)*qr/30.0d0
795 : : else
796 : 188 : bessel=bes1_psp5(qr)
797 : : end if
798 : : ! ztor1=(-bes1(arg*rad(1))*two_pi*rad(1)*r(1)*dvwf)*rad(1)/5.0d0
799 : 280 : ztor1=(-bessel*two_pi*rad(1)*rad(1)*dvwf)*rad(1)/5.0d0
800 [ + + ]: 128567 : do ir=1,mmax
801 : 128287 : qr=arg*rad(ir)
802 [ + + ]: 128287 : if(qr<1.d-3)then
803 : 7802 : bessel=(10.d0-qr*qr)*qr/30.0d0
804 : : else
805 : 120485 : bessel=bes1_psp5(qr)
806 : : end if
807 : 128287 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
808 : : ! work(ir)=rad(ir)*(-bes1(arg*rad(ir))*two_pi*rad(ir)*rad(ir)*dvwf)
809 : 128567 : work1(ir)=rad(ir)*(-bessel*two_pi*rad(ir)*rad(ir)*dvwf)
810 : : end do
811 : 280 : call ctrap(mmax,work1,al,result)
812 : 280 : ypn=ztor1+result
813 : :
814 : : ! Fit spline to get second derivatives by spline fit
815 : 280 : call spline(qgrid,ffspl(1,1,1),mqgrid,yp1,ypn,ffspl(1,2,1))
816 : :
817 : : else
818 : : ! or else put nonlocal correction at l=0 to 0
819 [ + + + + ]: 253308 : ffspl(:,:,1)=0.0d0
820 : : end if
821 : :
822 : : ! Finished if lmax=0 (highest nonlocal correction)
823 : : ! Do l=1 form factor if ekb(2) not 0 and lmax>=1
824 [ + - ]: 322 : if (lmax>0)then
825 [ + + ]: 322 : if(ekb(2)/=0.0d0) then
826 : :
827 : 283 : lp1=2
828 : : ! do q=0 separately: f_1(q=0) vanishes !
829 : 283 : ffspl(1,1,2)=0.0d0
830 : :
831 : : ! do rest of q points
832 [ + + ]: 848332 : do iq=2,mqgrid
833 : 848049 : arg=two_pi*qgrid(iq)
834 : 848049 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
835 : 848049 : qr=arg*rad(1)
836 [ + + ]: 848049 : if(qr<1.d-3)then
837 : 588535 : bessel=(10.d0-qr*qr)*qr/30.0d0
838 : : else
839 : 259514 : bessel=bes1_psp5(qr)
840 : : end if
841 : : ! ztor1=(bes1(arg*rad(1))*rad(1)*dvwf)*rad(1)/5.0d0
842 : 848049 : ztor1=(bessel*rad(1)*dvwf)*rad(1)/5.0d0
843 : :
844 [ + + ]: 401004969 : do ir=1,mmax
845 : 400156920 : qr=arg*rad(ir)
846 [ + + ]: 400156920 : if(qr<1.d-3)then
847 : 53985365 : bessel=(10.d0-qr*qr)*qr/30.0d0
848 : : else
849 : 346171555 : bessel=bes1_psp5(qr)
850 : : end if
851 : 400156920 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
852 : 401004969 : work2(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
853 : : end do
854 : :
855 : 848049 : call ctrap(mmax,work2,al,result)
856 : 848332 : ffspl(iq,1,2)=ztor1+result
857 : : end do
858 : :
859 : : ! Compute yp1,ypn for l=1
860 : : ! yp1=$\displaystyle \int [2\pi r^2 wf(r) dV(r)]/3$
861 : 283 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
862 : 283 : ztor1=((two_pi*rad(1)**2)*dvwf)*rad(1)/(3.0d0*5.0d0)
863 [ + + ]: 134187 : do ir=1,mmax
864 : 133904 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
865 : 134187 : work2(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf/3.0d0)
866 : : end do
867 : 283 : call ctrap(mmax,work2,al,result)
868 : 283 : yp1=ztor1+result
869 : : ! ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$
870 : : ! where x=2 Pi qgrid(mqgrid) r
871 : 283 : arg=two_pi*qgrid(mqgrid)
872 : 283 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
873 : 283 : qr=arg*rad(1)
874 [ + + ]: 283 : if(qr<1.d-3)then
875 : 96 : bessel=(10.d0-3.0d0*qr*qr)/30.0d0
876 : : else
877 : 187 : bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr
878 : : end if
879 : : ! ztor1=( (two_pi*rad(1)**2)*dvwf* (bes0(arg*rad(1))-
880 : : ! 2.0d0*bes1(arg*rad(1))/(arg*rad(1))) ) * rad(1)/5.0d0
881 : 283 : ztor1=( (two_pi*rad(1)**2)*dvwf*bessel)* rad(1)/5.0d0
882 : :
883 [ + + ]: 134187 : do ir=1,mmax
884 : 133904 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
885 : 133904 : qr=arg*rad(ir)
886 [ + + ]: 133904 : if(qr<1.d-3)then
887 : 11581 : bessel=(10.d0-3.0d0*qr*qr)/30.0d0
888 : : else
889 : 122323 : bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr
890 : : end if
891 : : ! work(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
892 : : ! (bes0(arg*rad(ir))-2.d0*bes1(arg*rad(ir))/(arg*rad(ir))) )
893 : 134187 : work2(ir)=rad(ir)*(two_pi*rad(ir)**2)*dvwf*bessel
894 : : end do
895 : 283 : call ctrap(mmax,work2,al,result)
896 : 283 : ypn=ztor1+result
897 : :
898 : : ! Fit spline for l=1 Kleinman-Bylander form factor
899 : 283 : call spline(qgrid,ffspl(1,1,2),mqgrid,yp1,ypn,ffspl(1,2,2))
900 : :
901 : : else
902 : : ! or else put form factor to 0 for l=1
903 [ + + + + ]: 232325 : ffspl(:,:,2)=0.0d0
904 : : end if
905 : : ! Endif condition of lmax>0
906 : : end if
907 : :
908 : : ! Finished if lmax=1 (highest nonlocal correction)
909 : : ! Do l=2 nonlocal form factor if eta(3) not 0 and lmax>=2
910 [ + + ]: 322 : if (lmax>1)then
911 [ + + ]: 284 : if(ekb(3)/=0.0d0) then
912 : :
913 : 43 : lp1=3
914 : : ! do q=0 separately; f_2(q=0) vanishes
915 : 43 : ffspl(1,1,3)=0.0d0
916 : :
917 : : ! do rest of q points
918 [ + + ]: 127909 : do iq=2,mqgrid
919 : 127866 : arg=two_pi*qgrid(iq)
920 : 127866 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
921 : 127866 : qr=arg*rad(1)
922 [ + + ]: 127866 : if(qr<1.d-3)then
923 : 118980 : bessel=qr*qr/15.0d0-qr**4/210.0d0
924 : : else
925 : 8886 : bessel=bes2_psp5(qr)
926 : : end if
927 : : ! ztor1=(bes2(arg*rad(1))*rad(1)*dvwf)*rad(1)/7.0d0
928 : 127866 : ztor1=(bessel*rad(1)*dvwf)*rad(1)/7.0d0
929 [ + + ]: 93120186 : do ir=1,mmax
930 : 92992320 : qr=arg*rad(ir)
931 [ + + ]: 92992320 : if(qr<1.d-3)then
932 : 22338050 : bessel=qr*qr/15.0d0-qr**4/210.0d0
933 : : else
934 : 70654270 : bessel=bes2_psp5(qr)
935 : : end if
936 : 92992320 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
937 : : ! work(ir)=rad(ir)*(r(ir)*bes2(arg*rad(ir))*dvwf)
938 : 93120186 : work3(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
939 : : end do
940 : 127866 : call ctrap(mmax,work3,al,result)
941 : 127909 : ffspl(iq,1,3)=ztor1+result
942 : : end do
943 : :
944 : : ! Compute yp1,ypn for l=2
945 : : ! yp1=0 for l=2
946 : 43 : yp1=0.0d0
947 : : ! ypn=$\int [2 \pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$
948 : : ! where x=2 Pi qgrid(mqgrid) r
949 : 43 : arg=two_pi*qgrid(mqgrid)
950 : 43 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
951 : 43 : qr=arg*rad(1)
952 [ + + ]: 43 : if(qr<1.d-3)then
953 : 37 : bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0
954 : : else
955 : 6 : bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr
956 : : end if
957 : : ! ztor1=( (two_pi*rad(1)**2)*dvwf* (bes1(arg*rad(1))-
958 : : ! 3.0d0*bes2(arg*rad(1))/(arg*rad(1))) ) * rad(1)/7.0d0
959 : 43 : ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/7.0d0
960 [ + + ]: 31594 : do ir=1,mmax
961 : 31551 : qr=arg*rad(ir)
962 [ + + ]: 31551 : if(qr<1.d-3)then
963 : 5753 : bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0
964 : : else
965 : 25798 : bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr
966 : : end if
967 : 31551 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
968 : : ! work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
969 : : ! (bes1(arg*rad(ir))-3.d0*bes2(arg*rad(ir))/(arg*rad(ir))) )
970 : 31594 : work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel)
971 : : end do
972 : 43 : call ctrap(mmax,work3,al,result)
973 : 43 : ypn=ztor1+result
974 : :
975 : : ! Fit spline for l=2 Kleinman-Bylander form factor
976 : 43 : call spline(qgrid,ffspl(1,1,3),mqgrid,yp1,ypn,ffspl(1,2,3))
977 : :
978 : : else
979 : : ! or else put form factor to 0 for l=1
980 [ + + + + ]: 1444205 : ffspl(:,:,3)=0.0d0
981 : : end if
982 : : ! Endif condition of lmax>1
983 : : end if
984 : :
985 : : ! Finished if lmax=2 (highest nonlocal correction)
986 : : ! Do l=3 nonlocal form factor if eta(4) not 0 and lmax>=3
987 [ + + ]: 322 : if (lmax>2)then
988 [ + - ]: 27 : if(ekb(4)/=0.0d0) then
989 : :
990 : 27 : lp1=4
991 : : ! do q=0 separately; f_3(q=0) vanishes
992 : 27 : ffspl(1,1,4)=0.0d0
993 : :
994 : : ! do rest of q points
995 [ + + ]: 79527 : do iq=2,mqgrid
996 : 79500 : arg=two_pi*qgrid(iq)
997 : 79500 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
998 : 79500 : qr=arg*rad(1)
999 [ + + ]: 79500 : if(qr<1.d-3)then
1000 : 74602 : bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0
1001 : : else
1002 : 4898 : bessel=bes3_psp5(qr)
1003 : : end if
1004 : : ! ztor1=(bes3(arg*rad(1))*rad(1)*dvwf)*rad(1)/9.0d0
1005 : 79500 : ztor1=(bessel*rad(1)*dvwf)*rad(1)/9.0d0
1006 [ + + ]: 48801000 : do ir=1,mmax
1007 : 48721500 : qr=arg*rad(ir)
1008 [ + + ]: 48721500 : if(qr<1.d-3)then
1009 : 7439225 : bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0
1010 : : else
1011 : 41282275 : bessel=bes3_psp5(qr)
1012 : : end if
1013 : 48721500 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1014 : : ! work(ir)=rad(ir)*(rad(ir)*bes3(arg*rad(ir))*dvwf)
1015 : 48801000 : work4(ir)=rad(ir)*(rad(ir)*bessel*dvwf)
1016 : : end do
1017 : 79500 : call ctrap(mmax,work4,al,result)
1018 : 79527 : ffspl(iq,1,4)=ztor1+result
1019 : : end do
1020 : :
1021 : : ! Compute yp1,ypn for l=3
1022 : : ! yp1=0 for l=3
1023 : 27 : yp1=0.0d0
1024 : : ! ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$
1025 : : ! where x=2 Pi qgrid(mqgrid) r
1026 : 27 : arg=two_pi*qgrid(mqgrid)
1027 : 27 : dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1)
1028 : 27 : qr=arg*rad(1)
1029 [ + + ]: 27 : if(qr<1.d-3)then
1030 : 24 : bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0
1031 : : else
1032 : 3 : bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr
1033 : : end if
1034 : : ! ztor1=( (two_pi*rad(1)**2)*dvwf* (bes2(arg*rad(1))-
1035 : : ! 3.0d0*bes3(arg*rad(1))/(arg*rad(1))) ) * rad(1)/9.0d0
1036 : 27 : ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/9.0d0
1037 [ + + ]: 16915 : do ir=1,mmax
1038 : 16888 : qr=arg*rad(ir)
1039 [ + + ]: 16888 : if(qr<1.d-3)then
1040 : 1548 : bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0
1041 : : else
1042 : 15340 : bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr
1043 : : end if
1044 : 16888 : dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)
1045 : : ! work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*
1046 : : ! (bes2(arg*rad(ir))-4.d0*bes3(arg*rad(ir))/(arg*rad(ir))) )
1047 : 16915 : work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel)
1048 : : end do
1049 : 27 : call ctrap(mmax,work4,al,result)
1050 : 27 : ypn=ztor1+result
1051 : :
1052 : : ! Fit spline for l=3 Kleinman-Bylander form factor
1053 : 27 : call spline(qgrid,ffspl(1,1,4),mqgrid,yp1,ypn,ffspl(1,2,4))
1054 : :
1055 : : else
1056 : : ! or else put form factor to 0 for l=3
1057 [ # # # # ]: 0 : ffspl(:,:,4)=0.0d0
1058 : : end if
1059 : : ! Endif condition of lmax>2
1060 : : end if
1061 : :
1062 : : ! Endif condition lmax/=-1
1063 : : end if
1064 : :
1065 : : !DEBUG
1066 : : !write(std_out,*) 'EKB=',(ekb(iq),iq=1,3)
1067 : : !write(std_out,*) 'COSKB=',(ckb(iq),iq=1,3)
1068 : : !ENDDEBUG
1069 : :
1070 [ - + ]: 322 : ABI_FREE(work1)
1071 [ - + ]: 322 : ABI_FREE(work2)
1072 [ - + ]: 322 : ABI_FREE(work3)
1073 [ - + ]: 322 : ABI_FREE(work4)
1074 : :
1075 : : contains
1076 : :
1077 : 382303960 : function bes0_psp5(arg)
1078 : : real(dp) :: bes0_psp5,arg
1079 : 382303960 : bes0_psp5=sin(arg)/arg
1080 : : end function bes0_psp5
1081 : :
1082 : 346700056 : function bes1_psp5(arg)
1083 : : real(dp) :: bes1_psp5,arg
1084 : 346700056 : bes1_psp5=(sin(arg)-arg*cos(arg))/arg**2
1085 : 346700056 : end function bes1_psp5
1086 : :
1087 : 70704303 : function bes2_psp5(arg)
1088 : : real(dp) :: bes2_psp5,arg
1089 : 70704303 : bes2_psp5=( (3.0d0-arg**2)*sin(arg)- 3.0d0*arg*cos(arg))/arg**3
1090 : 70704303 : end function bes2_psp5
1091 : :
1092 : 41302516 : function bes3_psp5(arg)
1093 : : real(dp) :: bes3_psp5, arg
1094 : : bes3_psp5=(15.d0*sin(arg)-15.d0*arg*cos(arg) &
1095 : 41302516 : -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4
1096 : 41302516 : end function bes3_psp5
1097 : :
1098 : : end subroutine psp5nl
1099 : : !!***
1100 : :
1101 : : !!****f* ABINIT/psp8lo
1102 : : !! NAME
1103 : : !! psp8lo
1104 : : !!
1105 : : !! FUNCTION
1106 : : !! Compute sine transform to transform from V(r) to q^2 V(q).
1107 : : !! Computes integrals on linear grid interpolated from the linear input
1108 : : !! grid with a spacing adjusted to ensure convergence at the maximum
1109 : : !! wavevector using corrected trapezoidal integration.
1110 : : !!
1111 : : !! INPUTS
1112 : : !! amesh=spacing for linear radial atomic grid.
1113 : : !! mmax=number of radial r grid points
1114 : : !! mqgrid=number of grid points in q from 0 to qmax.
1115 : : !! qgrid(mqgrid)=q grid values (bohr**-1).
1116 : : !! rad(mmax)=r grid values (bohr).
1117 : : !! vloc(mmax)=V(r) on radial grid.
1118 : : !! zion=nominal valence charge of atom.
1119 : : !!
1120 : : !! OUTPUT
1121 : : !! epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
1122 : : !!{{\\ \begin{equation}
1123 : : !! q2vq(mqgrid)
1124 : : !! =q^2 V(q)
1125 : : !! = -\frac{Zv}{\pi}
1126 : : !! + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
1127 : : !!\end{equation} }}
1128 : : !! yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax
1129 : : !! (needed for spline fitter).
1130 : : !!
1131 : : !! SOURCE
1132 : :
1133 : 510 : subroutine psp8lo(amesh, epsatm, mmax, mqgrid, qgrid, q2vq, rad, vloc, yp1, ypn, zion)
1134 : :
1135 : : !Arguments----------------------------------------------------------
1136 : : !scalars
1137 : : integer,intent(in) :: mmax,mqgrid
1138 : : real(dp),intent(in) :: amesh,zion
1139 : : real(dp),intent(out) :: epsatm,yp1,ypn
1140 : : !arrays
1141 : : real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax)
1142 : : real(dp),intent(out) :: q2vq(mqgrid)
1143 : :
1144 : : !Local variables-------------------------------
1145 : : !Following parameter controls accuracy of Fourier transform based on qmax
1146 : : !and represents the minimun number of integration points in one period scalars
1147 : : integer,parameter :: NPT_IN_2PI=200
1148 : : integer :: ider,iq,ir,irmu,irn,mesh_mult,mmax_new
1149 : : real(dp) :: amesh_new,arg,fp1,fpn,qmesh,result,ztor1
1150 : : !arrays
1151 : 510 : real(dp),allocatable :: rad_new(:),rvlpz(:),rvlpz_new(:),sprvlpz(:,:),work(:)
1152 : : ! *************************************************************************
1153 : :
1154 [ - + + - : 1530 : ABI_MALLOC(work,(mmax))
+ - + - -
+ - + ]
1155 [ - + + - : 1020 : ABI_MALLOC(rvlpz,(mmax))
- + - + ]
1156 : :
1157 : : !Do q=0 separately (compute epsatm)
1158 : 510 : ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2
1159 : : !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$
1160 [ + + ]: 293123 : do ir=1,mmax
1161 : 292613 : rvlpz(ir)=rad(ir)*vloc(ir)+zion
1162 : 293123 : work(ir)=rad(ir)*rvlpz(ir)
1163 : : end do
1164 : :
1165 : : !Do integral from zero to r(max)
1166 : 510 : call ctrap(mmax,work,amesh,result)
1167 : :
1168 : 510 : epsatm=4.d0*pi*result
1169 : 510 : q2vq(1)=-zion/pi
1170 : :
1171 : : !Find r mesh spacing necessary for accurate integration at qmax
1172 : 510 : amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1173 : :
1174 : : !Choose submultiple of input mesh
1175 : 510 : mesh_mult=int(amesh/amesh_new) + 1
1176 : : !mesh_mult = 1 ! DEBUG
1177 : 510 : mmax_new=mesh_mult*(mmax-1)+1
1178 : 510 : amesh_new=amesh/dble(mesh_mult)
1179 : :
1180 [ - + + - : 1530 : ABI_MALLOC(rad_new,(mmax_new))
+ - + - -
+ - + ]
1181 [ - + + - : 1020 : ABI_MALLOC(rvlpz_new,(mmax_new))
- + - + ]
1182 : :
1183 : : !print *, "in psp8lo with mesh_mult:", mesh_mult
1184 : : !print *, "in psp8lo with mmax_new:", mmax_new
1185 : : !print *, "in psp8lo with amesh_new:", amesh_new
1186 : :
1187 [ + + ]: 510 : if(mesh_mult==1) then
1188 [ + + ]: 248629 : rad_new(:)=rad(:)
1189 [ + + ]: 248629 : rvlpz_new(:)=rvlpz(:)
1190 : : else
1191 : : ! Set up spline and interpolate to finer mesh.
1192 : : ! First, compute derivatives at end points
1193 : : fp1=(-50.d0*rvlpz(1)+96.d0*rvlpz(2)-72.d0*rvlpz(3)+32.d0*rvlpz(4)&
1194 : 78 : & -6.d0*rvlpz(5))/(24.d0*amesh)
1195 : : fpn=(6.d0*rvlpz(mmax-4)-32.d0*rvlpz(mmax-3)+72.d0*rvlpz(mmax-2)&
1196 : 78 : & -96.d0*rvlpz(mmax-1)+50.d0*rvlpz(mmax))/(24.d0*amesh)
1197 [ - + + - : 234 : ABI_MALLOC(sprvlpz,(mmax,2))
+ - + - -
+ - + ]
1198 [ + + ]: 44494 : work(:)=zero
1199 : :
1200 : : ! Spline fit
1201 : 78 : call spline(rad, rvlpz,mmax,fp1,fpn,sprvlpz(:,2))
1202 [ + + ]: 44494 : sprvlpz(:,1)=rvlpz(:)
1203 : :
1204 : : ! Set up new radial mesh
1205 : : irn=1
1206 [ + + ]: 44416 : do ir=1,mmax-1
1207 [ + + ]: 133404 : do irmu=0,mesh_mult-1
1208 : 88988 : rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
1209 : 133326 : irn=irn+1
1210 : : end do
1211 : : end do
1212 : 78 : rad_new(mmax_new)=rad(mmax)
1213 : :
1214 : 78 : ider=0
1215 : 78 : call splfit(rad,work,sprvlpz,ider,rad_new,rvlpz_new,mmax,mmax_new)
1216 : :
1217 [ - + ]: 78 : ABI_FREE(sprvlpz)
1218 [ - + ]: 78 : ABI_FREE(work)
1219 [ - + + - : 156 : ABI_MALLOC(work,(mmax_new))
- + - + ]
1220 : : end if
1221 : :
1222 : : !Loop over q values
1223 [ + + ]: 1548893 : do iq=2,mqgrid
1224 : 1548383 : arg=2.d0*pi*qgrid(iq)
1225 : :
1226 : : ! Set up integrand
1227 [ + + ]: 1031190390 : do ir=1,mmax_new
1228 : 1031190390 : work(ir)=sin(arg*rad_new(ir))*rvlpz_new(ir)
1229 : : end do
1230 : :
1231 : : ! Do integral from zero to rad(mmax)
1232 : 1548383 : call ctrap(mmax_new,work,amesh_new,result)
1233 : :
1234 : : ! Store q^2 v(q)
1235 : 1548893 : q2vq(iq)=q2vq(1)+2.d0*qgrid(iq)*result
1236 : :
1237 : : end do
1238 : :
1239 : : !Compute derivatives of q^2 v(q) at ends of interval
1240 : 510 : qmesh=qgrid(2)-qgrid(1)
1241 : : yp1=(-50.d0*q2vq(1)+96.d0*q2vq(2)-72.d0*q2vq(3)+32.d0*q2vq(4)&
1242 : 510 : & -6.d0*q2vq(5))/(24.d0*qmesh)
1243 : : ypn=(6.d0*q2vq(mqgrid-4)-32.d0*q2vq(mqgrid-3)+72.d0*q2vq(mqgrid-2)&
1244 : 510 : & -96.d0*q2vq(mqgrid-1)+50.d0*q2vq(mqgrid))/(24.d0*qmesh)
1245 : :
1246 [ - + ]: 510 : ABI_FREE(work)
1247 [ - + ]: 510 : ABI_FREE(rad_new)
1248 [ - + ]: 510 : ABI_FREE(rvlpz_new)
1249 [ - + ]: 510 : ABI_FREE(rvlpz)
1250 : :
1251 : 510 : end subroutine psp8lo
1252 : : !!***
1253 : :
1254 : : !!****f* ABINIT/psp8nl
1255 : : !! NAME
1256 : : !! psp8nl
1257 : : !!
1258 : : !! FUNCTION
1259 : : !! Make Kleinman-Bylander/Bloechl form factors f_ln(q) for each
1260 : : !! projector n for each angular momentum l excepting an l corresponding to the local potential.
1261 : : !! Note that an arbitrary local potential can be used, so all l from 0 to lmax may be represented.
1262 : : !!
1263 : : !! INPUTS
1264 : : !! amesh=grid spacing for uniform (linear) radial grid
1265 : : !! indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
1266 : : !! or i=lmn (if useylm=1)
1267 : : !! lmax=maximum ang momentum for which nonlocal form factor is desired. lmax <= 2 allowed.
1268 : : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
1269 : : !! =if useylm=0, max number of (l,n) comp. over all type of psps
1270 : : !! lnmax=max. number of (l,n) components over all type of psps
1271 : : !! mmax=number of radial grid points for atomic grid
1272 : : !! mqgrid=number of grid points for q grid
1273 : : !! qgrid(mqgrid)=values at which form factors are returned
1274 : : !! rad(mmax)=radial grid values
1275 : : !! vpspll(mmax,lnmax)=nonlocal projectors for each (l,n) on linear radial grid.
1276 : : !! Here, these are the product of the reference
1277 : : !! wave functions and (v(l,n)-vloc), calculated in the psp generation
1278 : : !! program and normalized so that integral(0,rc(l)) vpsll^2 dr = 1,
1279 : : !! which leads to the the usual convention for the energies ekb(l,n)
1280 : : !! also calculated in the psp generation program.
1281 : : !!
1282 : : !! OUTPUT
1283 : : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_ln(q) and
1284 : : !! second derivative from spline fit for each (l,n).
1285 : : !!
1286 : : !! NOTES
1287 : : !! u_l(r) is reference state wavefunction (input as wf);
1288 : : !! j_l(q) is a spherical Bessel function;
1289 : : !! dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
1290 : : !! f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
1291 : : !! where dvms = $\int_0^{rmax} [(u_l(r) dV_l(r))^2 dr]$ is the mean
1292 : : !! square value of the nonlocal correction for angular momentum l.
1293 : : !! Xavier Gonze s E_KB = $ dvms/\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]$.
1294 : : !! This is the eigenvalue of the Kleinman-Bylander operator and sets
1295 : : !! the energy scale of the nonlocal psp corrections.
1296 : : !!
1297 : : !! SOURCE
1298 : :
1299 : 510 : subroutine psp8nl(amesh, ffspl, indlmn, lmax, lmnmax, lnmax, mmax, mqgrid, qgrid, rad, vpspll)
1300 : :
1301 : : !Arguments----------------------------------------------------------
1302 : : !scalars
1303 : : integer,intent(in) :: lmax,lmnmax,lnmax,mmax,mqgrid
1304 : : real(dp),intent(in) :: amesh
1305 : : !arrays
1306 : : integer,intent(in) :: indlmn(6,lmnmax)
1307 : : real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vpspll(mmax,lnmax)
1308 : : real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
1309 : :
1310 : : !Local variables-------------------------------
1311 : : !Following parameter controls accuracy of Fourier transform based on qmax
1312 : : !and represents the minimun number of integration points in one period.
1313 : : !scalars
1314 : : integer,parameter :: NPT_IN_2PI=200
1315 : : integer :: iln,iln0,ilmn,iq,ir,irmu,irn,ll,mesh_mult,mmax_new,mvpspll
1316 : : real(dp) :: amesh_new,arg,c1,c2,c3,c4,dri,qmesh,result,tv,xp,xpm1,xpm2,xpp1,yp1,ypn
1317 : : !arrays
1318 : : real(dp) :: sb_out(4)
1319 : 510 : real(dp),allocatable :: rad_new(:),vpspll_new(:,:),work(:,:),work2(:)
1320 : : ! *************************************************************************
1321 : :
1322 : : ! Find r mesh spacing necessary for accurate integration at qmax
1323 : 510 : amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1324 : :
1325 : : ! Choose submultiple of input mesh
1326 : 510 : mesh_mult=int(amesh/amesh_new) + 1
1327 : 510 : mmax_new=mesh_mult*(mmax-1)+1
1328 : 510 : amesh_new=amesh/dble(mesh_mult)
1329 : :
1330 [ - + + - : 1530 : ABI_MALLOC(rad_new,(mmax_new))
+ - + - -
+ - + ]
1331 [ - + + - : 2040 : ABI_MALLOC(vpspll_new,(mmax_new,lnmax))
- + + - +
- + - - +
- + ]
1332 : :
1333 [ + + ]: 510 : if (mesh_mult == 1) then
1334 [ + + ]: 253682 : rad_new(:)=rad(:)
1335 : : else
1336 : : ! Set up new radial mesh
1337 : : irn=1
1338 [ + + ]: 45331 : do ir=1,mmax-1
1339 [ + + ]: 136149 : do irmu=0,mesh_mult-1
1340 : 90818 : rad_new(irn)=rad(ir)+dble(irmu)*amesh_new
1341 : 136071 : irn=irn+1
1342 : : end do
1343 : : end do
1344 : 78 : rad_new(mmax_new)=rad(mmax)
1345 : : end if
1346 : :
1347 : : ! Interpolate projectors onto new grid if called for
1348 : : ! Cubic polynomial interpolation is used which is consistent
1349 : : ! with the original interpolation of these functions from
1350 : : ! a log grid to the input linear grid.
1351 : 510 : dri = one/amesh
1352 [ + + ]: 344656 : do irn=1,mmax_new
1353 : : ! index to find bracketing input mesh points
1354 [ + + ]: 344656 : if(mesh_mult>1) then
1355 : 90896 : ir = irn/mesh_mult + 1
1356 : 90896 : ir = max(ir,2)
1357 : 90896 : ir = min(ir,mmax-2)
1358 : : ! interpolation coefficients
1359 : 90896 : xp = dri * (rad_new(irn) - rad(ir))
1360 : 90896 : xpp1 = xp + one
1361 : 90896 : xpm1 = xp - one
1362 : 90896 : xpm2 = xp - two
1363 : 90896 : c1 = -xp * xpm1 * xpm2 * sixth
1364 : 90896 : c2 = xpp1 * xpm1 * xpm2 * half
1365 : 90896 : c3 = - xp * xpp1 * xpm2 * half
1366 : 90896 : c4 = xp * xpp1 * xpm1 * sixth
1367 : :
1368 : : ! Now do the interpolation on all projectors for this grid point
1369 : 90896 : iln0=0
1370 [ + + ]: 501103 : do ilmn=1,lmnmax
1371 : 410207 : iln=indlmn(5,ilmn)
1372 [ + + ]: 501103 : if (iln>iln0) then
1373 : 387032 : iln0=iln
1374 : : tv = c1 * vpspll(ir - 1, iln) &
1375 : : & + c2 * vpspll(ir , iln) &
1376 : : & + c3 * vpspll(ir + 1, iln) &
1377 : 387032 : & + c4 * vpspll(ir + 2, iln)
1378 [ + + ]: 387032 : if(abs(tv)>tol10) then
1379 : 93636 : vpspll_new(irn,iln)=tv
1380 : 93636 : mvpspll=irn
1381 : : else
1382 : 293396 : vpspll_new(irn,iln)=zero
1383 : : end if
1384 : : end if
1385 : : end do
1386 : :
1387 : : else
1388 : : ! With no mesh multiplication, just copy projectors
1389 : 3085814 : ir=irn
1390 : : iln0=0
1391 [ + + ]: 3085814 : do ilmn=1,lmnmax
1392 : 2832564 : iln=indlmn(5,ilmn)
1393 [ + + ]: 3085814 : if (iln>iln0) then
1394 : 1921400 : iln0=iln
1395 : 1921400 : tv = vpspll(ir,iln)
1396 [ + + ]: 1921400 : if(abs(tv)>tol10) then
1397 : 558923 : vpspll_new(irn,iln)=tv
1398 : 558923 : mvpspll=irn
1399 : : else
1400 : 1362477 : vpspll_new(irn,iln)=zero
1401 : : end if
1402 : : end if
1403 : : end do
1404 : :
1405 : : end if
1406 : : end do ! irn
1407 : :
1408 [ - + + - : 2040 : ABI_MALLOC(work, (mvpspll,lnmax))
- + + - +
- + - - +
+ - ]
1409 : :
1410 : : ! Loop over q values
1411 [ + + ]: 1549403 : do iq=1,mqgrid
1412 : 1548893 : arg=2.d0*pi*qgrid(iq)
1413 : :
1414 : : ! Set up integrands
1415 [ + + ]: 314288135 : do ir=1,mvpspll
1416 : 312739242 : call sbf8(lmax+1,arg*rad_new(ir),sb_out)
1417 : 312739242 : iln0=0
1418 [ + + ]: 3304690469 : do ilmn=1,lmnmax
1419 : 2990402334 : iln=indlmn(5,ilmn)
1420 [ + + ]: 3303141576 : if (iln>iln0) then
1421 : 2159656511 : iln0=iln
1422 : 2159656511 : ll=indlmn(1,ilmn)
1423 : 2159656511 : work(ir,iln)=sb_out(ll+1)*vpspll_new(ir,iln)*rad_new(ir)
1424 : : end if
1425 : : end do
1426 : : end do !ir
1427 : :
1428 : : ! Do integral from zero to rad_new(mvpspll)
1429 : : iln0=0
1430 [ + + ]: 16733936 : do ilmn=1,lmnmax
1431 : 15184533 : iln=indlmn(5,ilmn)
1432 [ + + ]: 16733426 : if (iln>iln0) then
1433 : 10632016 : iln0=iln
1434 : 10632016 : call ctrap(mvpspll,work(1,iln),amesh_new,result)
1435 : 10632016 : ffspl(iq,1,iln)=result
1436 : : end if
1437 : : end do
1438 : : end do ! iq mesh
1439 : :
1440 : : ! Fit splines for form factors
1441 [ - + + - : 1530 : ABI_MALLOC(work2,(mqgrid))
+ - + - -
+ - + ]
1442 : 510 : qmesh=qgrid(2)-qgrid(1)
1443 : :
1444 : 510 : iln0=0
1445 [ + + ]: 5545 : do ilmn=1,lmnmax
1446 : 5035 : iln=indlmn(5,ilmn)
1447 [ + + ]: 5545 : if (iln>iln0) then
1448 : 3518 : iln0=iln
1449 : : ! Compute derivatives of form factors at ends of interval
1450 : : yp1=(-50.d0*ffspl(1,1,iln)+96.d0*ffspl(2,1,iln)-72.d0*ffspl(3,1,iln)&
1451 : 3518 : & +32.d0*ffspl(4,1,iln)- 6.d0*ffspl(5,1,iln))/(24.d0*qmesh)
1452 : : ypn=(6.d0*ffspl(mqgrid-4,1,iln)-32.d0*ffspl(mqgrid-3,1,iln)&
1453 : : & +72.d0*ffspl(mqgrid-2,1,iln)-96.d0*ffspl(mqgrid-1,1,iln)&
1454 : 3518 : & +50.d0*ffspl(mqgrid,1,iln))/(24.d0*qmesh)
1455 : :
1456 : 3518 : call spline(qgrid,ffspl(1,1,iln),mqgrid,yp1,ypn,ffspl(1,2,iln))
1457 : : end if
1458 : : end do
1459 : :
1460 [ - + ]: 510 : ABI_FREE(rad_new)
1461 [ - + ]: 510 : ABI_FREE(vpspll_new)
1462 [ - + ]: 510 : ABI_FREE(work)
1463 [ - + ]: 510 : ABI_FREE(work2)
1464 : :
1465 : 510 : end subroutine psp8nl
1466 : : !!***
1467 : :
1468 : : !!****f* ABINIT/cc_derivatives
1469 : : !! NAME
1470 : : !! cc_derivatives
1471 : : !!
1472 : : !! FUNCTION
1473 : : !! subroutine to spline the core charge and get derivatives
1474 : : !! extracted from previous version of psp6cc_drh
1475 : : !! input on log grid, and splined to regular grid between 0 and rchrg
1476 : : !!
1477 : : !! INPUTS
1478 : : !! mmax=maximum number of points in real space grid in the psp file
1479 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
1480 : : !! rchrg=cut-off radius for the core density
1481 : : !! rad=radial grid points
1482 : : !! ff=core charge at points in rad
1483 : : !! ff1=first derivative of ff on log grid
1484 : : !! ff2=second derivative of ff on log grid
1485 : : !!
1486 : : !! OUTPUT
1487 : : !! xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
1488 : : !!
1489 : : !! NOTES
1490 : : !! Test version by DRH - requires very smooth model core charge
1491 : : !!
1492 : : !! SOURCE
1493 : :
1494 : 22 : subroutine cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d)
1495 : :
1496 : : !Arguments ------------------------------------
1497 : : ! scalars
1498 : : integer,intent(in) :: mmax,n1xccc
1499 : : real(dp),intent(in) :: rchrg
1500 : : !arrays
1501 : : real(dp),intent(in) :: rad(mmax),ff(mmax),ff1(mmax),ff2(mmax)
1502 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
1503 : :
1504 : : !Local variables-------------------------------
1505 : : ! scalars
1506 : : integer :: i1xccc
1507 : : real(dp) :: der1,dern
1508 : : !arrays
1509 : 22 : real(dp),allocatable :: ff3(:),ff4(:),gg(:),gg1(:),gg2(:)
1510 : 22 : real(dp),allocatable :: gg3(:),gg4(:),work(:),xx(:)
1511 : : ! *************************************************************************
1512 : :
1513 : : !write(std_out,*) 'cc_derivatives : enter'
1514 : :
1515 [ - + + - : 66 : ABI_MALLOC(ff3, (mmax))
+ - + - -
+ - + ]
1516 [ - + + - : 44 : ABI_MALLOC(ff4, (mmax))
- + - + ]
1517 [ - + + - : 66 : ABI_MALLOC(gg, (n1xccc))
+ - + - -
+ - + ]
1518 [ - + + - : 44 : ABI_MALLOC(gg1, (n1xccc))
- + - + ]
1519 [ - + + - : 44 : ABI_MALLOC(gg2, (n1xccc))
- + - + ]
1520 [ - + + - : 44 : ABI_MALLOC(gg3, (n1xccc))
- + - + ]
1521 [ - + + - : 44 : ABI_MALLOC(gg4, (n1xccc))
- + - + ]
1522 [ - + + - : 44 : ABI_MALLOC(work, (mmax))
- + - + ]
1523 [ - + + - : 44 : ABI_MALLOC(xx, (n1xccc))
- + - + ]
1524 : :
1525 : : ! calculate third derivative ff3 on logarithmic grid
1526 : 22 : der1=ff2(1)
1527 : 22 : dern=ff2(mmax)
1528 : 22 : call spline(rad,ff1,mmax,der1,dern,ff3)
1529 : :
1530 : : ! calculate fourth derivative ff4 on logarithmic grid
1531 : 22 : der1=0.d0
1532 : 22 : dern=0.d0
1533 : 22 : call spline(rad,ff2,mmax,der1,dern,ff4)
1534 : :
1535 : : ! generate uniform mesh xx in the box cut by rchrg:
1536 [ + + ]: 55044 : do i1xccc=1,n1xccc
1537 : 55044 : xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1)
1538 : : end do
1539 : :
1540 : : !now interpolate core charge and derivatives on the uniform grid
1541 : : !core charge, input=ff, output=gg
1542 : 22 : call splint(mmax,rad,ff,ff2,n1xccc,xx,gg)
1543 : :
1544 : : ! first derivative input=ff1, output=gg1
1545 : 22 : call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1)
1546 : :
1547 : : !normalize gg1
1548 : : !gg1(:)=gg1(:)*rchrg
1549 : :
1550 : : ! second derivative input=ff2, output=gg2
1551 : 22 : call splint(mmax,rad,ff2,ff4,n1xccc,xx,gg2)
1552 : :
1553 : : !normalize gg2
1554 : : !gg2(:)=gg2(:)*rchrg**2
1555 : :
1556 : : ! reallocate work otherwise the calls to spline crash (n1xccc /= mmax)
1557 [ - + ]: 22 : ABI_FREE(work)
1558 [ - + + - : 44 : ABI_MALLOC(work, (n1xccc))
- + - + ]
1559 : :
1560 : : !recalculate 3rd derivative consistent with spline fit to first derivative on linear grid
1561 : 22 : der1=gg2(1)
1562 : 22 : dern=gg2(n1xccc)
1563 : 22 : call spline(xx,gg1,n1xccc,der1,dern,gg3)
1564 : :
1565 : : !calculate 4th derivative consistent with spline fit to second derivative on linear grid
1566 : 22 : der1=0.0d0
1567 : 22 : dern=0.0d0
1568 : 22 : call spline(xx,gg2,n1xccc,der1,dern,gg4)
1569 : :
1570 : : !now calculate second to fourth derivative by forward differences
1571 : : !to avoid numerical noise uses a smoothing function
1572 : : !
1573 : : !call smooth(gg1,n1xccc,10)
1574 : :
1575 : : !gg2(n1xccc)=0.0
1576 : : !do i1xccc=1,n1xccc-1
1577 : : !gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1)
1578 : : !end do
1579 : :
1580 : : !call smooth(gg2,n1xccc,10)
1581 : :
1582 : : !gg3(n1xccc)=0.0
1583 : : !do i1xccc=1,n1xccc-1
1584 : : !gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1)
1585 : : !end do
1586 : :
1587 : : !call smooth(gg3,n1xccc,10)
1588 : :
1589 : : !gg4(n1xccc)=0.0
1590 : : !do i1xccc=1,n1xccc-1
1591 : : !gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1)
1592 : : !end do
1593 : :
1594 : : !call smooth(gg4,n1xccc,10)
1595 : :
1596 : : !write on xcc1d
1597 : : !normalize to unit range usage later in program
1598 [ + + ]: 55044 : xccc1d(:,1)=gg(:)
1599 [ + + ]: 55044 : xccc1d(:,2)=gg1(:)*rchrg
1600 [ + + ]: 55044 : xccc1d(:,3)=gg2(:)*rchrg**2
1601 [ + + ]: 55044 : xccc1d(:,4)=gg3(:)*rchrg**3
1602 [ + + ]: 55044 : xccc1d(:,5)=gg4(:)*rchrg**4
1603 : : !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc
1604 : :
1605 : : !DEBUG
1606 : : !note: the normalization condition is the following:
1607 : : !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg
1608 : : !
1609 : : !norm=0.d0
1610 : : !do i1xccc=1,n1xccc
1611 : : !norm = norm + 4.d0*pi*rchrg/dble(n1xccc-1)*&
1612 : : !& xx(i1xccc)**2*xccc1d(i1xccc,1)
1613 : : !end do
1614 : : !write(std_out,*) ' norm=',norm
1615 : : !
1616 : : !write(std_out,*)' psp6cc_drh : output of core charge density and derivatives '
1617 : : !write(std_out,*)' xx gg gg1 '
1618 : : !do i1xccc=1,n1xccc
1619 : : !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
1620 : : !end do
1621 : : !write(std_out,*)' xx gg2 gg3 '
1622 : : !do i1xccc=1,n1xccc
1623 : : !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
1624 : : !end do
1625 : : !write(std_out,*)' xx gg4 gg5 '
1626 : : !do i1xccc=1,n1xccc
1627 : : !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
1628 : : !end do
1629 : : !write(std_out,*)' psp1cc : debug done, stop '
1630 : : !stop
1631 : : !ENDDEBUG
1632 : :
1633 [ - + ]: 22 : ABI_FREE(ff3)
1634 [ - + ]: 22 : ABI_FREE(ff4)
1635 [ - + ]: 22 : ABI_FREE(gg)
1636 [ - + ]: 22 : ABI_FREE(gg1)
1637 : 22 : ABI_FREE(gg2)
1638 [ - + ]: 22 : ABI_FREE(gg3)
1639 [ - + ]: 22 : ABI_FREE(gg4)
1640 [ - + ]: 22 : ABI_FREE(work)
1641 [ - + ]: 22 : ABI_FREE(xx)
1642 : :
1643 : 22 : end subroutine cc_derivatives
1644 : : !!***
1645 : :
1646 : : end module m_psptk
1647 : : !!***
|
