Branch data 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 : : !! [i_sz]=Integral of the q --> 0 singularity, used for the q+G=0 component with icutcoul=6 and 7
132 : : !! (auxiliary function of Carrier et al. or Gaussian one, computed by carrier_isz_gmet or gb_isz
133 : : !! in fock_init). Required for icutcoul=6 and 7.
134 : : !!
135 : : !! OUTPUT
136 : : !! barev(nfft)=4pi/(G+q)**2, q+G=0 component is set carefully
137 : : !!
138 : : !! NOTES
139 : : !! This routine operates on the full FFT mesh. DO NOT PASS MPI_TYPE
140 : : !! One can easily implemente MPI-FFT by just calling this routine and then
141 : : !! extracting the G-vectors treated by the node.
142 : : !!
143 : : !! SOURCE
144 : :
145 : 214670 : subroutine barevcoul(rcut,icutcoul,qpoint,gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,izero,barev,shortrange,i_sz)
146 : :
147 : : !Arguments ------------------------------------
148 : : !scalars
149 : : integer,intent(in) :: icutcoul,nfft,nkpt_bz,izero
150 : : real(dp),intent(in) :: rcut,gsqcut,ucvol
151 : : logical,intent(in),optional:: shortrange
152 : : real(dp),intent(in),optional :: i_sz
153 : : !arrays
154 : : integer,intent(in) :: ngfft(18)
155 : : integer :: ng!!!!
156 : : real(dp),intent(in) :: qpoint(3)
157 : : real(dp),intent(inout) :: gmet(3,3)
158 : : real(dp),intent(inout) :: barev(nfft)
159 : : !real(dp) :: a1(3),a2(3),a3(3)
160 : : real(dp) :: b1(3),b2(3),b3(3),rmet(3,3) !,gprimd(3,3),
161 : 214670 : type(MPI_type) :: mpi_enreg !!!!
162 [ + + + + : 11162840 : type(crystal_t) :: Cryst !!!!
+ + + + +
+ + + + +
+ + + + ]
163 : : !type(gsphere_t) :: Gsph
164 : : !Local variables-------------------------------
165 : : !scalars
166 : : integer,parameter :: empty(3,3)=zero
167 : : integer,parameter :: cplex1=1
168 : : integer :: comm
169 : : integer :: ii1,i1,i2,i23,i3,id1,id2,id3,icutcoul_local
170 : : integer :: ig,ig0,ig1,ig2,ig3,ig1min,ig1max,ig2min,ig2max,ig3min,ig3max
171 : : integer :: ii,ing,n1,n2,n3,npar,npt
172 : : integer :: opt_cylinder,opt_slab,test
173 : : integer :: qeq0,qeq05
174 : : real(dp),parameter :: tolfix=1.000000001e0_dp ! Same value as the one used in hartre
175 : : real(dp) :: check,step
176 : : real(dp) :: cutoff,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,rcut0
177 : : real(dp) :: bz_plane,dx,integ,q0_vol,q0_volsph
178 : : character(len=500) :: msg
179 : : !arrays
180 : : integer :: id(3), gamma_pt(3,1)
181 : 214670 : real(dp),allocatable :: gq(:,:),gpq(:),gpq2(:)
182 : 214670 : real(dp),allocatable :: vcfit(:,:),xx(:),yy(:)
183 : 214670 : real(dp),allocatable :: cov(:,:),par(:),qfit(:,:),sigma(:),var(:),qcart(:,:)
184 : 214670 : type(vcut_t) :: vcut !!!!
185 : : !
186 : 214670 : comm=mpi_enreg%comm_world
187 : : !
188 : : ! === Save dimension and other useful quantities in vcut% ===
189 [ + + ]: 858680 : vcut%nfft = PRODUCT(ngfft(1:3)) ! Number of points in the FFT mesh.
190 : : ! ng and gvec are not used yet we don't want to use them without being defined
191 : 214670 : ng = -1
192 : 214670 : vcut%ng = ng ! Number of G-vectors in the Coulomb matrix elements.
193 : 214670 : vcut%rcut = rcut ! Cutoff radius for cylinder.
194 : 214670 : vcut%hcyl = zero ! Length of finite cylinder (Rozzi"s method, default is Beigi).
195 : 214670 : vcut%ucvol = ucvol ! Unit cell volume.
196 : :
197 : : !FBruneval: comment the definitions below since Cryst and dtset have never been initialized!
198 : : !vcut%rprimd = Cryst%rprimd(:,:) ! Dimensional direct lattice.
199 : : !vcut%boxcenter = dtset%boxcenter ! boxcenter at the moment is supposed to be at the origin.
200 : : !vcut%vcutgeo = dtset%vcutgeo(:) ! Info on the orientation and extension of the cutoff region.
201 : : !
202 : : ! === Define geometry and cutoff radius (if used) ===
203 : 214670 : vcut%mode='NONE'
204 : 214670 : icutcoul_local=icutcoul
205 : :
206 : : ! for short-range exchange (e.g. HSE06), enforce ERFC
207 [ + + ]: 214670 : if (PRESENT(shortrange)) then
208 [ - + ]: 81956 : if (shortrange) then
209 : : icutcoul_local=5
210 : : end if
211 : : end if
212 : : ! -------------------------------------
213 : :
214 [ + + ]: 132714 : if (icutcoul_local == 0) vcut%mode = 'SPHERE'
215 [ - + ]: 214670 : if (icutcoul_local == 1) vcut%mode = 'CYLINDER'
216 [ - + ]: 214670 : if (icutcoul_local == 2) vcut%mode = 'SLAB'
217 [ + + ]: 214670 : if (icutcoul_local == 3) vcut%mode = 'CRYSTAL'
218 [ - + ]: 214670 : if (icutcoul_local == 4) vcut%mode = 'ERF'
219 [ + + ]: 214670 : if (icutcoul_local == 5) vcut%mode = 'ERFC'
220 [ + + ]: 214670 : if (icutcoul_local == 6) vcut%mode = 'AUXILIARY_FUNCTION'
221 [ + + ]: 214670 : if (icutcoul_local == 7) vcut%mode = 'AUX_GB'
222 : :
223 : : !Initialize a few quantities
224 : 214670 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
225 : 214670 : cutoff = gsqcut * tolfix
226 [ + + ]: 909720770 : barev=zero
227 : :
228 : : !Some peculiar values of q: q=0 or q on the BZ edge
229 [ + + ]: 214670 : qeq0=0; if (qpoint(1)**2+qpoint(2)**2+qpoint(3)**2<1.d-15) qeq0=1
230 : 201633 : qeq05=0
231 : : if (qeq0==0) then
232 [ + + + + : 201633 : if (abs(abs(qpoint(1))-half)<tol12.or.abs(abs(qpoint(2))-half)<tol12.or. &
+ + ]
233 : 103089 : & abs(abs(qpoint(3))-half)<tol12) qeq05=1
234 : : end if
235 : :
236 : : !In order to speed the routine, precompute the components of g+q
237 : : !Also check if the booked space was large enough...
238 : :
239 [ - + + - : 644010 : ABI_MALLOC(gq,(3,max(n1,n2,n3)))
+ - + - -
+ - + ]
240 [ - + + - : 644010 : ABI_MALLOC(gpq,(nfft))
+ - + - -
+ - + ]
241 [ - + + - : 429340 : ABI_MALLOC(gpq2,(nfft))
- + - + ]
242 : :
243 [ + + ]: 858680 : do ii=1,3
244 : 644010 : id(ii)=ngfft(ii)/2+2
245 [ + + ]: 10771586 : do ing=1,ngfft(ii)
246 : 9912906 : ig=ing-(ing/id(ii))*ngfft(ii)-1
247 : 10556916 : gq(ii,ing)=ig+qpoint(ii)
248 : : end do
249 : : end do
250 : 214670 : ig1max=-1;ig2max=-1;ig3max=-1
251 : 214670 : ig1min=n1;ig2min=n2;ig3min=n3
252 : 214670 : ig0 = 0 ! Index of q+G = 0 (only if q is a reciprocal lattice vector)
253 : :
254 : 214670 : id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
255 : :
256 : : ! Triple loop on each dimension
257 [ + + ]: 3514028 : do i3=1,n3
258 : 3299358 : ig3=i3-(i3/id3)*n3-1
259 : : ! Precompute some products that do not depend on i2 and i1
260 : 3299358 : gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
261 : 3299358 : gqgm23=gq(3,i3)*gmet(2,3)*2
262 : 3299358 : gqgm13=gq(3,i3)*gmet(1,3)*2
263 [ + + ]: 55986182 : do i2=1,n2
264 : 52472154 : ig2=i2-(i2/id2)*n2-1
265 : 52472154 : i23=n1*(i2-1 +(n2)*(i3-1))
266 : 52472154 : gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
267 : 52472154 : gqgm12=gq(2,i2)*gmet(1,2)*2
268 : 52472154 : gqg2p3=gqgm13+gqgm12
269 : :
270 [ + + ]: 961978254 : do i1=1,n1
271 : 909506100 : ii=i1+i23
272 : 909506100 : gpq(ii)= gs2 + gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
273 : : ! q+G = 0 tested in reduced coordinates (no threshold on |q+G|: small q+G occur for large supercells)
274 [ + + + + : 961978254 : if (abs(gq(1,i1))<tol8 .and. abs(gq(2,i2))<tol8 .and. abs(gq(3,i3))<tol8) then
+ + ]
275 : 13037 : ig0 = ii
276 : 13037 : gpq2(ii) = zero
277 : : else
278 : : ! gpq2 contains 4*pi / |q+G|**2
279 : 909493063 : gpq2(ii) = piinv/gpq(ii)
280 : : end if
281 : : end do
282 : :
283 : : !
284 : : ! Next part looks for ig1min,ig1max that are needed by zerosym
285 : : ! Do the test that eliminates the Gamma point outside of the inner loop
286 : 52472154 : ii1=1
287 [ + + - + ]: 52472154 : if (i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0) then
288 : 52472154 : ii1=2
289 : : end if
290 : :
291 : : ! Final inner loop on the first dimension (note the lower limit)
292 [ + + ]: 965264575 : do i1=ii1,n1
293 : 909493063 : gs = gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
294 : :
295 : : !ii=i1+i23
296 : :
297 [ + + ]: 961965217 : if (gs<=cutoff) then
298 : : ! Identify min/max indexes (to cancel unbalanced contributions later)
299 : : ! Count (q+g)-vectors with similar norm
300 [ + + + + ]: 298177439 : if ((qeq05==1).and.(izero==1)) then
301 : 6790878 : ig1=i1-(i1/id1)*n1-1
302 : 6790878 : ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
303 : 6790878 : ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
304 : 6790878 : ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
305 : : end if
306 : :
307 : : end if ! Cut-off
308 : : end do ! End loop on i1
309 : : end do
310 : : end do
311 : :
312 : :
313 [ + + ]: 909720770 : barev(:)=zero
314 : :
315 : : ! MG: This triggers SIGFPE as cryst is not initialized
316 : : !a1=Cryst%rprimd(:,1); b1=two_pi*gprimd(:,1)
317 : : !a2=Cryst%rprimd(:,2); b2=two_pi*gprimd(:,2)
318 : : !a3=Cryst%rprimd(:,3); b3=two_pi*gprimd(:,3)
319 : :
320 : 216398 : select case(TRIM(vcut%mode))
321 : : case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB")
322 [ + + ]: 1728 : if (vcut%mode == "CRYSTAL") then
323 : : ! Analytic integration of 4pi/q^2 over the volume element:
324 : : ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
325 : : ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
326 : : ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (see gwa.pdf, appendix A.4)
327 : 1296 : q0_vol = two_pi**3 / (nkpt_bz * ucvol)
328 : 1296 : vcut%i_sz = four_pi*7.44*q0_vol**(-two_thirds)
329 : :
330 [ + + + - ]: 432 : else if (vcut%mode == "AUXILIARY_FUNCTION" .or. vcut%mode == "AUX_GB") then
331 : : ! Auxiliary function of Carrier et al. [[cite:Carrier2007]] (6) or Gygi-Baldereschi method [[cite:Gygi1986]]
332 : : ! with the Gaussian auxiliary function of [[cite:Massidda1993]] (7).
333 : : ! i_sz depends on the q-mesh (k-points of the Fock operator), it is computed once in fock_init
334 : : ! (carrier_isz_gmet, gb_isz).
335 [ - + ]: 432 : if (.not. present(i_sz)) then
336 : 0 : ABI_BUG(sjoin("icutcoul=6 or 7 (", vcut%mode, ") requires the argument i_sz"))
337 : : end if
338 : 432 : vcut%i_sz = i_sz
339 : :
340 : : else
341 : 0 : ABI_ERROR(sjoin("Need treatment of 1/q^2 singularity! for mode", vcut%mode))
342 : : end if
343 : :
344 [ + + ]: 7079616 : do ig=1,nfft
345 [ + + ]: 7079616 : if (ig == ig0) then
346 : 216 : barev(ig) = vcut%i_sz
347 [ + + ]: 7077672 : else if (gpq(ig)<=cutoff) then
348 : 2099520 : barev(ig) = gpq2(ig)
349 : : end if
350 : : end do
351 : :
352 : : case ('SPHERE') ! Spherical cutoff
353 : :
354 : : !
355 : : ! Treatment of the divergence at q+g=zero
356 : : !
357 : : ! rcut is not set (rcut<=0), use the default Spencer-Alavi definition:
358 [ + + ]: 130986 : if ( rcut < tol8 ) then
359 : 129594 : rcut0= (three*nkpt_bz*ucvol/four_pi)**(one/three)
360 : : else
361 : : rcut0 = rcut
362 : : end if
363 : :
364 [ + + ]: 574623786 : do ig=1,nfft
365 [ + + ]: 574623786 : if (ig == ig0) then
366 : 8429 : barev(ig) = two_pi*rcut0**two
367 [ + + ]: 574484371 : else if (gpq(ig)<=cutoff) then
368 : 206086387 : barev(ig) = gpq2(ig) * (one - cos( rcut0*sqrt(four_pi/gpq2(ig)) ) )
369 : : end if
370 : : end do
371 : :
372 : : case ('CYLINDER')
373 : : !FBruneval: not working. For instance, Cryst is never initialized
374 : 0 : ABI_BUG("Cylinder cutoff coding is not finalized")
375 : :
376 [ # # # # ]: 0 : test=COUNT(ABS(vcut%vcutgeo)>tol6)
377 [ # # ]: 0 : ABI_CHECK(test==1,'Wrong cutgeo for cylinder')
378 : :
379 : : ! === Beigi method is the default one, i.e infinite cylinder of radius rcut ===
380 : : ! * Negative values to use Rozzi method with finite cylinder of extent hcyl.
381 [ # # ]: 0 : opt_cylinder=1; vcut%hcyl=zero; vcut%pdir(:)=0
382 [ # # ]: 0 : do ii=1,3
383 : 0 : check=vcut%vcutgeo(ii)
384 [ # # ]: 0 : if (ABS(check)>tol6) then
385 : 0 : vcut%pdir(ii)=1
386 [ # # ]: 0 : if (check<zero) then ! use Rozzi's method.
387 [ # # # # : 0 : vcut%hcyl=ABS(check)*NORM2(Cryst%rprimd(:,ii))
# # ]
388 : 0 : opt_cylinder=2
389 : : end if
390 : : end if
391 : : end do
392 : :
393 [ # # # # ]: 0 : test=COUNT(vcut%pdir==1)
394 [ # # ]: 0 : ABI_CHECK((test==1),'Wrong pdir for cylinder')
395 [ # # ]: 0 : if (vcut%pdir(3)/=1) then
396 : 0 : ABI_ERROR("The cylinder must be along the z-axis")
397 : : end if
398 : :
399 : 0 : ABI_BUG("cutoff cylinder API has changed!")
400 : :
401 : : ! call cutoff_cylinder(nfft,gq,ng,Gsph%gvec,vcut%rcut,vcut%hcyl,vcut%pdir,&
402 : : !& vcut%boxcenter,Cryst%rprimd,barev,opt_cylinder,comm)
403 : :
404 : : ! === If Beigi, treat the limit q--> 0 ===
405 [ # # ]: 0 : if (opt_cylinder==1) then
406 : 0 : npar=8; npt=100 ; gamma_pt=RESHAPE((/0,0,0/),(/3,1/))
407 [ # # ]: 0 : ABI_MALLOC(qfit,(3,npt))
408 [ # # ]: 0 : ABI_MALLOC(vcfit,(1,npt))
409 [ # # ]: 0 : if (nfft==1) then
410 : 0 : ABI_ERROR("nfft == 1 not supported when Beigi's method is used")
411 : : endif
412 [ # # # # ]: 0 : qfit(:,:)=zero
413 [ # # ]: 0 : step=half/(npt*(nfft-1)) ; qfit(3,:)=arth(tol6,step,npt)
414 : :
415 : : !call cutoff_cylinder(npt,qfit,1,gamma_pt,vcut%rcut,vcut%hcyl,vcut%pdir,&
416 : : ! vcut%boxcenter,Cryst%rprimd,vcfit,opt_cylinder,comm)
417 : :
418 [ # # ]: 0 : ABI_MALLOC(xx,(npt))
419 [ # # ]: 0 : ABI_MALLOC(yy,(npt))
420 [ # # ]: 0 : ABI_MALLOC(sigma,(npt))
421 [ # # ]: 0 : ABI_MALLOC(par,(npar))
422 [ # # ]: 0 : ABI_MALLOC(var,(npar))
423 [ # # ]: 0 : ABI_MALLOC(cov,(npar,npar))
424 [ # # ]: 0 : do ii=1,npt
425 : 0 : xx(ii)=normv(qfit(:,ii),gmet,'G')
426 : : end do
427 [ # # ]: 0 : ABI_FREE(qfit)
428 [ # # # # ]: 0 : sigma=one ; yy(:)=vcfit(1,:)
429 [ # # ]: 0 : ABI_FREE(vcfit)
430 : :
431 : 0 : bz_plane=l2norm(b1.x.b2)
432 : 0 : dx=(xx(2)-xx(1))
433 : 0 : integ=yy(2)*dx*3.0/2.0
434 [ # # ]: 0 : integ=integ + SUM(yy(3:npt-2))*dx
435 : 0 : integ=integ+yy(npt-1)*dx*3.0/2.0
436 : 0 : write(std_out,*)' simple integral',integ
437 : 0 : q0_volsph=(two_pi)**3/(nkpt_bz*ucvol)
438 : 0 : q0_vol=bz_plane*two*xx(npt)
439 : 0 : write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
440 : 0 : vcut%i_sz=bz_plane*two*integ/q0_vol
441 : 0 : write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
442 : 0 : write(std_out,*)' Cylindrical cutoff value ',vcut%i_sz
443 : :
444 : 0 : ABI_FREE(xx)
445 : 0 : ABI_FREE(yy)
446 [ # # ]: 0 : ABI_FREE(sigma)
447 [ # # ]: 0 : ABI_FREE(par)
448 [ # # ]: 0 : ABI_FREE(var)
449 [ # # ]: 0 : ABI_FREE(cov)
450 : :
451 : : else
452 : : ! In Rozzi"s method the lim q+G --> 0 is finite.
453 : 0 : vcut%i_sz=barev(1)
454 : : end if
455 : :
456 : : CASE ('SLAB')
457 : : !FBruneval: not working. For instance, Cryst is never initialized
458 : 0 : ABI_BUG("Slab cutoff coding is not finalized")
459 : :
460 [ # # # # ]: 0 : test=COUNT(vcut%vcutgeo/=zero)
461 [ # # ]: 0 : ABI_CHECK(test==2,"Wrong vcutgeo")
462 : : !
463 : : ! Two methods available
464 : : !
465 : : ! === Default is Beigi"s method ===
466 [ # # ]: 0 : opt_slab=1; vcut%alpha(:)=zero
467 [ # # # # : 0 : if (ANY(vcut%vcutgeo<zero)) opt_slab=2
# # ]
468 [ # # ]: 0 : vcut%pdir(:)=zero
469 [ # # ]: 0 : do ii=1,3
470 : 0 : check=vcut%vcutgeo(ii)
471 [ # # ]: 0 : if (ABS(check)>zero) then ! Use Rozzi"s method with a finite slab along x-y
472 : 0 : vcut%pdir(ii)=1
473 [ # # # # ]: 0 : if (check<zero) vcut%alpha(ii)=normv(check*Cryst%rprimd(:,ii),rmet,'R')
474 : : end if
475 : : end do
476 : :
477 : : ! Beigi"s method: the slab must be along x-y and R must be L_Z/2.
478 [ # # ]: 0 : if (opt_slab==1) then
479 [ # # # # : 0 : ABI_CHECK(ALL(vcut%pdir == (/1,1,0/)),"Surface must be in the x-y plane")
# # ]
480 : : !vcut%rcut = half*SQRT(DOT_PRODUCT(a3,a3))
481 : : end if
482 : :
483 : 0 : ABI_BUG("cutoff surface API has changed!")
484 : : !call cutoff_slab(nfft,gq,ng,Gsph%gvec,gprimd,vcut%rcut,&
485 : : ! vcut%boxcenter,vcut%pdir,vcut%alpha,barev,opt_slab)
486 : :
487 : : !
488 : : ! === If Beigi, treat the limit q--> 0 ===
489 [ # # ]: 0 : if (opt_slab==1) then
490 : : ! Integrate numerically in the plane close to 0
491 : 0 : npt=100 ! Number of points in 1D
492 : : gamma_pt=RESHAPE((/0,0,0/),(/3,1/)) ! Gamma point
493 [ # # ]: 0 : ABI_MALLOC(qfit,(3,npt))
494 [ # # ]: 0 : ABI_MALLOC(qcart,(3,npt))
495 [ # # ]: 0 : ABI_MALLOC(vcfit,(1,npt))
496 [ # # ]: 0 : if (nfft==1) then
497 : 0 : ABI_ERROR("nfft == 1 not supported when Beigi's method is used")
498 : : endif
499 [ # # # # ]: 0 : qfit(:,:)=zero
500 [ # # # # ]: 0 : qcart(:,:)=zero
501 : : ! Size of the third vector
502 : 0 : bz_plane=l2norm(b3)
503 : 0 : q0_volsph=(two_pi)**3/(nkpt_bz*ucvol)
504 : : ! radius that gives the same volume as q0_volsph
505 : : ! Let's assume that c is perpendicular to the plane
506 : : ! We also assume isotropic BZ around gamma
507 : 0 : step=sqrt((q0_volsph/bz_plane)/pi)/npt
508 : :
509 : : !step=half/(npt*(Qmesh%nibz-1))
510 : : ! Let's take qpoints along 1 line, the vcut does depend only on the norm
511 [ # # ]: 0 : qcart(1,:)=arth(tol6,step,npt)
512 : :
513 [ # # ]: 0 : do ii = 1,npt
514 [ # # ]: 0 : qfit(:,ii) = MATMUL(TRANSPOSE(Cryst%rprimd),qcart(:,ii))/(2*pi)
515 : : end do
516 : :
517 : 0 : ABI_BUG("cutoff surface API has changed!")
518 : :
519 : : ! call cutoff_slab(npt,qfit,1,gamma_pt,gprimd,vcut%rcut,&
520 : : ! vcut%boxcenter,vcut%pdir,vcut%alpha,vcfit,opt_slab)
521 : :
522 [ # # ]: 0 : ABI_MALLOC(xx,(npt))
523 [ # # ]: 0 : ABI_MALLOC(yy,(npt))
524 [ # # ]: 0 : ABI_MALLOC(sigma,(npt))
525 [ # # ]: 0 : do ii=1,npt
526 : : !xx(ii)=qfit(1,:)
527 : 0 : xx(ii)=normv(qfit(:,ii),gmet,'G')
528 : : end do
529 [ # # ]: 0 : ABI_FREE(qfit)
530 [ # # ]: 0 : sigma=one
531 [ # # ]: 0 : yy(:)=vcfit(1,:)
532 : : !yy(:)=one
533 [ # # ]: 0 : ABI_FREE(vcfit)
534 : 0 : dx=(xx(2)-xx(1))
535 : : ! integ = \int dr r f(r)
536 : 0 : integ=xx(2)*yy(2)*dx*3.0/2.0
537 [ # # ]: 0 : integ=integ + DOT_PRODUCT(xx(3:npt-2),yy(3:npt-2))*dx
538 : 0 : integ=integ+xx(npt-1)*yy(npt-1)*dx*3.0/2.0
539 : 0 : write(std_out,*)' simple integral',integ
540 : 0 : q0_vol=bz_plane*pi*xx(npt)**2
541 : 0 : write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
542 : 0 : vcut%i_sz=bz_plane*2*pi*integ/q0_vol
543 : 0 : write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
544 : 0 : write(std_out,*)' Cylindrical cutoff value ',vcut%i_sz
545 : :
546 : 0 : ABI_FREE(xx)
547 : 0 : ABI_FREE(yy)
548 : : else
549 : : ! In Rozzi"s method the lim q+G --> 0 is finite.
550 : 0 : vcut%i_sz=barev(1)
551 : : end if
552 : :
553 : : CASE ('ERF')
554 : :
555 [ # # ]: 0 : do ig=1,nfft
556 [ # # ]: 0 : if (abs(gpq(ig))<tol4) then
557 : : !FIXME FBruneval check this value, ERFC value was wrong, so why not this one.
558 : 0 : barev(ig) = zero ! Stupid definition to remember something should be done here.
559 [ # # ]: 0 : else if (gpq(ig)<=cutoff) then
560 : : !FIXME FBruneval shortrange does not make sense here (it is an optional argument that may not be present)
561 : : ! and ERF is long range any way
562 [ # # ]: 0 : if (shortrange) then
563 : 0 : barev(ig) = + gpq2(ig) * exp( -pi * rcut**2 /gpq2(ig) )
564 : : end if
565 : : end if
566 : : end do
567 : :
568 : : CASE ('ERFC')
569 : :
570 [ + + ]: 328017368 : do ig=1,nfft
571 [ + + ]: 328017368 : if (ig == ig0) then
572 : : !FBruneval there was a wrong value here
573 : 4392 : barev(ig) = pi * rcut**2
574 [ + + ]: 327931020 : else if (gpq(ig)<=cutoff) then
575 : : ! gpq2 is 4 pi / (q+G)**2
576 : : ! 4 pi / (q+G)**2 * [ 1 - exp( -1/4 * Rc**2 * (q+G)**2 ) ]
577 : 89991532 : barev(ig) = gpq2(ig) * ( one - exp( -pi * rcut**2 / gpq2(ig) ) )
578 : : end if
579 : : end do
580 : :
581 : : case default
582 : 0 : write(msg,'(3a)')'No cut-off applied to the Coulomb Potential.', ch10, &
583 : 0 : 'Either icutcoul value not allowed or not defined.'
584 [ + + - - : 214670 : ABI_WARNING(msg)
- + - ]
585 : : end select
586 : :
587 [ + + ]: 214670 : if (izero==1) then
588 : : ! Set contribution of unbalanced components to zero
589 [ + + ]: 11980 : if (qeq0==1) then !q=0
590 : 2101 : call zerosym(barev,cplex1,n1,n2,n3)
591 [ + + ]: 9879 : else if (qeq05==1) then
592 : : !q=1/2; this doesn't work in parallel
593 [ + - ]: 5607 : ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
594 [ + - ]: 5607 : ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
595 [ + - ]: 5607 : ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
596 [ + + ]: 5607 : if (abs(abs(qpoint(1))-half)<tol12) then
597 [ + + ]: 3204 : if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
598 [ + + ]: 3204 : if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
599 : : end if
600 [ + + ]: 5607 : if (abs(abs(qpoint(2))-half)<tol12) then
601 [ + + ]: 3204 : if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
602 [ + + ]: 3204 : if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
603 : : end if
604 [ + + ]: 5607 : if (abs(abs(qpoint(3))-half)<tol12) then
605 [ + - ]: 3204 : if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
606 [ - + ]: 3204 : if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
607 : : end if
608 : 5607 : call zerosym(barev,cplex1,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
609 : : end if
610 : : end if
611 : :
612 [ - + ]: 214670 : ABI_FREE(gq)
613 [ - + ]: 214670 : ABI_FREE(gpq)
614 [ - + ]: 214670 : ABI_FREE(gpq2)
615 : :
616 [ - + - + : 214670 : end subroutine barevcoul
- + - + ]
617 : : !!***
618 : :
619 [ # # # # : 0 : end module m_barevcoul
# # # # #
# # # #
# ]
620 : : !!***
|
