Branch data 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 : : !!***
|
