Line data Source code
1 : !!****m* ABINIT/m_cutoff_cylinder
2 : !! NAME
3 : !! m_cutoff_cylinder
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! SOURCE
8 :
9 : #if defined HAVE_CONFIG_H
10 : #include "config.h"
11 : #endif
12 :
13 : #include "abi_common.h"
14 :
15 : module m_cutoff_cylinder
16 :
17 : use defs_basis
18 : use m_abicore
19 : use m_errors
20 : use m_xmpi
21 : use m_splines
22 : use m_sort
23 :
24 : use m_fstrings, only : sjoin, itoa
25 : use m_geometry, only : normv, metric
26 : use m_bessel, only : CALJY0, CALJY1, CALCK0, CALCK1
27 : use m_numeric_tools, only : OPERATOR(.x.), quadrature
28 : use m_paw_numeric, only : paw_jbessel
29 :
30 : implicit none
31 :
32 : private
33 : !!***
34 :
35 : public :: cutoff_cylinder, K0cos
36 :
37 : !integer,public,parameter :: CYLINDER_BEIGI = 1
38 : !integer,public,parameter :: CYLINDER_ROZZI = 2
39 : !!***
40 :
41 : ! private variables used for the integration needed by the cylindrical case.
42 : integer,save :: npts_,ntrial_,qopt_
43 : real(dp),save :: ha_,hb_,r0_
44 : real(dp),save :: qpg_perp_,qpg_para_,qpgx_,qpgy_
45 : real(dp),save :: zz_,xx_, rho_
46 : real(dp),save :: hcyl_,rcut_,accuracy_
47 :
48 : CONTAINS
49 : !!***
50 :
51 : !----------------------------------------------------------------------
52 :
53 : !!****f* ABINIT/cutoff_cylinder
54 : !! NAME
55 : !! cutoff_cylinder
56 : !!
57 : !! FUNCTION
58 : !! Calculate the Fourier components of an effective Coulomb interaction
59 : !! zeroed outside a finite cylindrical region. Two methods are implemented:
60 : !!
61 : !! method==1: The interaction in the (say) x-y plane is truncated outside the Wigner-Seitz
62 : !! cell centered on the wire in the x-y plane. The interaction has infinite
63 : !! extent along the z axis and the Fourier transform is singular only at the Gamma point.
64 : !! Only orthorombic Bravais lattices are supported.
65 : !! method==2: The interaction is truncated outside a cylinder of radius rcut. The cylinder has finite
66 : !! extent along z. No singularity occurs.
67 : !!
68 : !! INPUTS
69 : !! boxcenter(3)= center of the wire in the x-y axis
70 : !! qpt(3)= q-point
71 : !! ng=number of G vectors
72 : !! gvec(3,ng)=G vectors in reduced coordinates
73 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
74 : !! method=1 for Beigi approach (infinite cylinder with interaction truncated outside the W-S cell)
75 : !! 2 for Rozzi method (finite cylinder)
76 : !! comm=MPI communicator.
77 : !!
78 : !! OUTPUT
79 : !! vc_cut(ng)= Fourier components of the effective Coulomb interaction
80 : !!
81 : !! SOURCE
82 :
83 112 : subroutine cutoff_cylinder(qpt, ng, gvec, rcut, hcyl, pdir, boxcenter, rprimd, vc_cut, method, comm)
84 :
85 : !Arguments ------------------------------------
86 : !scalars
87 : integer,intent(in) :: ng,method,comm
88 : real(dp),intent(in) :: rcut,hcyl
89 : !arrays
90 : integer,intent(in) :: gvec(3,ng),pdir(3)
91 : real(dp),intent(in) :: boxcenter(3),qpt(3),rprimd(3,3)
92 : real(dp),intent(out) :: vc_cut(ng)
93 :
94 : !Local variables-------------------------------
95 : !scalars
96 : integer,parameter :: N0=1000
97 : integer :: ig,igs,ierr, my_rank, nproc
98 : real(dp) :: j0,j1,k0,k1,qpg2,qpg_xy,tmp
99 : real(dp) :: qpg_z,quad,rcut2,hcyl2,c1,c2,ucvol,SMALL
100 : logical :: q_is_zero
101 : character(len=500) :: msg
102 : !arrays
103 : real(dp) :: qpg(3),b1(3),b2(3),b3(3),gmet(3,3),rmet(3,3),gprimd(3,3),qc(3),gcart(3)
104 : !************************************************************************
105 :
106 : ABI_UNUSED(pdir)
107 : ABI_UNUSED(boxcenter)
108 :
109 : ! ===================================================
110 : ! === Setup for the quadrature of matrix elements ===
111 : ! ===================================================
112 112 : qopt_ =6 ! Quadrature method, see quadrature routine.
113 112 : ntrial_ =30 ! Max number of attempts.
114 112 : accuracy_=0.001 ! Fractional accuracy required.
115 112 : npts_ =6 ! Initial number of point (only for Gauss-Legendre method).
116 112 : SMALL =tol4 ! Below this value (q+G)_i is treated as zero.
117 112 : rcut_ =rcut ! Radial cutoff, used only if method==2
118 112 : hcyl_ =hcyl ! Lenght of cylinder along z, only if method==2
119 :
120 : !write(msg,'(3a,2(a,i5,a),a,f8.5)')ch10,&
121 : ! ' cutoff_cylinder: Info on the quadrature method : ',ch10,&
122 : ! ' Quadrature scheme = ',qopt_,ch10,&
123 : ! ' Max number of attempts = ',ntrial_,ch10,&
124 : ! ' Fractional accuracy = ',accuracy_
125 : !call wrtout(std_out, msg)
126 :
127 : ! From reduced to Cartesian coordinates.
128 112 : call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
129 448 : b1(:) =two_pi*gprimd(:,1)
130 448 : b2(:) =two_pi*gprimd(:,2)
131 448 : b3(:) =two_pi*gprimd(:,3)
132 :
133 448 : qc = b1(:)*qpt(1) + b2(:)*qpt(2) + b3(:)*qpt(3)
134 :
135 112 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
136 :
137 : ! ================================================
138 : ! === Different approaches according to method ===
139 : ! ================================================
140 1640 : vc_cut = zero
141 :
142 112 : select case (method)
143 :
144 : case (1)
145 : ! Infinite cylinder, interaction is zeroed outside the Wigner-Seitz cell.
146 : ! NB: Beigi's expression holds only if the BZ is sampled only along z.
147 : !call wrtout(std_out, 'cutoff_cylinder: Using Beigi''s Infinite cylinder')
148 :
149 318 : if (ANY(qc(1:2) > SMALL)) then
150 : write(msg,'(5a)')&
151 0 : ' found q-points with non zero components in the X-Y plane. ',ch10,&
152 0 : ' This is not allowed, see Notes in cutoff_cylinder.F90. ',ch10,&
153 0 : ' ACTION: Modify the q-point sampling. '
154 0 : ABI_ERROR(msg)
155 : end if
156 :
157 : ! Check if Bravais lattice is orthorombic and parallel to the Cartesian versors.
158 : ! In this case the intersection of the WS cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_
159 : if ((ANY(ABS(rprimd(2:3, 1)) > tol6)) .or. &
160 954 : (ANY(ABS(rprimd(1:3:2,2)) > tol6)) .or. &
161 : (ANY(ABS(rprimd(1:2, 3)) > tol6))) then
162 0 : ABI_ERROR('Bravais lattice should be orthorombic and parallel to the Cartesian versors')
163 : end if
164 :
165 424 : ha_ = half*NORM2(rprimd(:,1))
166 424 : hb_ = half*NORM2(rprimd(:,2))
167 106 : r0_ = MIN(ha_,hb_)/N0
168 :
169 : ! For each (q,G) pair evaluate the integral defining the Coulomb cutoff.
170 : ! NB: the code assumes that all q-vectors are non zero and q_xy/=0.
171 106 : igs=1
172 : ! Skip singularity at Gamma, it will be treated "by hand" in csigme.
173 106 : q_is_zero = (normv(qpt, gmet, 'G') < tol4)
174 :
175 920 : do ig=igs,ng
176 814 : if (mod(ig, nproc) /= my_rank) cycle ! MPI parallelism
177 :
178 3256 : gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
179 3256 : qpg(:)=qc(:)+gcart(:)
180 814 : qpgx_=qpg(1); qpgy_=qpg(2); qpg_para_=ABS(qpg(3))
181 : !write(std_out,*)"qpgx_=",qpgx_, "qpgy_=",qpgy_, "qpg_para=",qpg_para_
182 :
183 : ! Avoid singularity in K_0{qpg_para_\rho) by using a small q along the periodic dimension.
184 814 : if (q_is_zero .and. qpg_para_ < tol6) qpg_para_ = tol6
185 :
186 : ! Calculate $ 2\int_{WS} dxdy K_0{qpg_para_\rho) cos(x.qpg_x + y.qpg_y) $
187 : ! where WS is the Wigner-Seitz cell.
188 814 : tmp=zero
189 :
190 : ! Difficult part, integrate on a small cirle of radius r0 using spherical coordinates
191 : !call quadrature(K0cos_dth_r0,zero,r0_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
192 : !ABI_CHECK(ierr == 0, "Accuracy not reached")
193 : !write(std_out,'(i8,a,es14.6)')ig,' 1 ',quad
194 : !tmp=tmp+quad
195 : ! Add region with 0<=x<=r0 and y>=+-(SQRT(r0^2-x^2))since WS is rectangular
196 : !call quadrature(K0cos_dy_r0,zero,r0_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
197 : !ABI_CHECK(ierr == 0, "Accuracy not reached")
198 : !write(std_out,'(i8,a,es14.6)')ig,' 2 ',quad
199 : !tmp=tmp+quad
200 : ! Get the in integral in the rectangle with x>=r0, should be the easiest but sometimes has problems to converge
201 : !call quadrature(K0cos_dy,r0_,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
202 : !ABI_CHECK(ierr == 0, "Accuracy not reached")
203 : !write(std_out,'(i8,a,es14.6)')ig,' 3 ',quad
204 : !
205 : ! More stable method: midpoint integration with Romberg extrapolation ===
206 814 : call quadrature(K0cos_dy,zero,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
207 : !write(std_out,'(i8,a,es14.6)')ig,' 3 ',quad
208 814 : ABI_CHECK(ierr == 0, "Accuracy not reached in quadrature!")
209 :
210 : ! Store final result
211 : ! Factor two comes from the replacement WS -> (1,4) quadrant thanks to symmetries of the integrad.
212 814 : tmp = tmp+quad
213 920 : vc_cut(ig) = two*(tmp*two)
214 : end do ! ig
215 :
216 : case (2)
217 : ! Finite cylinder of length hcyl from Rozzi et al.
218 : ! TODO add check on hcyl value that should be smaller that 1/deltaq
219 6 : if (hcyl_ < zero) then
220 0 : write(msg,'(a,f8.4)')' Negative value for cylinder length hcyl_=',hcyl_
221 0 : ABI_BUG(msg)
222 : end if
223 :
224 6 : if (ABS(hcyl_) > tol12) then
225 : !write(std_out,'(2(a,f8.4))')' cutoff_cylinder: using finite cylinder of length= ',hcyl_,' rcut= ',rcut_
226 6 : hcyl2=hcyl_**2
227 6 : rcut2=rcut_**2
228 :
229 : ! No singularity occurs in finite cylinder, thus start from 1.
230 720 : do ig=1,ng
231 714 : if (mod(ig, nproc) /= my_rank) cycle ! MPI parallelism
232 :
233 2856 : gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
234 2856 : qpg(:)=qc(:)+gcart(:)
235 714 : qpg_para_=ABS(qpg(3)) ; qpg_perp_=SQRT(qpg(1)**2+qpg(2)**2)
236 :
237 720 : if (qpg_perp_ /= zero .and. qpg_para_ /= zero) then
238 : ! $ 4\pi\int_0^{R_c} d\rho\rho j_o(qpg_perp_.\rho)\int_0^hcyl dz\cos(qpg_para_*z)/sqrt(\rho^2+z^2) $
239 655 : call quadrature(F2,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
240 655 : ABI_CHECK(ierr == 0, "Accuracy not reached")
241 655 : vc_cut(ig) = four_pi*quad
242 :
243 59 : else if (qpg_perp_ == zero .and. qpg_para_ /= zero) then
244 : ! $ \int_0^h sin(qpg_para_.z)/\sqrt(rcut^2+z^2)dz $
245 14 : call quadrature(F3,zero,hcyl_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
246 14 : ABI_CHECK(ierr == 0, "Accuracy not reached")
247 :
248 14 : c1=one/qpg_para_**2-COS(qpg_para_*hcyl_)/qpg_para_**2-hcyl_*SIN(qpg_para_*hcyl_)/qpg_para_
249 14 : c2=SIN(qpg_para_*hcyl_)*SQRT(hcyl2+rcut2)
250 14 : vc_cut(ig) = four_pi*c1+four_pi*(c2-quad)/qpg_para_
251 :
252 45 : else if (qpg_perp_ /= zero .and. qpg_para_ == zero) then
253 : ! $ 4pi\int_0^rcut d\rho \rho J_o(qpg_perp_.\rho) ln((h+\sqrt(h^2+\rho^2))/\rho) $
254 44 : call quadrature(F4,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
255 44 : ABI_CHECK(ierr == 0, "Accuracy not reached")
256 44 : vc_cut(ig) = four_pi*quad
257 :
258 1 : else if (qpg_perp_ == zero .and. qpg_para_ == zero) then
259 : ! Use lim q+G --> 0
260 1 : vc_cut(ig) = two_pi*(-hcyl2+hcyl_*SQRT(hcyl2+rcut2)+rcut2*LOG((hcyl_+SQRT(hcyl_+SQRT(hcyl2+rcut2)))/rcut_))
261 :
262 : else
263 0 : ABI_BUG('You should not be here!')
264 : end if
265 :
266 : end do !ig
267 :
268 : else
269 : ! Infinite cylinder.
270 : !call wrtout(std_out, ' cutoff_cylinder: using Rozzi''s method with infinite cylinder ')
271 :
272 0 : do ig=1,ng
273 0 : if (mod(ig, nproc) /= my_rank) cycle ! MPI parallelism
274 :
275 0 : gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig)
276 0 : qpg(:)=qc(:)+gcart(:)
277 0 : qpg2 =DOT_PRODUCT(qpg,qpg)
278 0 : qpg_z =ABS(qpg(3)) ; qpg_xy=SQRT(qpg(1)**2+qpg(2)**2)
279 :
280 0 : if (qpg_z > SMALL) then
281 : ! Analytic expression.
282 0 : call CALCK0(qpg_z *rcut_, k0, 1)
283 0 : call CALJY1(qpg_xy*rcut_, j1, 0)
284 0 : call CALJY0(qpg_xy*rcut_, j0, 0)
285 0 : call CALCK1(qpg_z *rcut_, k1, 1)
286 0 : vc_cut(ig) = (four_pi/qpg2)*(one+rcut_*qpg_xy*j1*k0-qpg_z*rcut_*j0*k1)
287 : else
288 0 : if (qpg_xy > SMALL) then
289 : ! Integrate r*Jo(G_xy r)log(r) from 0 up to rcut_
290 0 : call quadrature(F5,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
291 0 : ABI_CHECK(ierr == 0, "Accuracy not reached")
292 0 : vc_cut(ig)=-four_pi*quad
293 : else
294 : ! Analytic expression
295 0 : vc_cut(ig)=-pi*rcut_**2*(two*LOG(rcut_)-one)
296 : end if
297 : end if
298 : end do ! ig
299 : end if !finite/infinite
300 :
301 : case default
302 112 : ABI_BUG(sjoin('Wrong value for method:',itoa(method)))
303 : end select
304 :
305 : ! Collect vc_cut on each core
306 112 : call xmpi_sum(vc_cut, comm, ierr)
307 :
308 112 : end subroutine cutoff_cylinder
309 : !!***
310 :
311 : !----------------------------------------------------------------------
312 :
313 4297455 : real(dp) function F1(rho)
314 :
315 : real(dp),intent(in) :: rho
316 :
317 : !Local variables-------------------------------
318 : integer,parameter :: order=0,ll=0
319 : real(dp) :: arg,bes,besp,bespp
320 : !************************************************************************
321 :
322 : ! F1(\rho;z)= \rho*j_o(qpg_perp_*\rho)/sqrt(\rho**2+z**2)
323 4297455 : arg=rho*qpg_perp_
324 4297455 : call paw_jbessel(bes,besp,bespp,ll,order,arg)
325 :
326 4297455 : if (zz_==zero) then
327 0 : F1=bes
328 : else
329 4297455 : F1=bes*rho/SQRT(rho**2+zz_**2)
330 : end if
331 :
332 4297455 : end function F1
333 : !!***
334 :
335 : !----------------------------------------------------------------------
336 :
337 53055 : function F2(xx)
338 :
339 : real(dp),intent(in) :: xx
340 : real(dp) :: F2
341 :
342 : !Local variables-------------------------------
343 : !scalars
344 : integer :: ierr
345 : real(dp) :: intr
346 : !************************************************************************
347 :
348 53055 : zz_=xx
349 53055 : call quadrature(F1,zero,rcut_,qopt_,intr,ierr,ntrial_,accuracy_,npts_)
350 53055 : ABI_CHECK(ierr == 0, "Accuracy not reached")
351 :
352 53055 : F2=intr*COS(qpg_para_*xx)
353 :
354 53055 : end function F2
355 : !!***
356 :
357 : !----------------------------------------------------------------------
358 :
359 1134 : real(dp) pure function F3(xx)
360 :
361 : real(dp),intent(in) :: xx
362 :
363 : !************************************************************************
364 :
365 : ! F3(z)=z*\sin(qpg_para_*z)/\sqrt(rcut^2+z^2)
366 1134 : F3=xx*SIN(qpg_para_*xx)/SQRT(rcut_**2+xx**2)
367 :
368 1134 : end function F3
369 : !!***
370 :
371 : !----------------------------------------------------------------------
372 :
373 3564 : real(dp) function F4(rho)
374 :
375 : real(dp),intent(in) :: rho
376 :
377 : !Local variables-------------------------------
378 : integer,parameter :: order=0,ll=0
379 : real(dp) :: arg,bes,besp,bespp
380 : !************************************************************************
381 :
382 : ! $F4(rho)=\rho*j_o(qpg_perp_.\rho) \ln((hcyl+\sqrt(rho^2+hcyl^2))/\rho)$
383 3564 : if (ABS(rho)<tol12) then
384 : F4=zero
385 : else
386 3564 : arg=rho*qpg_perp_
387 3564 : call paw_jbessel(bes,besp,bespp,ll,order,arg)
388 3564 : F4=bes*rho*LOG((hcyl_+SQRT(rho**2+hcyl_**2))/rho)
389 : end if
390 :
391 3564 : end function F4
392 : !!***
393 :
394 : !----------------------------------------------------------------------
395 :
396 0 : real(dp) function F5(rho)
397 :
398 : real(dp),intent(in) :: rho
399 :
400 : !Local variables-------------------------------
401 : integer,parameter :: order = 0, ll = 0
402 : real(dp) :: arg,bes,besp,bespp
403 : !************************************************************************
404 :
405 : ! $F5(\rho)=\rho*j_o(G_perp\rho)log(\rho)$
406 0 : if (rho==0) then
407 : F5=zero
408 : else
409 0 : arg=rho*qpg_perp_
410 0 : call paw_jbessel(bes,besp,bespp,ll,order,arg)
411 0 : F5=bes*rho*LOG(rho)
412 : end if
413 :
414 0 : end function F5
415 : !!***
416 :
417 : !----------------------------------------------------------------------
418 :
419 5345190 : real(dp) function K0cos(yy)
420 :
421 : real(dp),intent(in) :: yy
422 :
423 : !Local variables-------------------------------
424 : real(dp) :: k0,rho,arg
425 : !************************************************************************
426 :
427 : ! K0cos(y)=K0(\rho*|qpg_z|)*COS(x.qpg_x+y*qpg_y)
428 5345190 : rho=SQRT(xx_**2+yy**2) ; arg=qpg_para_*rho
429 5345190 : call CALCK0(arg,k0,1)
430 5345190 : K0cos=k0*COS(qpgx_*xx_+qpgy_*yy)
431 :
432 5345190 : end function K0cos
433 : !!***
434 :
435 : !----------------------------------------------------------------------
436 :
437 65934 : real(dp) function K0cos_dy(xx)
438 :
439 : real(dp),intent(in) :: xx
440 :
441 : !Local variables-------------------------------
442 : integer :: ierr
443 : real(dp) :: quad
444 : !************************************************************************
445 :
446 : !! K0cos_dy(x)=\int_{-b/2}^{b/2} K0(|qpg_z|\rho)cos(x.qpg_x+y.qpg_y)dy$
447 65934 : xx_=xx
448 65934 : call quadrature(K0cos,-hb_,+hb_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
449 65934 : ABI_CHECK(ierr == 0, "Accuracy not reached")
450 :
451 65934 : K0cos_dy=quad
452 :
453 65934 : end function K0cos_dy
454 : !!***
455 :
456 :
457 :
458 :
459 : !----------------------------------------------------------------------
460 :
461 : real(dp) function K0cos_dy_r0(xx)
462 :
463 : real(dp),intent(in) :: xx
464 :
465 : !Local variables-------------------------------
466 : !scalars
467 : integer :: ierr
468 : real(dp) :: quad,yx
469 : !************************************************************************
470 :
471 : ! $ K0cos_dy_r0(x)= \int_{-b/2}^{-y(x)} K0(|qpg_z|\rho) cos(x.qpg_x+y.qpg_y)dy
472 : ! +\int_{y(x)}^{b/2} K0(|qpg_z|\rho)cos(x.qpg_x+y.qpg_y)dy$
473 : ! where y(x)=SQRT(r0^2-x^2) and x<=r0
474 : !
475 : xx_=xx; yx=SQRT(r0_**2-xx**2)
476 : call quadrature(K0cos,-hb_,-yx,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
477 : ABI_CHECK(ierr == 0, "Accuracy not reached in quadrature")
478 : K0cos_dy_r0=quad
479 :
480 : call quadrature(K0cos,+yx,+hb_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
481 : ABI_CHECK(ierr == 0, "Accuracy not reached in quadrature")
482 :
483 : K0cos_dy_r0=quad+K0cos_dy_r0
484 :
485 : end function K0cos_dy_r0
486 : !!***
487 :
488 : !----------------------------------------------------------------------
489 :
490 : real(dp) function K0cos_dth_r0(rho)
491 :
492 : real(dp),intent(in) :: rho
493 :
494 : !Local variables-------------------------------
495 : !scalars
496 : integer :: ierr
497 : real(dp) :: quad,arg,k0,tmp
498 :
499 : !************************************************************************
500 :
501 : ! $ K0cos_dth_r0(\rho)=
502 : ! \int_{0}^{2pi)} K0(|qpg_z|\rho)cos(\rho.cos(\theta).qpg_x+\rho.sin(\theta).qpg_y) d\theta $
503 : !
504 : ! where y(x)=SQRT(r0^2-x^2) and x<=r0
505 : !
506 : rho_=rho
507 : call quadrature(Fcos_th,zero,two_pi,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
508 : ABI_CHECK(ierr == 0, "Accuracy not reached in quadrature")
509 :
510 : arg=qpg_para_*rho_
511 : tmp=zero
512 : if (arg>tol6) then
513 : call CALCK0(arg,k0,1)
514 : tmp=k0*rho_
515 : end if
516 : K0cos_dth_r0=quad*tmp
517 :
518 : end function K0cos_dth_r0
519 : !!***
520 :
521 : !----------------------------------------------------------------------
522 :
523 : pure real(dp) function Fcos_th(theta)
524 :
525 : real(dp),intent(in) :: theta
526 :
527 : !************************************************************************
528 :
529 : ! $ Fcos_th(\theta)=rho*K0(\rho*|qpg_z|)*COS(\rho.COS(\theta).qpg_x+\rho.SIN/(\theta)*qpg_y) $
530 :
531 : !arg=qpg_para_*rho_
532 : !call CALCK0(arg,k0,1)
533 : !tmp=k0*rho_
534 : Fcos_th=COS(rho_*COS(theta)*qpgx_+rho_*SIN(theta)*qpgy_)
535 :
536 : end function Fcos_th
537 : !!***
538 :
539 : !----------------------------------------------------------------------
540 :
541 : !the following functions should be used to deal with the singularity in the Cylindrical cutoff
542 : !TODO Not yet used and indeed are still private
543 :
544 : function K0fit(mq,nn) result(vals)
545 :
546 : integer,intent(in) :: nn
547 : real(dp),intent(in) :: mq
548 : real(dp) :: vals(nn)
549 :
550 : !Local variables-------------------------------
551 : !scalars
552 : integer :: ii
553 : real(dp) :: mqh
554 : !arrays
555 : real(dp),parameter :: cc(7)=(/-0.57721566,0.42278420,0.23069756, &
556 : 0.03488590,0.00262698,0.00010750,0.00000740/)
557 : ! *************************************************************************
558 :
559 : if (nn>8.or.nn<1) then
560 : ABI_ERROR("nn>8.or.nn<1 not implemented")
561 : end if
562 :
563 : ! === Eq 9.8.5 in Abramovitz ===
564 : vals(1)=-LOG(mq*half)*I0(mq)
565 : mqh=mq*half
566 : do ii=2,nn
567 : vals(ii)=cc(ii-1)*mqh**(2*(ii-2))
568 : end do
569 :
570 : end function K0fit
571 :
572 : real(dp) function K0fit_int(mq,par,nn) result(integ)
573 :
574 : integer,intent(in) :: nn
575 : real(dp),intent(in) :: mq
576 : real(dp),intent(in) :: par(nn)
577 :
578 : !Local variables-------------------------------
579 : !scalars
580 : integer :: ii,aa
581 : real(dp) :: mqh
582 : !arrays
583 : real(dp),parameter :: cc(7)=(/-0.57721566,0.42278420,0.23069756,&
584 : & 0.03488590,0.00262698,0.00010750,0.00000740/)
585 : ! *************************************************************************
586 :
587 : if (nn>8.or.nn<1) then
588 : ABI_ERROR("nn>8.or.nn<1 not implemented")
589 : end if
590 :
591 : mqh=mq*half
592 : integ=-par(1)*int_I0ln(mqh)
593 : ! primitive of polynomial \sum_0^{N/2} cc_{2i} (x/2)^{2*i}
594 : do ii=2,nn
595 : aa=(2*(ii-1)+1)
596 : integ=integ+par(ii)*two*cc(ii-1)*(mqh**aa)/aa
597 : end do
598 :
599 : end function K0fit_int
600 :
601 : real(dp) function I0(xx)
602 :
603 : real(dp),intent(in) :: xx
604 :
605 : !Local variables-------------------------------
606 : real(dp) :: tt
607 :
608 : ! *************************************************************************
609 :
610 : ! Eq 9.8.1 of Abramovitz, entering the expansion of K0 -->0
611 : ! Expansion holds for |x|<3.75, Error<1.6*10D-07
612 : tt=xx/3.75
613 : I0=one+3.5156229*tt**2+3.0899424*tt**4 +1.2067492*tt**6 &
614 : +0.2659732*tt**8+0.0360768*tt**10+0.0045813*tt**12
615 : end function I0
616 :
617 : ! Primitive of x^m Ln(x) for m/=-1
618 : real(dp) function int_xmln(xx,mm) result(res)
619 :
620 : integer,intent(in) :: mm
621 : real(dp),intent(in) :: xx
622 :
623 : ! *********************************************************************
624 :
625 : if (mm==-1) then
626 : ABI_BUG('invalid value for mm')
627 : end if
628 :
629 : if (xx<=zero) then
630 : ABI_BUG(' invalid value for xx')
631 : end if
632 :
633 : res= (xx**(mm+1))/(mm+1) * (LOG(xx) - one/(mm+1))
634 :
635 : end function int_xmln
636 :
637 : ! Primitive function of ln(x/2)*I0(x) = sum_0^{N/2} 2^{2s+1} c_{2s} T(x/2,2s)
638 : ! where T(x,s)=\int x^s ln(x)dx
639 : real(dp) function int_I0ln(xx) result(res)
640 :
641 : !Arguments ------------------------------------
642 : real(dp),intent(in) :: xx
643 :
644 : !Local variables-------------------------------
645 : real(dp) :: yy
646 : ! *********************************************************************
647 :
648 : yy=xx*half
649 : res = ( one*2 *int_xmln(yy,0) &
650 : & +3.5156229*2**3 *int_xmln(yy,2) &
651 : & +3.0899424*2**5 *int_xmln(yy,4) &
652 : & +1.2067492*2**7 *int_xmln(yy,6) &
653 : & +0.2659732*2**9 *int_xmln(yy,8) &
654 : & +0.0360768*2**11*int_xmln(yy,10) &
655 : & +0.0045813*2**13*int_xmln(yy,12) &
656 : & )
657 :
658 : end function int_I0ln
659 : !!***
660 :
661 : end module m_cutoff_cylinder
662 : !!***
|
