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