Line data Source code
1 : !!****m* ABINIT/m_barevcoul
2 : !! NAME
3 : !! m_barevcoul
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_barevcoul
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 :
29 : use m_fstrings, only : sjoin
30 : use defs_abitypes, only : MPI_type
31 : use m_numeric_tools, only : arth, l2norm, OPERATOR(.x.)
32 : use m_geometry, only : normv
33 : use m_crystal, only : crystal_t
34 : use m_fft, only : zerosym
35 : !use m_gsphere, only : gsphere_t
36 :
37 : ! Cut-off methods modules
38 : !use m_cutoff_sphere, only : cutoff_sphere
39 : !use m_cutoff_slab, only : cutoff_slab
40 : !use m_cutoff_cylinder, only : cutoff_cylinder
41 :
42 : implicit none
43 :
44 : private
45 : !!***
46 :
47 : !!****t* m_barevcoul/vcut_t
48 : !! NAME
49 : !! vcoul_t
50 : !!
51 : !! FUNCTION
52 : !!
53 : !! SOURCE
54 :
55 : type,public :: vcut_t
56 :
57 : integer :: nfft
58 : ! Number of points in FFT grid
59 :
60 : integer :: ng
61 : ! Number of G-vectors
62 :
63 : real(dp) :: alpha(3)
64 : ! Lenght of the finite slab
65 :
66 : real(dp) :: rcut
67 : ! Cutoff radius
68 :
69 : real(dp) :: i_sz
70 : ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2
71 :
72 : real(dp) :: hcyl
73 : ! Length of the finite cylinder along the periodic dimension
74 :
75 : real(dp) :: ucvol
76 : ! Volume of the unit cell
77 :
78 : character(len=50) :: mode
79 : ! String defining the cutoff mode, possible values are: sphere,cylinder,slab,crystal
80 :
81 : integer :: pdir(3)
82 : ! 1 if the system is periodic along this direction
83 :
84 : real(dp) :: boxcenter(3)
85 : ! 1 if the point in inside the cutoff region 0 otherwise
86 : ! Reduced coordinates of the center of the box (input variable)
87 :
88 : real(dp) :: vcutgeo(3)
89 : ! For each reduced direction gives the length of the finite system
90 : ! 0 if the system is infinite along that particular direction
91 : ! negative value to indicate that a finite size has to be used
92 :
93 : real(dp) :: rprimd(3,3)
94 : ! Lattice vectors in real space.
95 :
96 : real(dp),allocatable :: qibz(:,:)
97 : ! qibz(3,nqibz)
98 : ! q-points in the IBZ.
99 :
100 : real(dp),allocatable :: barev(:)
101 : ! barev(nfft)
102 : ! Bare Coulomb potential on the FFT grid
103 : ! A cut might be applied.
104 :
105 : end type vcut_t
106 :
107 : public :: barevcoul
108 : !!***
109 :
110 : contains
111 : !!***
112 :
113 : !----------------------------------------------------------------------
114 :
115 : !!****f* ABINIT/barevcoul
116 : !! NAME
117 : !! barevcoul
118 : !!
119 : !! FUNCTION
120 : !! Compute bare coulomb term in G-space on the FFT mesh i.e. 4pi/(G+q)**2 for a specified q-point
121 : !!
122 : !! INPUTS
123 : !! qpoint(3)=reduced coordinates for the phonon wavelength
124 : !! gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
125 : !! icutcoul=Option for the Coulomb potential cutoff technique
126 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
127 : !! izero=if 1, unbalanced components of V(q,g) are set to zero # Used by the PAW library
128 : !! nfft=Total number of FFT grid points.
129 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
130 : !! comm=MPI communicator.
131 : !!
132 : !! OUTPUT
133 : !! barev(nfft)=4pi/(G+q)**2, q+G=0 component is set carefully
134 : !!
135 : !! NOTES
136 : !! This routine operates on the full FFT mesh. DO NOT PASS MPI_TYPE
137 : !! One can easily implemente MPI-FFT by just calling this routine and then
138 : !! extracting the G-vectors treated by the node.
139 : !!
140 : !! SOURCE
141 :
142 214238 : subroutine barevcoul(rcut,icutcoul,qpoint,gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,izero,barev,shortrange)
143 :
144 : !Arguments ------------------------------------
145 : !scalars
146 : integer,intent(in) :: icutcoul,nfft,nkpt_bz,izero
147 : real(dp),intent(in) :: rcut,gsqcut,ucvol
148 : logical,intent(in),optional:: shortrange
149 : !arrays
150 : integer,intent(in) :: ngfft(18)
151 : integer :: ng!!!!
152 : real(dp),intent(in) :: qpoint(3)
153 : real(dp),intent(inout) :: gmet(3,3)
154 : real(dp),intent(inout) :: barev(nfft)
155 : !real(dp) :: a1(3),a2(3),a3(3)
156 : real(dp) :: b1(3),b2(3),b3(3),rmet(3,3) !,gprimd(3,3),
157 214238 : type(MPI_type) :: mpi_enreg !!!!
158 11140376 : type(crystal_t) :: Cryst !!!!
159 : !type(gsphere_t) :: Gsph
160 : !Local variables-------------------------------
161 : !scalars
162 : integer,parameter :: empty(3,3)=zero
163 : integer,parameter :: cplex1=1
164 : integer :: comm
165 : integer :: ii1,i1,i2,i23,i3,id1,id2,id3,icutcoul_local
166 : integer :: ig,ig1,ig2,ig3,ig1min,ig1max,ig2min,ig2max,ig3min,ig3max
167 : integer :: ii,ing,n1,n2,n3,npar,npt
168 : integer :: opt_cylinder,opt_slab,test
169 : integer :: qeq0,qeq05
170 : real(dp),parameter :: tolfix=1.000000001e0_dp ! Same value as the one used in hartre
171 : real(dp) :: check,step
172 : real(dp) :: cutoff,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,rcut0
173 : real(dp) :: bz_plane,dx,integ,q0_vol,q0_volsph
174 : character(len=500) :: msg
175 : !arrays
176 : integer :: id(3), gamma_pt(3,1)
177 214238 : real(dp),allocatable :: gq(:,:),gpq(:),gpq2(:)
178 214238 : real(dp),allocatable :: vcfit(:,:),xx(:),yy(:)
179 214238 : real(dp),allocatable :: cov(:,:),par(:),qfit(:,:),sigma(:),var(:),qcart(:,:)
180 214238 : type(vcut_t) :: vcut !!!!
181 : !
182 214238 : comm=mpi_enreg%comm_world
183 : !
184 : ! === Save dimension and other useful quantities in vcut% ===
185 856952 : vcut%nfft = PRODUCT(ngfft(1:3)) ! Number of points in the FFT mesh.
186 : ! ng and gvec are not used yet we don't want to use them without being defined
187 214238 : ng = -1
188 214238 : vcut%ng = ng ! Number of G-vectors in the Coulomb matrix elements.
189 214238 : vcut%rcut = rcut ! Cutoff radius for cylinder.
190 214238 : vcut%hcyl = zero ! Length of finite cylinder (Rozzi"s method, default is Beigi).
191 214238 : vcut%ucvol = ucvol ! Unit cell volume.
192 :
193 : !FBruneval: comment the definitions below since Cryst and dtset have never been initialized!
194 : !vcut%rprimd = Cryst%rprimd(:,:) ! Dimensional direct lattice.
195 : !vcut%boxcenter = dtset%boxcenter ! boxcenter at the moment is supposed to be at the origin.
196 : !vcut%vcutgeo = dtset%vcutgeo(:) ! Info on the orientation and extension of the cutoff region.
197 : !
198 : ! === Define geometry and cutoff radius (if used) ===
199 214238 : vcut%mode='NONE'
200 214238 : icutcoul_local=icutcoul
201 :
202 : ! for short-range exchange (e.g. HSE06), enforce ERFC
203 214238 : if (PRESENT(shortrange)) then
204 81956 : if (shortrange) then
205 : icutcoul_local=5
206 : end if
207 : end if
208 : ! -------------------------------------
209 :
210 132282 : if (icutcoul_local == 0) vcut%mode = 'SPHERE'
211 214238 : if (icutcoul_local == 1) vcut%mode = 'CYLINDER'
212 214238 : if (icutcoul_local == 2) vcut%mode = 'SLAB'
213 214238 : if (icutcoul_local == 3) vcut%mode = 'CRYSTAL'
214 214238 : if (icutcoul_local == 4) vcut%mode = 'ERF'
215 214238 : if (icutcoul_local == 5) vcut%mode = 'ERFC'
216 214238 : if (icutcoul_local == 6) vcut%mode = 'AUXILIARY_FUNCTION'
217 214238 : if (icutcoul_local == 7) vcut%mode = 'AUX_GB'
218 :
219 : !Initialize a few quantities
220 214238 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
221 214238 : cutoff = gsqcut * tolfix
222 907950866 : barev=zero
223 :
224 : !Some peculiar values of q: q=0 or q on the BZ edge
225 214238 : qeq0=0; if (qpoint(1)**2+qpoint(2)**2+qpoint(3)**2<1.d-15) qeq0=1
226 201255 : qeq05=0
227 : if (qeq0==0) then
228 201255 : if (abs(abs(qpoint(1))-half)<tol12.or.abs(abs(qpoint(2))-half)<tol12.or. &
229 102711 : & abs(abs(qpoint(3))-half)<tol12) qeq05=1
230 : end if
231 :
232 : !In order to speed the routine, precompute the components of g+q
233 : !Also check if the booked space was large enough...
234 :
235 642714 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
236 642714 : ABI_MALLOC(gpq,(nfft))
237 428476 : ABI_MALLOC(gpq2,(nfft))
238 :
239 856952 : do ii=1,3
240 642714 : id(ii)=ngfft(ii)/2+2
241 10749122 : do ing=1,ngfft(ii)
242 9892170 : ig=ing-(ing/id(ii))*ngfft(ii)-1
243 10534884 : gq(ii,ing)=ig+qpoint(ii)
244 : end do
245 : end do
246 214238 : ig1max=-1;ig2max=-1;ig3max=-1
247 214238 : ig1min=n1;ig2min=n2;ig3min=n3
248 :
249 214238 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
250 :
251 : ! Triple loop on each dimension
252 3506684 : do i3=1,n3
253 3292446 : ig3=i3-(i3/id3)*n3-1
254 : ! Precompute some products that do not depend on i2 and i1
255 3292446 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
256 3292446 : gqgm23=gq(3,i3)*gmet(2,3)*2
257 3292446 : gqgm13=gq(3,i3)*gmet(1,3)*2
258 55868246 : do i2=1,n2
259 52361562 : ig2=i2-(i2/id2)*n2-1
260 52361562 : i23=n1*(i2-1 +(n2)*(i3-1))
261 52361562 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
262 52361562 : gqgm12=gq(2,i2)*gmet(1,2)*2
263 52361562 : gqg2p3=gqgm13+gqgm12
264 :
265 960098190 : do i1=1,n1
266 907736628 : ii=i1+i23
267 907736628 : gpq(ii)= gs2 + gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
268 960098190 : if (gpq(ii)>=tol4) then
269 : ! gpq2 contains 4*pi / |q+G|**2
270 907723645 : gpq2(ii) = piinv/gpq(ii)
271 : end if
272 : end do
273 :
274 : !
275 : ! Next part looks for ig1min,ig1max that are needed by zerosym
276 : ! Do the test that eliminates the Gamma point outside of the inner loop
277 52361562 : ii1=1
278 52361562 : if (i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0) then
279 52361562 : ii1=2
280 : end if
281 :
282 : ! Final inner loop on the first dimension (note the lower limit)
283 963377653 : do i1=ii1,n1
284 907723645 : gs = gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
285 :
286 : !ii=i1+i23
287 :
288 960085207 : if (gs<=cutoff) then
289 : ! Identify min/max indexes (to cancel unbalanced contributions later)
290 : ! Count (q+g)-vectors with similar norm
291 297652559 : if ((qeq05==1).and.(izero==1)) then
292 6333066 : ig1=i1-(i1/id1)*n1-1
293 6333066 : ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
294 6333066 : ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
295 6333066 : ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
296 : end if
297 :
298 : end if ! Cut-off
299 : end do ! End loop on i1
300 : end do
301 : end do
302 :
303 :
304 907950866 : barev(:)=zero
305 :
306 : ! MG: This triggers SIGFPE as cryst is not initialized
307 : !a1=Cryst%rprimd(:,1); b1=two_pi*gprimd(:,1)
308 : !a2=Cryst%rprimd(:,2); b2=two_pi*gprimd(:,2)
309 : !a3=Cryst%rprimd(:,3); b3=two_pi*gprimd(:,3)
310 :
311 215534 : select case(TRIM(vcut%mode))
312 : case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB")
313 1296 : if (vcut%mode == "CRYSTAL") then
314 : ! Analytic integration of 4pi/q^2 over the volume element:
315 : ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
316 : ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
317 : ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (see gwa.pdf, appendix A.4)
318 1296 : q0_vol = two_pi**3 / (nkpt_bz * ucvol)
319 1296 : vcut%i_sz = four_pi*7.44*q0_vol**(-two_thirds)
320 :
321 : !TODO FBruneval: cryst is not available here, find a workaround!
322 : !else if (vcut%mode == "AUXILIARY_FUNCTION") then
323 : ! ! Numerical integration of the exact-exchange divergence through the
324 : ! ! auxiliary function of Carrier et al. PRB 75, 205126 (2007) [[cite:Carrier2007]].
325 : ! vcut%i_sz = carrier_isz(cryst, 1, qpoint, rcut, comm)
326 :
327 : !else if (vcut%mode == "AUX_GB") then
328 : ! ! We use the auxiliary function of a Gygi-Baldereschi variant [[cite:Gigy1986]]
329 : ! vcut%i_sz = gygi_baldereschi_isz(cryst, 1, qpoint, vc_ecut, ng, gvec)
330 :
331 : else
332 0 : ABI_ERROR(sjoin("Need treatment of 1/q^2 singularity! for mode", vcut%mode))
333 : end if
334 :
335 5309712 : do ig=1,nfft
336 5309712 : if (abs(gpq(ig))<tol4) then
337 162 : barev(ig) = vcut%i_sz
338 5308254 : else if (gpq(ig)<=cutoff) then
339 1574640 : barev(ig) = gpq2(ig)
340 : end if
341 : end do
342 :
343 :
344 : case('SPHERE') ! Spherical cutoff
345 :
346 : !
347 : ! Treatment of the divergence at q+g=zero
348 : !
349 : ! rcut is not set (rcut<=0), use the default Spencer-Alavi definition:
350 130986 : if ( rcut < tol8 ) then
351 129594 : rcut0= (three*nkpt_bz*ucvol/four_pi)**(one/three)
352 : else
353 : rcut0 = rcut
354 : end if
355 :
356 574623786 : do ig=1,nfft
357 574623786 : if (abs(gpq(ig))<tol4) then
358 8429 : barev(ig) = two_pi*rcut0**two
359 574484371 : else if (gpq(ig)<=cutoff) then
360 206086387 : barev(ig) = gpq2(ig) * (one - cos( rcut0*sqrt(four_pi/gpq2(ig)) ) )
361 : end if
362 : end do
363 :
364 : case('CYLINDER')
365 : !FBruneval: not working. For instance, Cryst is never initialized
366 0 : ABI_BUG("Cylinder cutoff coding is not finalized")
367 :
368 0 : test=COUNT(ABS(vcut%vcutgeo)>tol6)
369 0 : ABI_CHECK(test==1,'Wrong cutgeo for cylinder')
370 :
371 : ! === Beigi method is the default one, i.e infinite cylinder of radius rcut ===
372 : ! * Negative values to use Rozzi method with finite cylinder of extent hcyl.
373 0 : opt_cylinder=1; vcut%hcyl=zero; vcut%pdir(:)=0
374 0 : do ii=1,3
375 0 : check=vcut%vcutgeo(ii)
376 0 : if (ABS(check)>tol6) then
377 0 : vcut%pdir(ii)=1
378 0 : if (check<zero) then ! use Rozzi's method.
379 0 : vcut%hcyl=ABS(check)*NORM2(Cryst%rprimd(:,ii))
380 0 : opt_cylinder=2
381 : end if
382 : end if
383 : end do
384 :
385 0 : test=COUNT(vcut%pdir==1)
386 0 : ABI_CHECK((test==1),'Wrong pdir for cylinder')
387 0 : if (vcut%pdir(3)/=1) then
388 0 : ABI_ERROR("The cylinder must be along the z-axis")
389 : end if
390 :
391 0 : ABI_BUG("cutoff cylinder API has changed!")
392 :
393 : ! call cutoff_cylinder(nfft,gq,ng,Gsph%gvec,vcut%rcut,vcut%hcyl,vcut%pdir,&
394 : !& vcut%boxcenter,Cryst%rprimd,barev,opt_cylinder,comm)
395 :
396 : ! === If Beigi, treat the limit q--> 0 ===
397 0 : if (opt_cylinder==1) then
398 0 : npar=8; npt=100 ; gamma_pt=RESHAPE((/0,0,0/),(/3,1/))
399 0 : ABI_MALLOC(qfit,(3,npt))
400 0 : ABI_MALLOC(vcfit,(1,npt))
401 0 : if (nfft==1) then
402 0 : ABI_ERROR("nfft == 1 not supported when Beigi's method is used")
403 : endif
404 0 : qfit(:,:)=zero
405 0 : step=half/(npt*(nfft-1)) ; qfit(3,:)=arth(tol6,step,npt)
406 :
407 : !call cutoff_cylinder(npt,qfit,1,gamma_pt,vcut%rcut,vcut%hcyl,vcut%pdir,&
408 : ! vcut%boxcenter,Cryst%rprimd,vcfit,opt_cylinder,comm)
409 :
410 0 : ABI_MALLOC(xx,(npt))
411 0 : ABI_MALLOC(yy,(npt))
412 0 : ABI_MALLOC(sigma,(npt))
413 0 : ABI_MALLOC(par,(npar))
414 0 : ABI_MALLOC(var,(npar))
415 0 : ABI_MALLOC(cov,(npar,npar))
416 0 : do ii=1,npt
417 0 : xx(ii)=normv(qfit(:,ii),gmet,'G')
418 : end do
419 0 : ABI_FREE(qfit)
420 0 : sigma=one ; yy(:)=vcfit(1,:)
421 0 : ABI_FREE(vcfit)
422 :
423 0 : bz_plane=l2norm(b1.x.b2)
424 0 : dx=(xx(2)-xx(1))
425 0 : integ=yy(2)*dx*3.0/2.0
426 0 : integ=integ + SUM(yy(3:npt-2))*dx
427 0 : integ=integ+yy(npt-1)*dx*3.0/2.0
428 0 : write(std_out,*)' simple integral',integ
429 0 : q0_volsph=(two_pi)**3/(nkpt_bz*ucvol)
430 0 : q0_vol=bz_plane*two*xx(npt)
431 0 : write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
432 0 : vcut%i_sz=bz_plane*two*integ/q0_vol
433 0 : write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
434 0 : write(std_out,*)' Cylindrical cutoff value ',vcut%i_sz
435 :
436 0 : ABI_FREE(xx)
437 0 : ABI_FREE(yy)
438 0 : ABI_FREE(sigma)
439 0 : ABI_FREE(par)
440 0 : ABI_FREE(var)
441 0 : ABI_FREE(cov)
442 :
443 : else
444 : ! In Rozzi"s method the lim q+G --> 0 is finite.
445 0 : vcut%i_sz=barev(1)
446 : end if
447 :
448 : CASE('SLAB')
449 : !FBruneval: not working. For instance, Cryst is never initialized
450 0 : ABI_BUG("Slab cutoff coding is not finalized")
451 :
452 0 : test=COUNT(vcut%vcutgeo/=zero)
453 0 : ABI_CHECK(test==2,"Wrong vcutgeo")
454 : !
455 : ! Two methods available
456 : !
457 : ! === Default is Beigi"s method ===
458 0 : opt_slab=1; vcut%alpha(:)=zero
459 0 : if (ANY(vcut%vcutgeo<zero)) opt_slab=2
460 0 : vcut%pdir(:)=zero
461 0 : do ii=1,3
462 0 : check=vcut%vcutgeo(ii)
463 0 : if (ABS(check)>zero) then ! Use Rozzi"s method with a finite slab along x-y
464 0 : vcut%pdir(ii)=1
465 0 : if (check<zero) vcut%alpha(ii)=normv(check*Cryst%rprimd(:,ii),rmet,'R')
466 : end if
467 : end do
468 :
469 : ! Beigi"s method: the slab must be along x-y and R must be L_Z/2.
470 0 : if (opt_slab==1) then
471 0 : ABI_CHECK(ALL(vcut%pdir == (/1,1,0/)),"Surface must be in the x-y plane")
472 : !vcut%rcut = half*SQRT(DOT_PRODUCT(a3,a3))
473 : end if
474 :
475 0 : ABI_BUG("cutoff surface API has changed!")
476 : !call cutoff_slab(nfft,gq,ng,Gsph%gvec,gprimd,vcut%rcut,&
477 : ! vcut%boxcenter,vcut%pdir,vcut%alpha,barev,opt_slab)
478 :
479 : !
480 : ! === If Beigi, treat the limit q--> 0 ===
481 0 : if (opt_slab==1) then
482 : ! Integrate numerically in the plane close to 0
483 0 : npt=100 ! Number of points in 1D
484 : gamma_pt=RESHAPE((/0,0,0/),(/3,1/)) ! Gamma point
485 0 : ABI_MALLOC(qfit,(3,npt))
486 0 : ABI_MALLOC(qcart,(3,npt))
487 0 : ABI_MALLOC(vcfit,(1,npt))
488 0 : if (nfft==1) then
489 0 : ABI_ERROR("nfft == 1 not supported when Beigi's method is used")
490 : endif
491 0 : qfit(:,:)=zero
492 0 : qcart(:,:)=zero
493 : ! Size of the third vector
494 0 : bz_plane=l2norm(b3)
495 0 : q0_volsph=(two_pi)**3/(nkpt_bz*ucvol)
496 : ! radius that gives the same volume as q0_volsph
497 : ! Let's assume that c is perpendicular to the plane
498 : ! We also assume isotropic BZ around gamma
499 0 : step=sqrt((q0_volsph/bz_plane)/pi)/npt
500 :
501 : !step=half/(npt*(Qmesh%nibz-1))
502 : ! Let's take qpoints along 1 line, the vcut does depend only on the norm
503 0 : qcart(1,:)=arth(tol6,step,npt)
504 :
505 0 : do ii = 1,npt
506 0 : qfit(:,ii) = MATMUL(TRANSPOSE(Cryst%rprimd),qcart(:,ii))/(2*pi)
507 : end do
508 :
509 0 : ABI_BUG("cutoff surface API has changed!")
510 :
511 : ! call cutoff_slab(npt,qfit,1,gamma_pt,gprimd,vcut%rcut,&
512 : ! vcut%boxcenter,vcut%pdir,vcut%alpha,vcfit,opt_slab)
513 :
514 0 : ABI_MALLOC(xx,(npt))
515 0 : ABI_MALLOC(yy,(npt))
516 0 : ABI_MALLOC(sigma,(npt))
517 0 : do ii=1,npt
518 : !xx(ii)=qfit(1,:)
519 0 : xx(ii)=normv(qfit(:,ii),gmet,'G')
520 : end do
521 0 : ABI_FREE(qfit)
522 0 : sigma=one
523 0 : yy(:)=vcfit(1,:)
524 : !yy(:)=one
525 0 : ABI_FREE(vcfit)
526 0 : dx=(xx(2)-xx(1))
527 : ! integ = \int dr r f(r)
528 0 : integ=xx(2)*yy(2)*dx*3.0/2.0
529 0 : integ=integ + DOT_PRODUCT(xx(3:npt-2),yy(3:npt-2))*dx
530 0 : integ=integ+xx(npt-1)*yy(npt-1)*dx*3.0/2.0
531 0 : write(std_out,*)' simple integral',integ
532 0 : q0_vol=bz_plane*pi*xx(npt)**2
533 0 : write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
534 0 : vcut%i_sz=bz_plane*2*pi*integ/q0_vol
535 0 : write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
536 0 : write(std_out,*)' Cylindrical cutoff value ',vcut%i_sz
537 :
538 0 : ABI_FREE(xx)
539 0 : ABI_FREE(yy)
540 : else
541 : ! In Rozzi"s method the lim q+G --> 0 is finite.
542 0 : vcut%i_sz=barev(1)
543 : end if
544 :
545 : CASE('ERF')
546 :
547 0 : do ig=1,nfft
548 0 : if (abs(gpq(ig))<tol4) then
549 : !FIXME FBruneval check this value, ERFC value was wrong, so why not this one.
550 0 : barev(ig) = zero ! Stupid definition to remember something should be done here.
551 0 : else if (gpq(ig)<=cutoff) then
552 : !FIXME FBruneval shortrange does not make sense here (it is an optional argument that may not be present)
553 : ! and ERF is long range any way
554 0 : if (shortrange) then
555 0 : barev(ig) = + gpq2(ig) * exp( -pi * rcut**2 /gpq2(ig) )
556 : end if
557 : end if
558 : end do
559 :
560 : CASE('ERFC')
561 :
562 328017368 : do ig=1,nfft
563 328017368 : if (abs(gpq(ig))<tol4) then
564 : !FBruneval there was a wrong value here
565 4392 : barev(ig) = pi * rcut**2
566 327931020 : else if (gpq(ig)<=cutoff) then
567 : ! gpq2 is 4 pi / (q+G)**2
568 : ! 4 pi / (q+G)**2 * [ 1 - exp( -1/4 * Rc**2 * (q+G)**2 ) ]
569 89991532 : barev(ig) = gpq2(ig) * ( one - exp( -pi * rcut**2 / gpq2(ig) ) )
570 : end if
571 : end do
572 :
573 : case default
574 0 : write(msg,'(3a)')'No cut-off applied to the Coulomb Potential.', ch10, &
575 0 : 'Either icutcoul value not allowed or not defined.'
576 214238 : ABI_WARNING(msg)
577 : end select
578 :
579 214238 : if (izero==1) then
580 : ! Set contribution of unbalanced components to zero
581 11548 : if (qeq0==1) then !q=0
582 2047 : call zerosym(barev,cplex1,n1,n2,n3)
583 9501 : else if (qeq05==1) then
584 : !q=1/2; this doesn't work in parallel
585 5229 : ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
586 5229 : ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
587 5229 : ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
588 5229 : if (abs(abs(qpoint(1))-half)<tol12) then
589 2988 : if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
590 2988 : if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
591 : end if
592 5229 : if (abs(abs(qpoint(2))-half)<tol12) then
593 2988 : if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
594 2988 : if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
595 : end if
596 5229 : if (abs(abs(qpoint(3))-half)<tol12) then
597 2988 : if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
598 2988 : if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
599 : end if
600 5229 : call zerosym(barev,cplex1,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
601 : end if
602 : end if
603 :
604 214238 : ABI_FREE(gq)
605 214238 : ABI_FREE(gpq)
606 214238 : ABI_FREE(gpq2)
607 :
608 214238 : end subroutine barevcoul
609 : !!***
610 :
611 0 : end module m_barevcoul
612 : !!***
|
