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 149117 : 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 149117 : 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 149117 : cutoff=gsqcut*tolfix
167 149117 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
168 149117 : nfft=n1*n2*n3
169 149117 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
170 :
171 : ! Initialize container
172 447351 : ABI_MALLOC(gvec,(3,MAX(n1,n2,n3)))
173 447351 : ABI_MALLOC(gpq,(nfft))
174 298234 : ABI_MALLOC(gpq2,(nfft))
175 298234 : ABI_MALLOC(gcutoff,(nfft))
176 10654454576 : gcart(:) = zero ; gpq = zero ; gpq2 = zero ; gcutoff = zero
177 :
178 : !Set the q point for calls from linear-response routines
179 149117 : 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 596468 : do ii=1,3
184 447351 : id(ii)=ngfft(ii)/2+2
185 10570430 : do ing=1,ngfft(ii)
186 10421313 : 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 149117 : mode='NONE'
193 :
194 149117 : if (icutcoul==0) mode='SPHERE'
195 149117 : if (icutcoul==1) mode='CYLINDER'
196 149117 : if (icutcoul==2) mode='SLAB'
197 149117 : if (icutcoul==3) mode='CRYSTAL'
198 149117 : if (icutcoul==4) mode='ERF'
199 149117 : if (icutcoul==5) mode='ERFC'
200 149117 : if (icutcoul==22) mode='SLAB_SR'
201 149117 : if (icutcoul==55) mode='PCM'
202 :
203 : !Print in log info about the cut-off method at every call:
204 149117 : enough = enough + 1
205 149117 : if (enough < 5) then
206 5385 : write(msg,'(3a)')ch10,' 1/G**2 cut-off applied in the following step : cutoff-mode = ',TRIM(mode)
207 5385 : 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 149117 : if (icutcoul==55) mode='CRYSTAL'
213 :
214 3535123 : do i3=1,n3
215 : ! Precompute some products that do not depend on i2 and i1
216 3386006 : gs3=gvec(3,i3)*gvec(3,i3)*gmet(3,3)
217 3386006 : gvecgm23=gvec(3,i3)*gmet(2,3)*2
218 3386006 : gvecgm13=gvec(3,i3)*gmet(1,3)*2
219 :
220 93403905 : do i2=1,n2
221 89868782 : i23=n1*(i2-1 + n2*(i3-1))
222 89868782 : gs2=gs3+ gvec(2,i2)*(gvec(2,i2)*gmet(2,2)+gvecgm23)
223 89868782 : gvecgm12=gvec(2,i2)*gmet(1,2)*2
224 89868782 : gvecg2p3=gvecgm13+gvecgm12
225 3644689941 : do i1=1,n1
226 3551435153 : ii=i1+i23
227 3551435153 : gpq(ii)=gs2+gvec(1,i1)*(gvec(1,i1)*gmet(1,1)+gvecg2p3)
228 3641303935 : if(gpq(ii)>=tol4) then
229 3551295829 : gpq2(ii) = piinv/gpq(ii)
230 : end if
231 : end do
232 : end do
233 : end do
234 :
235 149134 : 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 3532416050 : 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 149117 : 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 149117 : ABI_FREE(gvec)
650 149117 : ABI_FREE(gpq)
651 149117 : ABI_FREE(gpq2)
652 : ! ABI_FREE(gcutoff)
653 :
654 149117 : 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 : !!***
|
