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 587 : 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 587 : real(dp),allocatable :: ff(:),ff2(:),work(:),yy(:)
104 : ! *************************************************************************
105 :
106 1761 : ABI_MALLOC(ff,(n1xccc))
107 1174 : ABI_MALLOC(ff2,(n1xccc))
108 1174 : ABI_MALLOC(work,(n1xccc))
109 1174 : ABI_MALLOC(yy,(n1xccc))
110 :
111 587 : if(n1xccc > 1)then
112 587 : factor=one/dble(n1xccc-1)
113 1468674 : do i1xccc=1,n1xccc
114 1468674 : 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 7631 : xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero
123 :
124 : !Take care of each derivative separately
125 2348 : do ider=0,2
126 :
127 1761 : if(ider==0)then
128 : ! Generate spline fitting for the function gg
129 1468674 : do i1xccc=1,n1xccc
130 1468087 : xx=three*yy(i1xccc)
131 1468087 : call gg1cc(gg1cc_xx,xx)
132 1468674 : ff(i1xccc)=fchrg*gg1cc_xx
133 : end do
134 : ! Complete with derivatives at end points
135 587 : der1=zero
136 587 : call gp1cc(gp1cc_xx,three)
137 587 : dern=three*fchrg*gp1cc_xx
138 1174 : else if(ider==1)then
139 : ! Generate spline fitting for the function gp
140 1468674 : do i1xccc=1,n1xccc
141 1468087 : xx=three*yy(i1xccc)
142 1468087 : call gp1cc(gp1cc_xx,xx)
143 1468674 : ff(i1xccc)=three*fchrg*gp1cc_xx
144 : end do
145 : ! Complete with derivatives at end points, already estimated
146 587 : der1=xccc1d(1,ider+2)
147 587 : 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 1468674 : do i1xccc=1,n1xccc
154 1468087 : xx=three*yy(i1xccc)
155 1468087 : call gpp1cc(gpp1cc_xx,xx)
156 1468674 : ff(i1xccc)=9.0_dp*fchrg*gpp1cc_xx
157 : end do
158 : ! Complete with derivatives of end points
159 587 : der1=xccc1d(1,ider+2)
160 587 : dern=xccc1d(n1xccc,ider+2)
161 : end if
162 :
163 : ! Produce second derivative numerically, for use with splines
164 1761 : call spline(yy,ff,n1xccc,der1,dern,ff2)
165 4406022 : xccc1d(:,ider+1)=ff(:)
166 4406609 : xccc1d(:,ider+3)=ff2(:)
167 : end do
168 :
169 1468674 : 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 587 : ABI_FREE(ff)
190 587 : ABI_FREE(ff2)
191 587 : ABI_FREE(work)
192 587 : ABI_FREE(yy)
193 :
194 587 : 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 1468087 : 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 1468087 : 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 1468087 : else if (abs(xx)<=1.d-09) then
233 587 : gg1cc_xx=1.d0
234 1467500 : 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 1467500 : 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 1467500 : & (1.d0-4.0d0*xx**2)*(1.d0-xx**2) ) )**2
245 : end if
246 :
247 1468087 : 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 1468674 : 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 1468674 : 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 978529 : denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
299 978529 : phi=denom*sin(two_pi*xx)*two_pim1
300 978529 : phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
301 978529 : gp1cc_xx=2.d0*phi*phip
302 : ! Handle limits where denominator vanishes
303 490145 : else if (abs(xx)<1.d-03) then
304 587 : gp1cc_xx=xx*(c11+xx**2*c12)
305 489558 : else if (abs(xx-0.5d0)<=1.d-03) then
306 1174 : gp1cc_xx=c21+(xx-0.5d0)*(c22+(xx-0.5d0)*(c23+(xx-0.5d0)*(c24+(xx-0.5d0)*c25)))
307 488384 : else if (abs(xx-1.d0)<=1.d-03) then
308 1174 : 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 487210 : denom=1.d0/(xx*(1.d0-4.d0*xx**2)*(1.d0-xx**2))
312 487210 : phi=denom*sin(two_pi*xx)*two_pim1
313 487210 : phip=denom*(cos(two_pi*xx)-(1.d0-xx**2*(15.d0-xx**2*20))*phi)
314 487210 : gp1cc_xx=2.d0*phi*phip
315 : end if
316 :
317 1468674 : 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 1468087 : 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 1468087 : 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 1468087 : else if (abs(xx)<=1.d-09) then
364 587 : gpp1cc_xx=c1
365 1467500 : 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 1467500 : 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 1467500 : t1 = xx**2
374 1467500 : t2 = 1/t1
375 1467500 : t3 = 1/Pi
376 1467500 : t4 = 2*xx
377 1467500 : t5 = t4-1
378 1467500 : t6 = t5**2
379 1467500 : t7 = 1/t6
380 1467500 : t8 = t4+1
381 1467500 : t9 = t8**2
382 1467500 : t10 = 1/t9
383 1467500 : t11 = xx-1
384 1467500 : t12 = t11**2
385 1467500 : t14 = 1/t12/t11
386 1467500 : t15 = xx+1
387 1467500 : t16 = t15**2
388 1467500 : t17 = 1/t16
389 1467500 : t18 = Pi*xx
390 1467500 : t19 = sin(t18)
391 1467500 : t20 = cos(t18)
392 1467500 : t21 = t20**2
393 1467500 : t23 = t19*t21*t20
394 1467500 : t24 = t17*t23
395 1467500 : t31 = t19**2
396 1467500 : t33 = t31*t19*t20
397 1467500 : t34 = t17*t33
398 1467500 : t41 = Pi**2
399 1467500 : t42 = 1/t41
400 1467500 : t44 = 1/t16/t15
401 1467500 : t45 = t31*t21
402 1467500 : t46 = t44*t45
403 1467500 : t54 = 1/t1/xx
404 1467500 : t55 = 1/t12
405 1467500 : t56 = t55*t46
406 1467500 : t57 = t10*t56
407 1467500 : t62 = t9**2
408 1467500 : t64 = t17*t45
409 1467500 : t65 = t55*t64
410 1467500 : t72 = 1/t9/t8
411 1467500 : t78 = t14*t64
412 1467500 : t79 = t10*t78
413 1467500 : t85 = t12**2
414 1467500 : t93 = t21**2
415 1467500 : t100 = t31**2
416 1467500 : t120 = 1/t6/t5
417 1467500 : t121 = t55*t34
418 1467500 : t122 = t10*t121
419 1467500 : t127 = t16**2
420 1467500 : t138 = t6**2
421 1467500 : t140 = t10*t65
422 1467500 : t152 = t72*t65
423 1467500 : t157 = t7*t140
424 1467500 : t160 = t1**2
425 1467500 : t174 = t55*t24
426 1467500 : 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 1467500 : & 6/t160*t42*t157
439 : end if
440 :
441 1468087 : 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 504 : 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 504 : real(dp),allocatable :: rad_new(:),rvlpz(:),rvlpz_new(:),sprvlpz(:,:),work(:)
1152 : ! *************************************************************************
1153 :
1154 1512 : ABI_MALLOC(work,(mmax))
1155 1008 : ABI_MALLOC(rvlpz,(mmax))
1156 :
1157 : !Do q=0 separately (compute epsatm)
1158 504 : 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 289517 : do ir=1,mmax
1161 289013 : rvlpz(ir)=rad(ir)*vloc(ir)+zion
1162 289517 : work(ir)=rad(ir)*rvlpz(ir)
1163 : end do
1164 :
1165 : !Do integral from zero to r(max)
1166 504 : call ctrap(mmax,work,amesh,result)
1167 :
1168 504 : epsatm=4.d0*pi*result
1169 504 : q2vq(1)=-zion/pi
1170 :
1171 : !Find r mesh spacing necessary for accurate integration at qmax
1172 504 : amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1173 :
1174 : !Choose submultiple of input mesh
1175 504 : mesh_mult=int(amesh/amesh_new) + 1
1176 : !mesh_mult = 1 ! DEBUG
1177 504 : mmax_new=mesh_mult*(mmax-1)+1
1178 504 : amesh_new=amesh/dble(mesh_mult)
1179 :
1180 1512 : ABI_MALLOC(rad_new,(mmax_new))
1181 1008 : 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 504 : if(mesh_mult==1) then
1188 245023 : rad_new(:)=rad(:)
1189 245023 : 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 1530887 : do iq=2,mqgrid
1224 1530383 : arg=2.d0*pi*qgrid(iq)
1225 :
1226 : ! Set up integrand
1227 1020372390 : do ir=1,mmax_new
1228 1020372390 : work(ir)=sin(arg*rad_new(ir))*rvlpz_new(ir)
1229 : end do
1230 :
1231 : ! Do integral from zero to rad(mmax)
1232 1530383 : call ctrap(mmax_new,work,amesh_new,result)
1233 :
1234 : ! Store q^2 v(q)
1235 1530887 : 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 504 : qmesh=qgrid(2)-qgrid(1)
1241 : yp1=(-50.d0*q2vq(1)+96.d0*q2vq(2)-72.d0*q2vq(3)+32.d0*q2vq(4)&
1242 504 : & -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 504 : & -96.d0*q2vq(mqgrid-1)+50.d0*q2vq(mqgrid))/(24.d0*qmesh)
1245 :
1246 504 : ABI_FREE(work)
1247 504 : ABI_FREE(rad_new)
1248 504 : ABI_FREE(rvlpz_new)
1249 504 : ABI_FREE(rvlpz)
1250 :
1251 504 : 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 504 : 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 504 : real(dp),allocatable :: rad_new(:),vpspll_new(:,:),work(:,:),work2(:)
1320 : ! *************************************************************************
1321 :
1322 : ! Find r mesh spacing necessary for accurate integration at qmax
1323 504 : amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid))
1324 :
1325 : ! Choose submultiple of input mesh
1326 504 : mesh_mult=int(amesh/amesh_new) + 1
1327 504 : mmax_new=mesh_mult*(mmax-1)+1
1328 504 : amesh_new=amesh/dble(mesh_mult)
1329 :
1330 1512 : ABI_MALLOC(rad_new,(mmax_new))
1331 2016 : ABI_MALLOC(vpspll_new,(mmax_new,lnmax))
1332 :
1333 504 : if (mesh_mult == 1) then
1334 250076 : 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 504 : dri = one/amesh
1352 341050 : do irn=1,mmax_new
1353 : ! index to find bracketing input mesh points
1354 341050 : 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 3037214 : ir=irn
1390 : iln0=0
1391 3037214 : do ilmn=1,lmnmax
1392 2787564 : iln=indlmn(5,ilmn)
1393 3037214 : if (iln>iln0) then
1394 1883600 : iln0=iln
1395 1883600 : tv = vpspll(ir,iln)
1396 1883600 : if(abs(tv)>tol10) then
1397 549809 : vpspll_new(irn,iln)=tv
1398 549809 : mvpspll=irn
1399 : else
1400 1333791 : vpspll_new(irn,iln)=zero
1401 : end if
1402 : end if
1403 : end do
1404 :
1405 : end if
1406 : end do ! irn
1407 :
1408 2016 : ABI_MALLOC(work, (mvpspll,lnmax))
1409 :
1410 : ! Loop over q values
1411 1531391 : do iq=1,mqgrid
1412 1530887 : arg=2.d0*pi*qgrid(iq)
1413 :
1414 : ! Set up integrands
1415 311188102 : do ir=1,mvpspll
1416 309657215 : call sbf8(lmax+1,arg*rad_new(ir),sb_out)
1417 309657215 : iln0=0
1418 3263786839 : do ilmn=1,lmnmax
1419 2952598737 : iln=indlmn(5,ilmn)
1420 3262255952 : if (iln>iln0) then
1421 2127794894 : iln0=iln
1422 2127794894 : ll=indlmn(1,ilmn)
1423 2127794894 : 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 16490849 : do ilmn=1,lmnmax
1431 14959458 : iln=indlmn(5,ilmn)
1432 16490345 : if (iln>iln0) then
1433 10442953 : iln0=iln
1434 10442953 : call ctrap(mvpspll,work(1,iln),amesh_new,result)
1435 10442953 : ffspl(iq,1,iln)=result
1436 : end if
1437 : end do
1438 : end do ! iq mesh
1439 :
1440 : ! Fit splines for form factors
1441 1512 : ABI_MALLOC(work2,(mqgrid))
1442 504 : qmesh=qgrid(2)-qgrid(1)
1443 :
1444 504 : iln0=0
1445 5464 : do ilmn=1,lmnmax
1446 4960 : iln=indlmn(5,ilmn)
1447 5464 : if (iln>iln0) then
1448 3455 : 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 3455 : & +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 3455 : & +50.d0*ffspl(mqgrid,1,iln))/(24.d0*qmesh)
1455 :
1456 3455 : call spline(qgrid,ffspl(1,1,iln),mqgrid,yp1,ypn,ffspl(1,2,iln))
1457 : end if
1458 : end do
1459 :
1460 504 : ABI_FREE(rad_new)
1461 504 : ABI_FREE(vpspll_new)
1462 504 : ABI_FREE(work)
1463 504 : ABI_FREE(work2)
1464 :
1465 504 : 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 : !!***
|
