Branch data Line data Source code
1 : : !!****m* ABINIT/m_gtermcutoff
2 : : !! NAME
3 : : !! m_gtermcutoff
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1999-2026 ABINIT group ()
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.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_gtermcutoff
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : ! use m_xmpi
28 : : use m_fstrings, only : sjoin, itoa
29 : : use m_profiling_abi, only : abimem_record
30 : : ! use defs_abitypes, only : MPI_type
31 : : use m_bessel, only : CALJY0, CALJY1, CALCK0, CALCK1
32 : : use m_numeric_tools, only : arth, l2norm, OPERATOR(.x.),quadrature
33 : : use m_paw_numeric, only : paw_jbessel
34 : : use m_geometry, only : normv, metric
35 : :
36 : : implicit none
37 : :
38 : : private
39 : : !!***
40 : :
41 : : !!****t* m_gtermcutoff/gtermcut_t
42 : : !! NAME
43 : : !! gtermcut_t
44 : : !!
45 : : !! FUNCTION
46 : : !!
47 : : !! SOURCE
48 : :
49 : : !!! type,public :: gtermcut_t
50 : :
51 : : !!! integer :: nfft
52 : : !!! ! Number of points in FFT grid
53 : :
54 : : !!! integer :: ng
55 : : !!! ! Number of G-vectors
56 : :
57 : : !!! real(dp) :: ucvol
58 : : !!! ! Volume of the unit cell
59 : :
60 : : !!! ! integer :: periodic_dir(3)
61 : : !!! ! 1 if the system is periodic along this direction
62 : :
63 : : !!! ! real(dp) :: boxcenter(3)
64 : : !!! ! 1 if the point in inside the cutoff region 0 otherwise
65 : : !!! ! Reduced coordinates of the center of the box (input variable)
66 : :
67 : : !!! real(dp) :: rprimd(3,3)
68 : : !!! ! Lattice vectors in real space.
69 : :
70 : : !!! real(dp),allocatable :: gtermcuoff(:)
71 : : !!! ! gtermcuoff(nfft)
72 : : !!! ! G cut-off array on the FFT grid
73 : :
74 : : !!! end type gtermcut_t
75 : :
76 : : public :: termcutoff
77 : : !!***
78 : : ! private variables used for the integration needed by the cylindrical case.
79 : : integer,save :: npts_,ntrial_,qopt_
80 : : real(dp),save :: ha_,hb_,hcyl_,r0_
81 : : real(dp),save :: gcart_para_,gcart_perp_,gcartx_,gcarty_
82 : : real(dp),save :: xx_,zz_,rcut_
83 : : real(dp),save :: accuracy_
84 : :
85 : :
86 : : contains
87 : : !!***
88 : :
89 : : !----------------------------------------------------------------------
90 : :
91 : : !!****f* ABINIT/termcutoff
92 : : !! NAME
93 : : !! termcutoff
94 : : !!
95 : : !! FUNCTION
96 : : !! Apply a cut-off term to the 1/G**2-like terms that appears throughout
97 : : !! the code at the ground-state level as follows: Ewald, NC-PSP, Hartree.
98 : : !!
99 : : !! INPUTS
100 : : !! gsqcut = cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential) (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
101 : : !! icutcoul = Information about the cut-off
102 : : !! ngfft(18) = Information on the (fine) FFT grid used for the density.
103 : : !! nkpt = Number of k-points in the Brillouin zone
104 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
105 : : !! vcutgeo(3)= Info on the orientation and extension of the cutoff region.
106 : : !!
107 : : !! OUTPUT
108 : : !! gcutoff = Cut-off term applied to 1/G**2 terms
109 : : !!
110 : : !! NOTES
111 : : !! 1. In order to incur minimal changes in some portions of the code
112 : : !! where a cut-off is needed to be applied, one can work only with
113 : : !! the cut-off part of the Coulomb potential, unlike what is done
114 : : !! in barevcoul module.
115 : : !! 2. Fock term has its own legacy cut-off for the moment.
116 : : !!
117 : : !! SOURCE
118 : :
119 : 150430 : subroutine termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo,qpt) !optional arguments
120 : :
121 : : !Arguments ------------------------------------
122 : : !scalars
123 : : integer,intent(in) :: icutcoul, nkpt
124 : : real(dp),intent(in) :: gsqcut,rcut
125 : :
126 : : !arrays
127 : : integer,intent(in) :: ngfft(18)
128 : : real(dp),intent(in) :: rprimd(3,3),vcutgeo(3)
129 : : real(dp),optional,intent(in) :: qpt(3)
130 : :
131 : : !Local variables-------------------------------
132 : : !scalars
133 : : integer,parameter :: N0=1000
134 : : integer,save :: enough
135 : : integer :: i1,i2,i23,i3,ierr,id(3),ii,ig,ing
136 : : integer :: c1,c2,opt_cylinder
137 : : integer :: n1,n2,n3,nfft
138 : : integer :: test,opt_slab !opt_cylinder
139 : : real(dp) :: alpha_fac, ap1sqrt, log_alpha
140 : : real(dp) :: cutoff,rcut_loc,rcut2,check,rmet(3,3)
141 : : real(dp) :: gvecg2p3,gvecgm12,gvecgm13,gvecgm23,gs2,gs3
142 : : real(dp) :: gcart_para,gcart_perp,gcart_x,gcart_y,gcart_z
143 : : real(dp) :: j0,j1,k0,k1
144 : : real(dp) :: odd2,quad,ucvol
145 : : real(dp) :: hcyl,hcyl2
146 : : real(dp),parameter :: tolfix=1.0000001_dp,tol999=999.0
147 : : character(len=50) :: mode
148 : : character(len=500) :: msg
149 : : ! type(gcut_t) :: gcut !
150 : :
151 : : !arrays
152 : : integer :: periodic_dir(3)
153 : : real(dp) :: a1(3),a2(3),a3(3),b1(3),b2(3),b3(3)
154 : : real(dp) :: gcart(3),gmet(3,3),gprimd(3,3)
155 : : real(dp) :: alpha(3),qpt_(3)
156 : 150430 : real(dp),allocatable :: gvec(:,:),gpq(:),gpq2(:)
157 : : real(dp),allocatable,intent(inout) :: gcutoff(:)
158 : :
159 : : ! === Save dimension and other useful quantities in vcut% ===
160 : : ! gcut%nfft = PRODUCT(ngfft(1:3)) ! Number of points in the FFT mesh.
161 : : ! gcut%ucvol = ucvol ! Unit cell volume.
162 : : ! gcut%rprimd = rprimd(:,:) ! Dimensional direct lattice.
163 : : ! gcut%vcutgeo = vcutgeo(:) ! Info on the orientation and extension of the cutoff region.
164 : : !
165 : : !Initialize a few quantities
166 : 150430 : cutoff=gsqcut*tolfix
167 : 150430 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
168 : 150430 : nfft=n1*n2*n3
169 : 150430 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
170 : :
171 : : ! Initialize container
172 [ - + + - : 451290 : ABI_MALLOC(gvec,(3,MAX(n1,n2,n3)))
+ - + - -
+ - + ]
173 [ - + + - : 451290 : ABI_MALLOC(gpq,(nfft))
+ - + - -
+ - + ]
174 [ - + + - : 300860 : ABI_MALLOC(gpq2,(nfft))
- + - + ]
175 [ - + + - : 300860 : ABI_MALLOC(gcutoff,(nfft))
- + - + -
+ ]
176 [ + + + + : 11332392280 : gcart(:) = zero ; gpq = zero ; gpq2 = zero ; gcutoff = zero
+ + ]
177 : :
178 : : !Set the q point for calls from linear-response routines
179 [ + + ]: 150430 : qpt_=zero; if (present(qpt)) qpt_=qpt
180 : :
181 : : !In order to speed the routine, precompute the components of gvectors
182 : : !Also check if the booked space was large enough...
183 [ + + ]: 601720 : do ii=1,3
184 : 451290 : id(ii)=ngfft(ii)/2+2
185 [ + + ]: 10722667 : do ing=1,ngfft(ii)
186 : 10572237 : gvec(ii,ing)=ing-(ing/id(ii))*ngfft(ii)-1 + qpt_(ii)
187 : : end do
188 : : end do
189 : :
190 : : ! Get the cut-off method info from the input file
191 : : ! Assign method to one of the available cases
192 : 150430 : mode='NONE'
193 : :
194 [ + + ]: 150430 : if (icutcoul==0) mode='SPHERE'
195 [ - + ]: 150430 : if (icutcoul==1) mode='CYLINDER'
196 [ + + ]: 150430 : if (icutcoul==2) mode='SLAB'
197 [ + + ]: 150430 : if (icutcoul==3) mode='CRYSTAL'
198 [ - + ]: 150430 : if (icutcoul==4) mode='ERF'
199 [ - + ]: 150430 : if (icutcoul==5) mode='ERFC'
200 [ + + ]: 150430 : if (icutcoul==22) mode='SLAB_SR'
201 [ + + ]: 150430 : if (icutcoul==55) mode='PCM'
202 : :
203 : : !Print in log info about the cut-off method at every call:
204 : 150430 : enough = enough + 1
205 [ + + ]: 150430 : if (enough < 5) then
206 : 5457 : write(msg,'(3a)')ch10,' 1/G**2 cut-off applied in the following step : cutoff-mode = ',TRIM(mode)
207 : 5457 : call wrtout(std_out,msg)
208 : : end if
209 : : !!!
210 : : !The Pick-Cohen-Martin cutoff is not a single factor on 1/G**2.
211 : : !It is therefore externally applied, not in this routine.
212 [ + + ]: 150430 : if (icutcoul==55) mode='CRYSTAL'
213 : :
214 [ + + ]: 3585431 : do i3=1,n3
215 : : ! Precompute some products that do not depend on i2 and i1
216 : 3435001 : gs3=gvec(3,i3)*gvec(3,i3)*gmet(3,3)
217 : 3435001 : gvecgm23=gvec(3,i3)*gmet(2,3)*2
218 : 3435001 : gvecgm13=gvec(3,i3)*gmet(1,3)*2
219 : :
220 [ + + ]: 96276276 : do i2=1,n2
221 : 92690845 : i23=n1*(i2-1 + n2*(i3-1))
222 : 92690845 : gs2=gs3+ gvec(2,i2)*(gvec(2,i2)*gmet(2,2)+gvecgm23)
223 : 92690845 : gvecgm12=gvec(2,i2)*gmet(1,2)*2
224 : 92690845 : gvecg2p3=gvecgm13+gvecgm12
225 [ + + ]: 3873539796 : do i1=1,n1
226 : 3777413950 : ii=i1+i23
227 : 3777413950 : gpq(ii)=gs2+gvec(1,i1)*(gvec(1,i1)*gmet(1,1)+gvecg2p3)
228 [ + + ]: 3870104795 : if(gpq(ii)>=tol4) then
229 : 3777273313 : gpq2(ii) = piinv/gpq(ii)
230 : : end if
231 : : end do
232 : : end do
233 : : end do
234 : :
235 : 150447 : SELECT CASE (TRIM(mode))
236 : :
237 : : CASE('SPHERE') ! Spherical cutoff
238 : :
239 : : ! Calculate rcut for each method
240 [ - + ]: 17 : if(rcut>tol4) then
241 : 0 : rcut_loc = rcut
242 : : else
243 : : ! Spencer-Alavi method
244 : 17 : rcut_loc = (three*nkpt*ucvol/four_pi)**(one/three)
245 : : endif
246 : :
247 [ + + ]: 1487177 : do ig=1,nfft
248 [ + + ]: 1487177 : if(abs(gpq(ig))<tol4) then
249 : 17 : gcutoff(ig)=0.0
250 : : else
251 : 1487143 : gcutoff(ig)=one-cos(rcut_loc*sqrt(four_pi/gpq2(ig)))
252 : : end if
253 : : end do
254 : :
255 : : CASE('CYLINDER')
256 : :
257 [ # # # # ]: 0 : test=COUNT(vcutgeo/=zero)
258 [ # # ]: 0 : ABI_CHECK(test==1,'Wrong cutgeo for cylinder')
259 : :
260 : : ! === From reduced to Cartesian coordinates ===
261 : 0 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
262 : :
263 [ # # # # ]: 0 : a1=rprimd(:,1); b1=two_pi*gprimd(:,1)
264 [ # # ]: 0 : a2=rprimd(:,2); b2=two_pi*gprimd(:,2)
265 [ # # ]: 0 : a3=rprimd(:,3); b3=two_pi*gprimd(:,3)
266 : :
267 : : !ntasks=nfft
268 : : !call xmpi_split_work(ntasks,comm,my_start,my_stop)
269 : :
270 : : !Calculate rcut for each method !
271 : : !
272 : :
273 : : ! * Check if Bravais lattice is orthorombic and parallel to the Cartesian versors.
274 : : ! In this case the intersection of the W-S cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_
275 : : if ( (ANY(ABS(rprimd(2:3, 1))>tol6)).or.&
276 [ # # # # : 0 : & (ANY(ABS(rprimd(1:3:2,2))>tol6)).or.&
# # # # #
# # # # #
# # ]
277 : : & (ANY(ABS(rprimd(1:2, 3))>tol6)) &
278 : : & ) then
279 : 0 : msg = ' Bravais lattice should be orthorombic and parallel to the cartesian versors '
280 : 0 : ABI_ERROR(msg)
281 : : end if
282 : :
283 : : ! === Beigi method is the default one, i.e infinite cylinder of radius rcut ===
284 : : ! * Negative values to use Rozzi method with finite cylinder of extent hcyl.
285 : 0 : opt_cylinder=1; hcyl=zero; periodic_dir(:)=0
286 [ # # ]: 0 : do ii=1,3
287 : 0 : check=vcutgeo(ii)
288 [ # # ]: 0 : if (ABS(check)>tol6) then
289 : 0 : periodic_dir(ii)=1
290 [ # # ]: 0 : if (check<zero) then ! use Rozzi's method.
291 [ # # ]: 0 : hcyl=ABS(check)*SQRT(SUM(rprimd(:,ii)**2))
292 : 0 : opt_cylinder=2
293 : : !Check to enter the infinite Rozzi treatment
294 [ # # ]: 0 : if(vcutgeo(3).le.-tol999) then
295 : 0 : hcyl=tol12
296 : : end if
297 : : end if
298 : : end if
299 : : end do
300 : :
301 : : ! Calculate rcut for each method
302 [ # # ]: 0 : if(rcut>tol4) then
303 : 0 : rcut_loc = rcut
304 : : else
305 [ # # ]: 0 : rcut_loc = half*SQRT(DOT_PRODUCT(a1,a1))
306 : : endif
307 : :
308 [ # # ]: 0 : if (opt_cylinder==1) then
309 [ # # # # : 0 : ABI_CHECK(ALL(periodic_dir == (/0,0,1/)),"The cylinder must be along the z-axis")
# # ]
310 : : end if
311 : :
312 : 0 : rcut_= rcut_loc
313 : :
314 : : ! ===================================================
315 : : ! === Setup for the quadrature of matrix elements ===
316 : : ! ===================================================
317 : 0 : qopt_ =6 ! Quadrature method, see quadrature routine.
318 : 0 : ntrial_ =30 ! Max number of attempts.
319 : 0 : accuracy_=0.001 ! Fractional accuracy required.
320 : 0 : npts_ =6 ! Initial number of point (only for Gauss-Legendre method).
321 : 0 : hcyl_ =hcyl ! Lenght of cylinder along z, only if method==2
322 : :
323 : 0 : write(msg,'(3a,2(a,i5,a),a,f8.5)')ch10,&
324 : 0 : ' cutoff_cylinder: Info on the quadrature method : ',ch10,&
325 : 0 : ' Quadrature scheme = ',qopt_,ch10,&
326 : 0 : ' Max number of attempts = ',ntrial_,ch10,&
327 : 0 : ' Fractional accuracy = ',accuracy_
328 : 0 : call wrtout(std_out,msg)
329 : :
330 : 0 : SELECT CASE (opt_cylinder)
331 : :
332 : : CASE(1)
333 : :
334 : : ! === Infinite cylinder, interaction is zeroed outside the Wigner-Seitz cell ===
335 : : ! * Beigi"s expression holds only if the BZ is sampled only along z.
336 : 0 : write(msg,'(2(a,f8.4))')' cutoff_cylinder: Using Beigi''s Infinite cylinder '
337 : 0 : call wrtout(std_out,msg)
338 : : ! * Check if Bravais lattice is orthorombic and parallel to the Cartesian versors.
339 : : ! In this case the intersection of the W-S cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_
340 : : if ( (ANY(ABS(rprimd(2:3, 1))>tol6)).or.&
341 [ # # # # : 0 : & (ANY(ABS(rprimd(1:3:2,2))>tol6)).or.&
# # # # #
# # # # #
# # ]
342 : : & (ANY(ABS(rprimd(1:2, 3))>tol6)) &
343 : : & ) then
344 : 0 : msg = ' Bravais lattice should be orthorhombic and parallel to the cartesian verctors '
345 : 0 : ABI_ERROR(msg)
346 : : end if
347 : :
348 [ # # ]: 0 : ha_=half*SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1)))
349 [ # # ]: 0 : hb_=half*SQRT(DOT_PRODUCT(rprimd(:,2),rprimd(:,2)))
350 : 0 : r0_=MIN(ha_,hb_)/N0
351 : :
352 [ # # ]: 0 : do i3=1,n3
353 [ # # ]: 0 : do i2=1,n2
354 : 0 : i23=n1*(i2-1 + n2*(i3-1))
355 [ # # ]: 0 : do i1=1,n1
356 : 0 : ii=i1+i23
357 : :
358 [ # # ]: 0 : gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3)
359 : 0 : gcartx_=gcart(1) ; gcarty_=gcart(2) ; gcart_para_=ABS(gcart(3))
360 [ # # ]: 0 : gpq(ii)=DOT_PRODUCT(gcart,gcart)
361 : :
362 : : ! Avoid singularity in K_0{gcart_para_\rho) by using a small g along the periodic dimension.
363 [ # # ]: 0 : if (gcart_para_<tol8) then
364 : 0 : gcart_para_ = tol8
365 : 0 : write(std_out,*)"setting gcart_para to=",gcart_para_
366 : : end if
367 : : !
368 : : ! * Calculate $ 2\int_{WS} dxdy K_0{gcart_para_\rho) cos(x.gcartx + y.gcarty) $
369 : : ! where WS is the Wigner-Seitz cell.
370 : : !tmp=zero
371 : : ! === More stable method: midpoint integration with Romberg extrapolation ===
372 : 0 : call quadrature(K0cos_dy,zero,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
373 : : !write(std_out,'(i8,a,es14.6)')ii,' 3 ',quad
374 [ # # ]: 0 : if (ierr/=0) then
375 : 0 : ABI_ERROR("Accuracy not reached")
376 : : end if
377 : : ! === Store final result ===
378 : : ! * Factor two comes from the replacement WS -> (1,4) quadrant thanks to symmetries of the integrad.
379 : : !tmp=tmp+quad
380 : 0 : gcutoff(ii)=quad*gpq(ii)/pi
381 : :
382 : : end do !i1
383 : : end do !i2
384 : : end do !i3
385 : :
386 : : CASE(2)
387 : :
388 : : ! === Finite cylinder of length hcyl, from Rozzi et al ===
389 : : ! TODO add check on hcyl value that should be smaller that 1/deltaq
390 [ # # ]: 0 : if (hcyl_<zero) then
391 : 0 : write(msg,'(a,f8.4)')' Negative value for cylinder length hcyl=',hcyl_
392 : 0 : ABI_BUG(msg)
393 : : end if
394 : :
395 [ # # ]: 0 : if (ABS(hcyl_)>tol12) then
396 : :
397 : 0 : write(msg,'(2(a,f8.4))')' cutoff_cylinder: using finite cylinder of length= ',hcyl,' rcut= ',rcut_loc
398 : 0 : call wrtout(std_out,msg)
399 : 0 : hcyl_=hcyl
400 : 0 : hcyl2=hcyl**2.0_dp
401 : 0 : rcut2=rcut_loc**2.0_dp
402 : :
403 [ # # ]: 0 : do i3=1,n3
404 [ # # ]: 0 : do i2=1,n2
405 : 0 : i23=n1*(i2-1 + n2*(i3-1))
406 [ # # ]: 0 : do i1=1,n1
407 : 0 : ii=i1+i23
408 : :
409 [ # # ]: 0 : gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3)
410 : 0 : gcart_para_=ABS(gcart(3)) ; gcart_perp_=SQRT(gcart(1)**2+gcart(2)**2)
411 [ # # ]: 0 : gpq(ii)=DOT_PRODUCT(gcart,gcart)
412 : :
413 [ # # # # ]: 0 : if (gcart_perp_/=zero.and.gcart_para_/=zero) then
414 : 0 : call quadrature(F2,zero,rcut_loc,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
415 [ # # ]: 0 : if (ierr/=0) then
416 : 0 : ABI_ERROR("Accuracy not reached")
417 : : end if
418 : :
419 : 0 : gcutoff(ii)=quad*gpq(ii)
420 : :
421 [ # # # # ]: 0 : else if (gcart_perp_==zero.and.gcart_para_/=zero) then
422 : :
423 : : ! $ \int_0^h sin(qpg_para_.z)/\sqrt(rcut^2+z^2)dz $
424 : 0 : call quadrature(F3,zero,hcyl,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
425 [ # # ]: 0 : if (ierr/=0) then
426 : 0 : ABI_ERROR("Accuracy not reached")
427 : : end if
428 : :
429 : 0 : c1=one/gcart_para_**2-COS(gcart_para_*hcyl_)/gcart_para_**2-hcyl_*SIN(gcart_para_*hcyl_)/gcart_para_
430 : 0 : c2=SIN(gcart_para_*hcyl_)*SQRT(hcyl2+rcut2)
431 : 0 : gcutoff(ii)=(c1+(c2-quad)/gcart_para_)*gpq(ii)
432 : :
433 [ # # # # ]: 0 : else if (gcart_perp_/=zero.and.gcart_para_==zero) then
434 : : ! $ 4pi\int_0^rcut d\rho \rho J_o(qpg_perp_.\rho) ln((h+\sqrt(h^2+\rho^2))/\rho) $
435 : 0 : call quadrature(F4,zero,rcut_loc,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
436 [ # # ]: 0 : if (ierr/=0) then
437 : 0 : ABI_ERROR("Accuracy not reached")
438 : : end if
439 : :
440 : 0 : gcutoff(ii)=quad*gpq(ii)
441 : :
442 [ # # # # ]: 0 : else if (gcart_perp_==zero.and.gcart_para_==zero) then
443 : : ! Use lim q+G --> 0
444 : 0 : gcutoff(ii)=zero
445 : : else
446 : 0 : ABI_BUG('You should not be here!')
447 : : end if
448 : :
449 : : end do !i1
450 : : end do !i2
451 : : end do !i3
452 : :
453 : : else
454 : :
455 : 0 : call wrtout(std_out,'Using Rozzi infinite cut-off cylinder method.')
456 : :
457 [ # # ]: 0 : do i3=1,n3
458 [ # # ]: 0 : do i2=1,n2
459 : 0 : i23=n1*(i2-1 + n2*(i3-1))
460 [ # # ]: 0 : do i1=1,n1
461 : 0 : ii=i1+i23
462 [ # # ]: 0 : gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3)
463 : 0 : gcart_x=gcart(1) ; gcart_y=gcart(2) ; gcart_z=ABS(gcart(3))
464 : 0 : gcart_perp_ = SQRT(gcart_x**2.0_dp+gcart_y**2.0_dp) ;
465 [ # # ]: 0 : gpq(ii)=DOT_PRODUCT(gcart,gcart)
466 : :
467 [ # # ]: 0 : if (gcart_z>tol4) then
468 : : ! === Analytic expression ===
469 : 0 : call CALJY1(gcart_perp_*rcut_loc,j1,0)
470 : 0 : call CALCK0(gcart_z*rcut_loc,k0,1)
471 : 0 : call CALJY0(gcart_perp_*rcut_loc,j0,0)
472 : 0 : call CALCK1(gcart_z*rcut_loc,k1,1)
473 : 0 : gcutoff(ii)=one+rcut_loc*gcart_perp_*j1*k0-rcut_loc*gcart_z*j0*k1
474 : : else
475 [ # # ]: 0 : if (gcart_perp_>tol4) then
476 : : ! === Integrate r*Jo(G_xy r)log(r) from 0 up to rcut_ ===
477 : 0 : call quadrature(F5,zero,rcut_loc,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
478 [ # # ]: 0 : if (ierr/=0) then
479 : 0 : ABI_ERROR("Accuracy not reached")
480 : : end if
481 : 0 : gcutoff(ii)= -quad*gpq(ii)
482 : : else
483 : 0 : gcutoff(ii)= zero !-pi*rcut_loc**2*(two*LOG(rcut_loc)-one)
484 : : end if
485 : : end if
486 : :
487 : : end do !i1
488 : : end do !i2
489 : : end do !i3
490 : : end if ! case 2 - selecting Rozzi
491 : :
492 : : CASE DEFAULT
493 [ # # # ]: 0 : ABI_BUG(sjoin('Wrong value for cylinder method:',itoa(opt_cylinder)))
494 : : END SELECT
495 : :
496 : : CASE('SLAB')
497 : :
498 [ + + + + ]: 1624 : test=COUNT(vcutgeo/=zero)
499 [ - + ]: 406 : ABI_CHECK(test==2,"Wrong vcutgeo")
500 : :
501 : : ! === From reduced to cartesian coordinates ===
502 [ + + ]: 1624 : a1=rprimd(:,1); b1=two_pi*gprimd(:,1)
503 [ + + ]: 1624 : a2=rprimd(:,2); b2=two_pi*gprimd(:,2)
504 [ + + + + ]: 2842 : a3=rprimd(:,3); b3=two_pi*gprimd(:,3)
505 : :
506 : : !SLAB Default - Beigi
507 : 406 : opt_slab=1; alpha(:)=zero
508 : : ! Otherwise use Rozzi's method
509 [ + + + + : 1573 : if (ANY(vcutgeo<zero) .or. rcut>tol8) opt_slab=2
+ + - + ]
510 : : periodic_dir(:)=0
511 [ + + ]: 1624 : do ii=1,3
512 : 1218 : check=vcutgeo(ii)
513 [ + + ]: 1624 : if (ABS(check)>zero) then
514 : : periodic_dir(ii)=1
515 : : !For Rozzi"s method
516 [ + + + + ]: 863 : if (check<zero) alpha(ii)=normv(check*rprimd(:,ii),rmet,'R')
517 : : end if
518 : : end do
519 : :
520 : : SELECT CASE (opt_slab)
521 : :
522 : : !CASE SLAB 1 - Beigi
523 : : CASE(1)
524 : :
525 : : ! Calculate rcut for each method !
526 [ + + ]: 1556 : rcut_loc = half*SQRT(DOT_PRODUCT(a3,a3))
527 : :
528 [ + + ]: 40319 : do i3=1,n3
529 [ + + ]: 675899 : do i2=1,n2
530 : 635580 : i23=n1*(i2-1 + n2*(i3-1))
531 [ + + ]: 14015670 : do i1=1,n1
532 : 13340160 : ii=i1+i23
533 [ + + ]: 53360640 : gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3)
534 : 13340160 : gcart_para=SQRT(gcart(1)**2+gcart(2)**2) ; gcart_perp = gcart(3)
535 [ + + + + ]: 13975740 : if(gcart_para<tol4.and.ABS(gcart_perp)<tol4) then
536 : : !if(gcart_para<tol12.and.ABS(gcart_perp)<tol12) then
537 : 267 : gcutoff(ii)=zero
538 : : else
539 : 13339893 : gcutoff(ii)=one-EXP(-gcart_para*rcut_loc)*COS(gcart_perp*rcut_loc)
540 : : end if
541 : : end do !i1
542 : : end do !i2
543 : : end do !i3
544 : :
545 : : !CASE SLAB 2 - Rozzi
546 : : CASE(2)
547 : :
548 : : !Set the cut-off radius
549 [ - + ]: 17 : if(rcut>tol4) then
550 : 0 : rcut_loc = rcut
551 : : else
552 [ + + ]: 68 : rcut_loc = half*SQRT(DOT_PRODUCT(a3,a3))
553 : : endif
554 : :
555 : : !In the case of finite, Rozzi's method provide another parameter
556 : : !for the cut-off: alpha
557 : : !!! ATT: alpha = L_x/L_y --> in-plane geometry dependence
558 : 17 : alpha_fac=SQRT(DOT_PRODUCT(a1,a1))/SQRT(DOT_PRODUCT(a2,a2))
559 : : ap1sqrt=SQRT(one+alpha_fac**2)
560 : : log_alpha=LOG((alpha_fac+ap1sqrt)*(one+ap1sqrt)/alpha_fac)
561 : :
562 [ + + ]: 2057 : do i3=1,n3
563 [ + + ]: 57137 : do i2=1,n2
564 : 55080 : i23=n1*(i2-1 + n2*(i3-1))
565 [ + + ]: 1544280 : do i1=1,n1
566 : 1487160 : ii=i1+i23
567 [ + + ]: 5948640 : gcart(:)=b1(:)*gvec(1,i1)+b2(:)*gvec(2,i2)+b3(:)*gvec(3,i3)
568 : 1487160 : gcart_para=SQRT(gcart(1)**2+gcart(2)**2) ; gcart_perp = gcart(3)
569 [ + + ]: 1542240 : if(gcart_para>tol4) then
570 : : gcutoff(ii)=one+EXP(-gcart_para*rcut_loc)*(gcart_perp/gcart_para*&
571 : 1485120 : & SIN(gcart_perp*rcut_loc)-COS(gcart_perp*rcut_loc))
572 : : else
573 [ + + ]: 2040 : if (ABS(gcart_perp)>tol4) then
574 : 2023 : gcutoff(ii)=one-COS(-gcart_perp*rcut_loc)-gcart_perp*rcut_loc*SIN(gcart_perp*rcut_loc)
575 : : ! gcutoff(ii)=one-COS(-gcart_perp*rcut_loc)-SIN(gcart_perp*rcut_loc) - Altered Rozzi's
576 : : else
577 : 17 : gcutoff(ii)=zero
578 : : endif
579 : : endif
580 : : end do !i1
581 : : end do !i2
582 : : end do !i3
583 : :
584 : : CASE DEFAULT
585 : 0 : write(msg,'(a,i3)')' Wrong value of slab method: ',opt_slab
586 [ + + - ]: 406 : ABI_BUG(msg)
587 : : END SELECT
588 : :
589 : : CASE('SLAB_SR')
590 : :
591 [ + + + + ]: 1268 : test=COUNT(vcutgeo/=zero)
592 [ - + ]: 317 : ABI_CHECK(test==2,"Wrong vcutgeo")
593 : :
594 [ + + ]: 28847 : do i3=1,n3
595 : 28530 : odd2=1-(-1)**(i3-1)
596 [ + + ]: 314147 : do i2=1,n2
597 : 285300 : i23=n1*(i2-1 + n2*(i3-1))
598 [ + + ]: 3166830 : do i1=1,n1
599 : 2853000 : ii=i1+i23
600 : 3138300 : gcutoff(ii)=odd2
601 : : end do
602 : : end do
603 : : end do
604 : :
605 : : CASE('ERF')
606 : :
607 : : ! Calculate rcut for each method ! Same as SPHERE
608 [ # # ]: 0 : if(rcut>tol4) then
609 : 0 : rcut_loc = rcut
610 : : else
611 : 0 : rcut_loc= (three*nkpt*ucvol/four_pi)**(one/three)
612 : : endif
613 : :
614 [ # # ]: 0 : do ig=1,nfft
615 [ # # ]: 0 : if(abs(gpq(ig))<tol4) then
616 : 0 : gcutoff(ig)=zero ! @Gamma: initialize quantity in each requiered routine
617 : : else !if(gpq(ig)<=cutoff) then
618 : 0 : gcutoff(ig)=exp(-pi/(gpq2(ig)*rcut_loc**2))
619 : : end if
620 : : end do !ig
621 : :
622 : : CASE('ERFC')
623 : :
624 : : ! Calculate rcut for each method ! Same as SPHERE
625 [ # # ]: 0 : if(rcut>tol4) then
626 : 0 : rcut_loc = rcut
627 : : else
628 : 0 : rcut_loc= (three*nkpt*ucvol/four_pi)**(one/three)
629 : : endif
630 : :
631 [ # # ]: 0 : do ig=1,nfft
632 [ # # ]: 0 : if(abs(gpq(ig))<tol4) then
633 : 0 : gcutoff(ig)=zero ! @Gamma: initialize quantity in each requiered routine
634 : : else
635 : 0 : gcutoff(ig)=one-exp(-pi/(gpq2(ig)*rcut_loc**2))
636 : : end if
637 : : end do !ig
638 : :
639 : : CASE('CRYSTAL')
640 [ + + ]: 3758396160 : gcutoff(:)=one ! Neutral cut-off
641 : : !write(msg,'(a)')'CRYSTAL method: no cut-off applied to G**2 while CRYSTAL method is implied!'
642 : : !ABI_WARNING(msg)
643 : : CASE DEFAULT
644 [ + - + + : 150430 : gcutoff=one ! Neutral cut-off
- - + - -
- ]
645 : : !write(msg,'(a)')'No cut-off applied to G**2!'
646 : : !ABI_WARNING(msg)
647 : : END SELECT
648 : :
649 [ - + ]: 150430 : ABI_FREE(gvec)
650 [ - + ]: 150430 : ABI_FREE(gpq)
651 [ - + ]: 150430 : ABI_FREE(gpq2)
652 : : ! ABI_FREE(gcutoff)
653 : :
654 : 150430 : end subroutine termcutoff
655 : : !!***
656 : :
657 : : !----------------------------------------------------------------------
658 : :
659 : 0 : function K0cos(yy)
660 : :
661 : : real(dp),intent(in) :: yy
662 : : real(dp) :: K0cos
663 : :
664 : : !Local variables-------------------------------
665 : : !scalars
666 : : real(dp) :: k0,rho,arg
667 : : !************************************************************************
668 : :
669 : : ! K0cos(y)=K0(\rho*|qpg_z|)*COS(x.qpg_x+y*qpg_y)
670 : 0 : rho=SQRT(xx_**2+yy**2) ; arg=gcart_para_*rho
671 : 0 : call CALCK0(arg,k0,1)
672 : 0 : K0cos=k0*COS(gcartx_*xx_+gcarty_*yy)
673 : :
674 : 0 : end function K0cos
675 : : !!***
676 : :
677 : : !----------------------------------------------------------------------
678 : :
679 : 0 : function K0cos_dy(xx)
680 : :
681 : : real(dp),intent(in) :: xx
682 : : real(dp) :: K0cos_dy
683 : : !Local variables-------------------------------
684 : : !scalars
685 : : integer :: ierr
686 : : real(dp) :: quad
687 : : !************************************************************************
688 : :
689 : : !! K0cos_dy(x)=\int_{-b/2}^{b/2} K0(|qpg_z|\rho)cos(x.qpg_x+y.qpg_y)dy$
690 : 0 : xx_=xx
691 : 0 : call quadrature(K0cos,-hb_,+hb_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
692 [ # # ]: 0 : if (ierr/=0) then
693 : 0 : ABI_ERROR("Accuracy not reached")
694 : : end if
695 : :
696 : 0 : K0cos_dy=quad
697 : :
698 : 0 : end function K0cos_dy
699 : : !!***
700 : :
701 : : !----------------------------------------------------------------------
702 : :
703 : 0 : function F1(rho)
704 : :
705 : : real(dp),intent(in) :: rho
706 : : real(dp) :: F1
707 : :
708 : : !Local variables-------------------------------
709 : : !scalars
710 : : integer,parameter :: order=0,ll=0
711 : : real(dp) :: arg,bes,besp,bespp
712 : : !************************************************************************
713 : :
714 : : !F1(\rho;z)= \rho*j_o(qpg_perp_*\rho)/sqrt(\rho**2+z**2)
715 : 0 : arg=rho*gcart_perp_
716 : 0 : call paw_jbessel(bes,besp,bespp,ll,order,arg)
717 : :
718 [ # # ]: 0 : if (zz_==zero) then
719 : 0 : F1=bes
720 : : else
721 : 0 : F1=bes*rho/SQRT(rho**2+zz_**2)
722 : : end if
723 : :
724 : 0 : end function F1
725 : : !!***
726 : :
727 : : !----------------------------------------------------------------------
728 : :
729 : 0 : function F2(xx)
730 : :
731 : : real(dp),intent(in) :: xx
732 : : real(dp) :: F2
733 : :
734 : : !Local variables-------------------------------
735 : : !scalars
736 : : integer :: ierr
737 : : real(dp) :: intr
738 : : !************************************************************************
739 : :
740 : 0 : zz_=xx
741 : 0 : call quadrature(F1,zero,rcut_,qopt_,intr,ierr,ntrial_,accuracy_,npts_)
742 [ # # ]: 0 : if (ierr/=0) then
743 : 0 : ABI_ERROR("Accuracy not reached")
744 : : end if
745 : :
746 : 0 : F2=intr*COS(gcart_para_*xx)
747 : :
748 : 0 : end function F2
749 : : !!***
750 : :
751 : : !----------------------------------------------------------------------
752 : :
753 : 0 : pure function F3(xx)
754 : :
755 : : real(dp),intent(in) :: xx
756 : : real(dp) :: F3
757 : : !************************************************************************
758 : :
759 : : ! F3(z)=z*\sin(qpg_para_*z)/\sqrt(rcut^2+z^2)
760 : 0 : F3=xx*SIN(gcart_para_*xx)/SQRT(rcut_**2+xx**2)
761 : :
762 : 0 : end function F3
763 : : !!***
764 : :
765 : : !----------------------------------------------------------------------
766 : :
767 : 0 : function F4(rho)
768 : :
769 : : real(dp),intent(in) :: rho
770 : : real(dp) :: F4
771 : :
772 : : !Local variables-------------------------------
773 : : !scalars
774 : : integer,parameter :: order=0,ll=0
775 : : real(dp) :: arg,bes,besp,bespp
776 : : !************************************************************************
777 : :
778 : : ! $F4(rho)=\rho*j_o(qpg_perp_.\rho) \ln((hcyl+\sqrt(rho^2+hcyl^2))/\rho)$
779 [ # # ]: 0 : if (ABS(rho)<tol12) then
780 : : F4=zero
781 : : else
782 : 0 : arg=rho*gcart_perp_
783 : 0 : call paw_jbessel(bes,besp,bespp,ll,order,arg)
784 : 0 : F4=bes*rho*LOG((hcyl_+SQRT(rho**2+hcyl_**2))/rho)
785 : : end if
786 : :
787 : 0 : end function F4
788 : : !!***
789 : :
790 : : !----------------------------------------------------------------------
791 : :
792 : 0 : function F5(rho)
793 : :
794 : : real(dp),intent(in) :: rho
795 : : real(dp) :: F5
796 : :
797 : : !Local variables-------------------------------
798 : : !scalars
799 : : integer,parameter :: order=0,ll=0
800 : : real(dp) :: arg,bes,besp,bespp
801 : : !************************************************************************
802 : :
803 : : ! $F5(\rho)=\rho*j_o(G_perp\rho)log(\rho)$
804 [ # # ]: 0 : if (rho==0) then
805 : : F5=zero
806 : : else
807 : 0 : arg=rho*gcart_perp_
808 : 0 : call paw_jbessel(bes,besp,bespp,ll,order,arg)
809 : 0 : F5=bes*rho*LOG(rho)
810 : : end if
811 : :
812 : 0 : end function F5
813 : : !!***
814 : :
815 : : end module m_gtermcutoff
816 : : !!***
|
