Branch data Line data Source code
1 : : !!****m* ABINIT/m_vcoul
2 : : !! NAME
3 : : !! m_vcoul
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the definition of the vcoul_t as well
7 : : !! as procedures to calculate the Coulomb interaction in reciprocal space
8 : : !! taking into account a possible cutoff in real space.
9 : : !! Procedures to deal with the singularity for q --> 0 are also provided.
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 1999-2026 ABINIT group (MG, FB)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : module m_vcoul
27 : :
28 : : use defs_basis
29 : : use m_abicore
30 : : use m_errors
31 : : use m_xmpi
32 : : use m_splines
33 : : use m_sort
34 : :
35 : : use m_fstrings, only : sjoin, itoa
36 : : use m_special_funcs, only : abi_derf
37 : : use m_bessel, only : calck0
38 : : use m_io_tools, only : open_file
39 : : use m_gwdefs, only : GW_TOLQ0
40 : : use m_numeric_tools, only : arth, geop, imin_loc, llsfit_svd, l2norm, OPERATOR(.x.), quadrature, isdiagmat
41 : : use m_hide_lapack, only : matrginv
42 : : use m_geometry, only : normv, metric
43 : : use m_qplusg, only : cmod_qpg
44 : : use m_crystal, only : crystal_t
45 : : use m_bz_mesh, only : kmesh_t
46 : : use m_gsphere, only : gsphere_t
47 : : use m_fftcore, only : get_kg
48 : : use m_dtfil, only : isfile
49 : :
50 : : ! Cut-off methods modules
51 : : use m_cutoff_sphere, only : cutoff_sphere
52 : : use m_cutoff_slab, only : cutoff_slab
53 : : use m_cutoff_cylinder, only : cutoff_cylinder
54 : :
55 : : implicit none
56 : :
57 : : public :: gw_icutcoul_to_mode
58 : : public :: carrier_isz
59 : :
60 : : private
61 : : !!***
62 : :
63 : : !!****t* m_vcoul/vcoul_t
64 : : !! NAME
65 : : !! vcoul_t
66 : : !!
67 : : !! FUNCTION
68 : : !! This data type contains the square root of the Fourier components of the Coulomb interaction
69 : : !! calculated taking into account a possible cutoff. It also stores info on the particular geometry
70 : : !! used for the cutoff as well as quantities required to deal with the Coulomb divergence for q --> 0.
71 : : !!
72 : : !! SOURCE
73 : :
74 : : type,public :: vcoul_t
75 : :
76 : : integer :: ng = -1
77 : : ! Number of G-vectors
78 : :
79 : : integer :: nqibz = -1
80 : : ! Number of irreducible q-points
81 : :
82 : : integer :: nqlwl = -1
83 : : ! Number of small q-points around Gamma
84 : :
85 : : real(dp) :: alpha(3) = -one
86 : : ! Length of the finite slab
87 : :
88 : : real(dp) :: rcut = -one
89 : : ! Cutoff radius.
90 : :
91 : : real(dp) :: i_sz = huge(one)
92 : : ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2
93 : :
94 : : real(dp) :: i_sz_resid = huge(one)
95 : : ! Residual difference between the i_sz in the sigma self-energy for exchange,
96 : : ! and the i_sz already present in the generalized Kohn-Sham eigenenergies
97 : : ! Initialized to the same value as i_sz
98 : :
99 : : real(dp) :: hcyl = -one
100 : : ! Length of the finite cylinder along the periodic dimension
101 : :
102 : : real(dp) :: ucvol = -one
103 : : ! Volume of the unit cell
104 : :
105 : : character(len=50) :: mode
106 : : ! String defining the cutoff mode, possible values are: sphere,cylinder,slab,crystal
107 : :
108 : : integer :: pdir(3)
109 : : ! 1 if the system is periodic along this direction
110 : :
111 : : real(dp) :: boxcenter(3) = -1
112 : : ! 1 if the point in inside the cutoff region 0 otherwise
113 : : ! Reduced coordinates of the center of the box (input variable)
114 : :
115 : : real(dp) :: vcutgeo(3) = huge(one)
116 : : ! For each reduced direction gives the length of the finite system
117 : : ! 0 if the system is infinite along this direction.
118 : : ! negative values indicate that a finite size has to be used.
119 : :
120 : : real(dp) :: rprimd(3,3) = zero
121 : : ! Lattice vectors in real space.
122 : :
123 : : real(dp),allocatable :: qibz(:,:)
124 : : ! (3, nqibz)
125 : : ! q-points in the IBZ.
126 : :
127 : : real(dp),allocatable :: qlwl(:,:)
128 : : ! (3, nqlwl)
129 : : ! q-points for the treatment of the Coulomb singularity.
130 : :
131 : : complex(gwp),allocatable :: vc_sqrt(:,:)
132 : : ! (ng, nqibz)
133 : : ! Square root of the Coulomb interaction in reciprocal space.
134 : : ! complex-valued to allow for a possible cutoff (Rozzi's method)
135 : :
136 : : complex(gwp),allocatable :: vcqlwl_sqrt(:,:)
137 : : ! (ng, nqlwl)
138 : : ! Square root of the Coulomb term calculated for small q-points
139 : :
140 : : complex(gwp),allocatable :: vc_sqrt_resid(:,:)
141 : : ! (ng, nqibz)
142 : : ! Square root of the residual difference between the Coulomb interaction in the sigma self-energy for exchange,
143 : : ! and the Coulomb interaction already present in the generalized Kohn-Sham eigenenergies (when they come from an hybrid)
144 : : ! Given in reciprocal space. At the call to vcoul_init, it is simply initialized at the value of vc_sqrt(:,:),
145 : : ! and only later modified. A cutoff might be applied.
146 : :
147 : : contains
148 : : procedure :: init => vcoul_init ! Main creation method.
149 : : procedure :: plot => vcoul_plot ! Plot vc in real and reciprocal space.
150 : : procedure :: print => vcoul_print ! Print info on the object.
151 : : procedure :: free => vcoul_free ! Free memory
152 : : end type vcoul_t
153 : : !!***
154 : :
155 : : !!****t* m_vcoul/mc_t
156 : : !! NAME
157 : : !! mc_t
158 : : !!
159 : : !! FUNCTION
160 : : !! Mimicking the BerkeleyGW technique
161 : : !! A Monte-Carlo sampling of each miniBZ surrounding each (q+G) point.
162 : : !! However:
163 : : !! - extended to multiple shifts
164 : : !! - with an adaptative number of MonteCarlo sampling points
165 : : !!
166 : : !! SOURCE
167 : :
168 : : type, public :: mc_t
169 : :
170 : : integer :: nmc_max = -1
171 : :
172 : : real(dp) :: q0sph = -one
173 : :
174 : : real(dp) :: ucvol = -one
175 : :
176 : : real(dp) :: gmet(3,3) = -one
177 : :
178 : : real(dp),allocatable :: qran(:,:)
179 : : ! (3, nmc_max)
180 : :
181 : : contains
182 : : procedure :: init => mc_init
183 : : procedure :: integrate => mc_integrate
184 : : procedure :: free => mc_free
185 : : end type mc_t
186 : : !!***
187 : :
188 : : !!****t* m_vcoul/vcgen_t
189 : : !! NAME
190 : : !! vcgen_t
191 : : !!
192 : : !! FUNCTION
193 : : !!
194 : : !! SOURCE
195 : :
196 : : type, public :: vcgen_t
197 : :
198 : : integer :: nkbz = -1
199 : : ! Number of k-points in full BZ.
200 : :
201 : : integer :: opt_cylinder
202 : :
203 : : integer :: opt_slab
204 : :
205 : : real(dp) :: alpha(3) = -one
206 : : ! Lenght of the finite surface.
207 : :
208 : : real(dp) :: rcut = -one
209 : : ! Cutoff radius.
210 : :
211 : : real(dp) :: hcyl = -one
212 : : ! Length of the finite cylinder along the periodic dimension
213 : :
214 : : character(len=50) :: mode
215 : : ! String defining the cutoff mode
216 : :
217 : : integer :: pdir(3)
218 : : ! 1 if the system is periodic along this direction
219 : :
220 : : real(dp) :: boxcenter(3) = -1
221 : : ! 1 if the point in inside the cutoff region 0 otherwise
222 : : ! Reduced coordinates of the center of the box (input variable)
223 : :
224 : : real(dp) :: vcutgeo(3) = huge(one)
225 : : ! For each reduced direction gives the length of the finite system
226 : : ! 0 if the system is infinite along this direction.
227 : : ! negative values indicate that a finite size has to be used.
228 : :
229 : : real(dp) :: i_sz = huge(one)
230 : : ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2
231 : :
232 : : type(mc_t) :: mc
233 : : ! Monte carlo integrator.
234 : :
235 : : contains
236 : : procedure :: init => vcgen_init ! Initialize the object
237 : : procedure :: get_vc_sqrt => vcgen_get_vc_sqrt ! Compute sqrt(vc(q,g))
238 : : procedure :: free => vcgen_free ! Free memory
239 : : procedure :: print => vcgen_print ! Print info on the object.
240 : : end type vcgen_t
241 : : !!***
242 : :
243 : : ! private stuff
244 : : real(dp),parameter :: TOLQ0 = 1.d-3
245 : :
246 : : CONTAINS !========================================================================================
247 : : !!***
248 : :
249 : : !!****f* m_vcoul/gw_icutcoul_to_mode
250 : : !! NAME
251 : : !! gw_icutcoul_to_mode
252 : : !!
253 : : !! FUNCTION
254 : : !! Convert gw_icutcoul_to_mode to mode string.
255 : : !!
256 : : !! SOURCE
257 : :
258 : 404 : subroutine gw_icutcoul_to_mode(gw_icutcoul, mode)
259 : :
260 : : !Arguments ------------------------------------
261 : : integer,intent(in) :: gw_icutcoul
262 : : character(len=*),intent(out) :: mode
263 : : ! *************************************************************************
264 : :
265 [ + - ]: 404 : mode = 'NONE'
266 [ + + + - : 404 : if (gw_icutcoul == 0) mode = 'SPHERE'
+ - ]
267 [ + + + - : 404 : if (gw_icutcoul == 1) mode = 'CYLINDER'
+ - ]
268 [ - + - - : 404 : if (gw_icutcoul == 2) mode = 'SLAB'
- - ]
269 [ + + + - : 404 : if (gw_icutcoul == 3) mode = 'CRYSTAL'
+ - ]
270 [ + + + - : 404 : if (gw_icutcoul == 4) mode = 'ERF'
+ - ]
271 [ + + + - : 404 : if (gw_icutcoul == 5) mode = 'ERFC'
+ - ]
272 [ + + + - : 404 : if (gw_icutcoul == 6) mode = 'AUXILIARY_FUNCTION'
+ - ]
273 [ + + + - : 404 : if (gw_icutcoul == 7) mode = 'AUX_GB'
+ - ]
274 [ + + + - : 404 : if (gw_icutcoul == 14) mode = 'MINIBZ-ERF'
+ - ]
275 [ + + + - : 404 : if (gw_icutcoul == 15) mode = 'MINIBZ-ERFC'
+ - ]
276 [ + + + - : 404 : if (gw_icutcoul == 16) mode = 'MINIBZ'
+ - ]
277 : :
278 [ + - ]: 404 : end subroutine gw_icutcoul_to_mode
279 : : !!***
280 : :
281 : : !!****f* m_vcoul/vcoul_init
282 : : !! NAME
283 : : !! vcoul_init
284 : : !!
285 : : !! FUNCTION
286 : : !! Perform general check and initialize the data type containing information on the cutoff technique
287 : : !! Note %vc_sqrt_resid and %i_sz_resid are simply initialized at the same value as %vc_sqrt and %i_sz
288 : : !!
289 : : !! INPUTS
290 : : !! Gsph=Info of the G sphere.
291 : : !! Qmesh=Info on the q-point sampling.
292 : : !! Kmesh=Info on the k-point sampling.
293 : : !! rcut=Cutoff radius for the cylinder.
294 : : !! gw_icutcoul=Option of the cutoff technique.
295 : : !! vcutgeo(3)= Info on the orientation and extension of the cutoff region.
296 : : !! ng=Number of G-vectors to be used to describe the Coulomb interaction
297 : : !! nqlwl=Number of point around Gamma for treatment of long-wavelength limit
298 : : !! qlwl(3,nqlwl)= The nqlwl "small" q-points
299 : : !! comm=MPI communicator.
300 : : !!
301 : : !! SOURCE
302 : :
303 : 404 : subroutine vcoul_init(vcp, Gsph, Cryst, Qmesh, Kmesh, rcut, gw_icutcoul, vcutgeo, vc_ecut, ng, nqlwl, qlwl, comm)
304 : :
305 : : !Arguments ------------------------------------
306 : : !scalars
307 : : class(vcoul_t),intent(out) :: vcp
308 : : integer,intent(in) :: ng,nqlwl, gw_icutcoul, comm
309 : : real(dp),intent(in) :: rcut, vc_ecut
310 : : type(kmesh_t),target,intent(in) :: Kmesh, Qmesh
311 : : type(gsphere_t),target,intent(in) :: Gsph
312 : : type(crystal_t),intent(in) :: Cryst
313 : : !arrays
314 : : real(dp),intent(in) :: qlwl(3,nqlwl),vcutgeo(3)
315 : :
316 : : !Local variables-------------------------------
317 : : !scalars
318 : : integer,parameter :: master=0
319 : : integer :: nqibz, nqbz, nkbz, iqlwl, iq_ibz, opt_cylinder,my_rank,nprocs, opt_slab
320 : : real(dp) :: bz_geometry_factor,q0_vol, rcut2
321 : : character(len=500) :: msg
322 [ + + + + ]: 5252 : type(mc_t) :: mc
323 : : !arrays
324 : : integer :: units(2)
325 : 404 : integer, contiguous, pointer :: gvec(:,:)
326 : : real(dp) :: a1(3),a2(3),a3(3),b1(3),b2(3),b3(3)
327 : 404 : real(dp),allocatable :: vcoul(:,:),vcoul_lwl(:,:)
328 : 404 : real(dp),contiguous, pointer :: qibz(:,:), qbz(:,:)
329 : : ! *************************************************************************
330 : :
331 : : !call wrtout(std_out, "in vcoul_init")
332 : 404 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
333 [ + + ]: 1212 : units = [std_out, ab_out]
334 : :
335 : 404 : nqibz = qmesh%nibz; nqbz = qmesh%nbz
336 : 404 : qibz => qmesh%ibz; qbz => qmesh%bz
337 : 404 : nkbz = kmesh%nbz
338 : :
339 : : ! Save dimension and other useful quantities in vcp
340 : 404 : vcp%ng = ng ! Number of G-vectors in the Coulomb matrix elements.
341 : 404 : vcp%nqibz = nqibz ! Number of irred q-point.
342 : 404 : vcp%nqlwl = nqlwl ! Number of small q-directions to deal with singularity and non Analytic behavior.
343 : 404 : vcp%rcut = rcut ! Cutoff radius for cylinder.
344 : 404 : vcp%hcyl = zero ! Length of finite cylinder (Rozzi"s method, default is Beigi).
345 : 404 : vcp%ucvol = cryst%ucvol ! Unit cell volume.
346 [ + + + + ]: 5252 : vcp%rprimd = Cryst%rprimd(:,:) ! Dimensional direct lattice.
347 [ + + ]: 1616 : vcp%boxcenter = zero ! Boxcenter at the moment is supposed to be at the origin.
348 [ + + ]: 1616 : vcp%vcutgeo = vcutgeo(:) ! Info on the orientation and extension of the cutoff region.
349 : :
350 : 404 : gvec => Gsph%gvec
351 : :
352 : 404 : call gw_icutcoul_to_mode(gw_icutcoul, vcp%mode)
353 : :
354 [ - + + - : 1212 : ABI_MALLOC(vcp%qibz, (3, nqibz))
+ - + - -
+ - + -
+ ]
355 [ + - + - : 11072 : vcp%qibz = Qmesh%ibz(:,:)
- + - - +
+ + + ]
356 [ - + + - : 1212 : ABI_MALLOC(vcp%qlwl, (3, nqlwl))
+ - + - -
+ - + -
+ ]
357 [ - + + + : 2424 : vcp%qlwl = qlwl(:,:)
+ + ]
358 : :
359 : : ! ===============================================
360 : : ! == Calculation of the FT of the Coulomb term ==
361 : : ! ===============================================
362 : : a1 = cryst%rprimd(:,1); a2 = cryst%rprimd(:,2); a3 = cryst%rprimd(:,3)
363 : : b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
364 : :
365 [ - + + - : 1616 : ABI_MALLOC(vcoul , (ng, nqibz))
- + + - +
- + - - +
- + ]
366 [ - + + - : 1616 : ABI_MALLOC(vcoul_lwl, (ng, nqlwl))
- + + - +
- + - - +
- + ]
367 : :
368 : 407 : select case (trim(vcp%mode))
369 : : case ('MINIBZ', 'MINIBZ-ERFC', 'MINIBZ-ERF')
370 : 3 : call mc%init(cryst%rprimd, cryst%ucvol, cryst%gprimd, cryst%gmet, kmesh%kptrlatt)
371 : :
372 : 3 : rcut2 = vcp%rcut**2
373 [ + + ]: 12 : do iq_ibz=1,nqibz
374 : 12 : call mc%integrate(vcp%mode, qibz(:, iq_ibz), ng, gvec, rcut2, nkbz, vcoul(:, iq_ibz), comm)
375 : : end do
376 : :
377 : : ! Treat the limit q --> 0
378 : 3 : vcp%i_sz = vcoul(1, 1)
379 : :
380 [ + + ]: 6 : do iqlwl=1,nqlwl
381 : 6 : call mc%integrate(vcp%mode, qlwl(:, iqlwl), ng, gvec, rcut2, nkbz, vcoul_lwl(:, iqlwl), comm)
382 : : end do
383 : :
384 : 3 : call mc%free()
385 : :
386 : : case ('SPHERE')
387 : : ! A non-positive value of rcut activates the recipe of Spencer & Alavi, PRB 77, 193110 (2008) [[cite:Spencer2008]].
388 [ + + ]: 46 : if (vcp%rcut < tol12) then
389 : 10 : vcp%rcut = (cryst%ucvol * nkbz * 3.d0 / four_pi) ** third
390 : 10 : write(msg,'(2a,2x,f8.4,a)')ch10,' Using calculated rcut: ',vcp%rcut,' to have same volume as the BvK crystal'
391 : 10 : call wrtout(std_out, msg)
392 : : end if
393 [ + + ]: 184 : vcp%vcutgeo = zero
394 : :
395 [ + + ]: 137 : do iq_ibz=1,nqibz
396 : 137 : call cutoff_sphere(qibz(:,iq_ibz), ng, gvec, cryst%gmet, vcp%rcut, vcoul(:,iq_ibz))
397 : : end do
398 : :
399 : : ! q-points for optical limit.
400 [ + + ]: 92 : do iqlwl=1,nqlwl
401 : 92 : call cutoff_sphere(qlwl(:,iqlwl), ng, gvec, cryst%gmet, vcp%rcut, vcoul_lwl(:,iqlwl))
402 : : end do
403 : :
404 : : ! Treat the limit q --> 0
405 : : ! The small cube is approximated by a sphere, while vc(q=0) = 2piR**2.
406 : : ! if a single q-point is used, the expression for the volume is exact.
407 : 46 : vcp%i_sz = two_pi * vcp%rcut**2
408 : :
409 : : case ('CYLINDER')
410 : 2 : call cylinder_setup(cryst, vcp%vcutgeo, vcp%hcyl, vcp%pdir, opt_cylinder)
411 : :
412 [ + + ]: 12 : do iq_ibz=1,nqibz
413 : : call cutoff_cylinder(qibz(:,iq_ibz), ng, gvec, vcp%rcut, vcp%hcyl, vcp%pdir,&
414 : 12 : vcp%boxcenter, Cryst%rprimd, vcoul(:,iq_ibz), opt_cylinder, comm)
415 : : end do
416 : :
417 : : ! q-points for optical limit.
418 [ + + ]: 4 : do iqlwl=1,nqlwl
419 : : call cutoff_cylinder(qlwl(:,iqlwl), ng, gvec, vcp%rcut, vcp%hcyl, vcp%pdir,&
420 : 4 : vcp%boxcenter, Cryst%rprimd, vcoul_lwl(:,iqlwl), opt_cylinder, comm)
421 : : end do
422 : :
423 : : ! If Beigi, treat the limit q --> 0.
424 [ + + ]: 2 : if (opt_cylinder == 1) then
425 : 1 : call beigi_cylinder_limit(opt_cylinder, cryst, nqibz, nkbz, vcp%rcut, vcp%hcyl, vcp%boxcenter, vcp%pdir, vcp%i_sz)
426 : : else
427 : : ! In Rozzi's method the lim q+G --> 0 is finite.
428 : 1 : vcp%i_sz = vcoul(1,1)
429 : : end if
430 : :
431 : : case ('SLAB')
432 : 0 : call surface_setup(cryst, vcp%vcutgeo, vcp%alpha, vcp%rcut, vcp%pdir, opt_slab)
433 : :
434 [ # # ]: 0 : do iq_ibz=1,nqibz
435 : : call cutoff_slab(qibz(:,iq_ibz), ng, gvec, cryst%gprimd, vcp%rcut, &
436 : 0 : vcp%boxcenter, vcp%pdir, vcp%alpha, vcoul(:,iq_ibz), opt_slab)
437 : : end do
438 : :
439 : : ! q-points for optical limit.
440 [ # # ]: 0 : do iqlwl=1,nqlwl
441 : : call cutoff_slab(qlwl(:,iq_ibz), ng, gvec, cryst%gprimd, vcp%rcut, &
442 : 0 : vcp%boxcenter, vcp%pdir, vcp%alpha, vcoul_lwl(:,iqlwl), opt_slab)
443 : : end do
444 : :
445 : : ! If Beigi, treat the limit q --> 0.
446 [ # # ]: 0 : if (opt_slab == 1) then
447 : : ! Integrate numerically in the plane close to 0
448 : : call beigi_surface_limit(opt_slab, cryst, nqibz, nkbz, vcp%rcut, vcp%alpha, &
449 : 0 : vcp%boxcenter, vcp%pdir, vcp%i_sz)
450 : : else
451 : : ! In Rozzi's method the lim q+G --> 0 is finite.
452 : 0 : vcp%i_sz=vcoul(1,1)
453 : : end if
454 : :
455 : : case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB")
456 [ + + ]: 2662 : do iq_ibz=1,nqibz
457 : 2328 : call cmod_qpg(nqibz, iq_ibz, qibz, ng, gvec, cryst%gprimd, vcoul(:,iq_ibz))
458 : :
459 [ + + ]: 2662 : if (iq_ibz == 1) then
460 : : ! The singularity is treated using vcoul_lwl.
461 : 334 : vcoul(1, iq_ibz) = zero
462 [ + + ]: 45452 : vcoul(2:,iq_ibz) = four_pi / vcoul(2:,iq_ibz)**2
463 : : else
464 [ + + ]: 219536 : vcoul(:,iq_ibz) = four_pi / vcoul(:,iq_ibz)**2
465 : : end if
466 : : end do ! iq_ibz
467 : :
468 : : ! q-points for optical limit.
469 [ + + ]: 668 : do iqlwl=1,nqlwl
470 : 668 : call cmod_qpg(nqlwl, iqlwl, qlwl, ng, gvec, cryst%gprimd, vcoul_lwl(:,iqlwl))
471 : : end do
472 [ + + + + ]: 46120 : vcoul_lwl = four_pi/vcoul_lwl**2
473 : :
474 : : ! Treatment of 1/q^2 singularity
475 : :
476 [ + + ]: 334 : if (vcp%mode == "CRYSTAL") then
477 : : ! Analytic integration of 4pi/q^2 over the volume element:
478 : : ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
479 : : ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
480 : : ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (see gwa.pdf, appendix A.4)
481 : 117 : q0_vol = (two_pi) **3 / (nkbz*cryst%ucvol); bz_geometry_factor=zero
482 : 117 : vcp%i_sz = four_pi*7.44*q0_vol**(-two_thirds)
483 : :
484 [ + + ]: 217 : else if (vcp%mode == "AUXILIARY_FUNCTION") then
485 : : ! Numerical integration of the exact-exchange divergence through the
486 : : ! auxiliary function of Carrier et al. PRB 75, 205126 (2007) [[cite:Carrier2007]].
487 : 209 : vcp%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, comm)
488 : :
489 [ + - ]: 8 : else if (vcp%mode == "AUX_GB") then
490 : : ! We use the auxiliary function of a Gygi-Baldereschi variant [[cite:Gigy1986]]
491 : 8 : vcp%i_sz = gygi_baldereschi_isz(cryst, nqbz, qbz, vc_ecut, ng, gvec)
492 : :
493 : : else
494 : 0 : ABI_ERROR(sjoin("Need treatment of 1/q^2 singularity! for mode", vcp%mode))
495 : : end if
496 : :
497 : : case ('ERF')
498 : : ! Modified long-range only Coulomb interaction thanks to the error function:
499 : : ! * Vc = erf(r/rcut)/r
500 : : ! * The singularity is treated using vcoul_lwl.
501 [ + + ]: 4 : do iq_ibz=1,nqibz
502 : 3 : call cmod_qpg(nqibz, iq_ibz, qibz, ng, gvec, cryst%gprimd, vcoul(:,iq_ibz))
503 : :
504 : : ! The Fourier transform of the error function reads
505 [ + + ]: 4 : if (iq_ibz == 1) then
506 : 1 : vcoul(1, iq_ibz) = zero
507 [ + + ]: 137 : vcoul(2:,iq_ibz) = four_pi/(vcoul(2:,iq_ibz)**2) * EXP( -0.25d0 * (vcp%rcut*vcoul(2:,iq_ibz))**2 )
508 : : else
509 [ + + ]: 276 : vcoul(:,iq_ibz) = four_pi/(vcoul(:, iq_ibz)**2) * EXP( -0.25d0 * (vcp%rcut*vcoul(: ,iq_ibz))**2 )
510 : : end if
511 : : end do
512 : :
513 : : ! q-points for optical limit.
514 [ + + ]: 2 : do iqlwl=1,nqlwl
515 : 2 : call cmod_qpg(nqlwl, iqlwl, qlwl, ng, gvec, cryst%gprimd, vcoul_lwl(:,iqlwl))
516 : : end do
517 [ + + + + ]: 139 : vcoul_lwl = four_pi/(vcoul_lwl**2) * EXP( -0.25d0 * (vcp%rcut*vcoul_lwl)**2 )
518 : :
519 : : ! === Treat 1/q^2 singularity ===
520 : : ! * We use the auxiliary function from PRB 75, 205126 (2007) [[cite:Carrier2007]]
521 : 1 : vcp%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, comm)
522 : :
523 : : case ('ERFC')
524 : : ! * Use a modified short-range only Coulomb interaction thanks to the complementary error function:
525 : : ! $ V_c = [1-erf(r/r_{cut})]/r $
526 : : ! * The Fourier transform of the error function reads
527 : : ! vcoul=four_pi/(vcoul**2) * ( 1.d0 - exp( -0.25d0 * (vcp%rcut*vcoul)**2 ) )
528 [ + + ]: 143 : do iq_ibz=1,nqibz
529 : 125 : call cmod_qpg(nqibz, iq_ibz, qibz,ng, gvec, cryst%gprimd, vcoul(:,iq_ibz))
530 : :
531 [ + + ]: 143 : if (iq_ibz == 1) then
532 : 18 : vcoul(1 ,iq_ibz) = zero
533 [ + + ]: 1602 : vcoul(2:,iq_ibz) = four_pi/(vcoul(2:,iq_ibz)**2) * ( one - EXP( -0.25d0 * (vcp%rcut*vcoul(2:,iq_ibz))**2 ) )
534 : : else
535 [ + + ]: 10446 : vcoul(:, iq_ibz) = four_pi/(vcoul(:, iq_ibz)**2) * ( one - EXP( -0.25d0 * (vcp%rcut*vcoul(:, iq_ibz))**2 ) )
536 : : end if
537 : : end do ! iq_ibz
538 : :
539 : : ! q-points for optical limit.
540 [ + + ]: 36 : do iqlwl=1,nqlwl
541 : 36 : call cmod_qpg(nqlwl, iqlwl, qlwl, ng, gvec, cryst%gprimd, vcoul_lwl(:,iqlwl))
542 : : end do
543 [ + + + + ]: 1638 : vcoul_lwl = four_pi/(vcoul_lwl**2) * ( one - EXP( -0.25d0 * (vcp%rcut*vcoul_lwl)**2 ) )
544 : :
545 : : ! === Treat 1/q^2 singularity ===
546 : : ! * There is NO singularity in this case.
547 : 18 : vcp%i_sz = pi * vcp%rcut**2 ! Final result stored here
548 : :
549 : : case default
550 [ + + + - : 404 : ABI_BUG(sjoin('Unsupported cutoff mode:', vcp%mode))
+ + + - ]
551 : : end select
552 : :
553 : : !call wrtout(std_out, sjoin("vcp%i_sz", ftoa(vcp%i_sz)))
554 : 404 : vcp%i_sz_resid = vcp%i_sz
555 : :
556 : : ! Store final results in complex array as Rozzi's cutoff can give real negative values
557 [ - + + - : 1616 : ABI_MALLOC(vcp%vc_sqrt, (ng, nqibz))
- + + - +
- + - - +
- + - + ]
558 [ - + + - : 1212 : ABI_MALLOC(vcp%vc_sqrt_resid, (ng, nqibz))
- + + - -
+ - + -
+ ]
559 [ - + + + : 356380 : vcp%vc_sqrt = CMPLX(vcoul, zero)
+ + ]
560 [ + + + + ]: 355976 : vcp%vc_sqrt = SQRT(vcp%vc_sqrt)
561 [ + - + - : 356380 : vcp%vc_sqrt_resid = vcp%vc_sqrt
- + - - +
+ + + ]
562 [ - + ]: 404 : ABI_FREE(vcoul)
563 : :
564 [ - + + - : 1616 : ABI_MALLOC(vcp%vcqlwl_sqrt, (ng, nqlwl))
- + + - +
- + - - +
- + - + ]
565 [ - + + + : 120488 : vcp%vcqlwl_sqrt = CMPLX(vcoul_lwl, zero)
+ + ]
566 [ + + + + ]: 120084 : vcp%vcqlwl_sqrt = SQRT(vcp%vcqlwl_sqrt)
567 [ - + ]: 404 : ABI_FREE(vcoul_lwl)
568 : :
569 : 404 : call vcp%print(units)
570 : : !call wrtout(std_out, " exiting vcoul_init")
571 : :
572 [ + - - + ]: 1212 : end subroutine vcoul_init
573 : : !!***
574 : :
575 : :
576 : : !!****f* m_vcoul/cylinder_setup
577 : : !! NAME
578 : : !! cylinder_setup
579 : : !!
580 : : !! FUNCTION
581 : : !!
582 : : !! SOURCE
583 : :
584 : 2 : subroutine cylinder_setup(cryst, vcutgeo, hcyl, pdir, opt_cylinder)
585 : :
586 : : type(crystal_t),intent(in) :: cryst
587 : : real(dp),intent(in) :: vcutgeo(3)
588 : : real(dp),intent(out) :: hcyl
589 : : integer,intent(out) :: pdir(3), opt_cylinder
590 : :
591 : : !Local variables-------------------------------
592 : : integer :: ii
593 : : real(dp),parameter :: tol999 = 999.0
594 : : real(dp) :: check
595 : : ! *************************************************************************
596 : :
597 [ + + + + : 8 : ABI_CHECK(count(abs(vcutgeo) > tol6) == 1, 'Wrong cutgeo for cylinder')
- + ]
598 : :
599 : : ! Beigi's method is the default one, i.e infinite cylinder of radius rcut.
600 : : ! Use negative values to use Rozzi's method with finite cylinder of extent hcyl.
601 : 2 : opt_cylinder = 1; hcyl = zero; pdir(:) = 0
602 [ + + ]: 8 : do ii=1,3
603 : 6 : check = vcutgeo(ii)
604 [ + + ]: 8 : if (abs(check) > tol6) then
605 : 2 : pdir(ii) = 1
606 [ + + ]: 2 : if (check < zero) then
607 : : ! use Rozzi's method.
608 [ + + + + : 4 : hcyl = ABS(check) * NORM2(cryst%rprimd(:,ii))
+ - ]
609 : 1 : opt_cylinder = 2
610 : : ! Check to enter the infinite Rozzi treatment
611 [ - + ]: 1 : if(vcutgeo(3) <= -tol999) hcyl = tol12
612 : : end if
613 : : end if
614 : : end do
615 : :
616 [ + + + + : 8 : ABI_CHECK((count(pdir == 1) == 1), 'Wrong pdir for cylinder')
- + ]
617 [ - + ]: 2 : if (pdir(3) /= 1) then
618 : 0 : ABI_ERROR("The cylinder must be along the z-axis")
619 : : end if
620 : :
621 : 2 : end subroutine cylinder_setup
622 : : !!***
623 : :
624 : : !!****f* m_vcoul/surface_setup
625 : : !! NAME
626 : : !! surface_setup
627 : : !!
628 : : !! FUNCTION
629 : : !!
630 : : !! SOURCE
631 : :
632 : 0 : subroutine surface_setup(cryst, vcutgeo, alpha, rcut, pdir, opt_slab)
633 : :
634 : : type(crystal_t),intent(in) :: cryst
635 : : real(dp),intent(in) :: vcutgeo(3)
636 : : real(dp),intent(out) :: alpha(3)
637 : : real(dp),intent(inout) :: rcut
638 : : integer,intent(out) :: pdir(3), opt_slab
639 : :
640 : : !Local variables-------------------------------
641 : : integer :: ii
642 : : real(dp) :: check
643 : : character(len=500) :: msg
644 : : ! *************************************************************************
645 : :
646 [ # # # # : 0 : ABI_CHECK(count(vcutgeo /= zero) == 2, "Wrong vcutgeo")
# # ]
647 : :
648 : : ! Default is Beigi's method.
649 [ # # # # : 0 : opt_slab = 1; if (any(vcutgeo < zero)) opt_slab = 2
# # ]
650 : 0 : pdir(:) = zero; alpha(:)=zero
651 [ # # ]: 0 : do ii=1,3
652 : 0 : check = vcutgeo(ii)
653 [ # # ]: 0 : if (abs(check) > zero) then
654 : : ! Use Rozzi's method with a finite surface along x-y
655 : 0 : pdir(ii) = 1
656 [ # # # # ]: 0 : if (check < zero) alpha(ii) = normv(check * cryst%rprimd(:,ii), cryst%rmet, 'R')
657 : : end if
658 : : end do
659 : :
660 : : ! In Beigi's method, the surface must be along x-y and R must be L_Z/2.
661 [ # # ]: 0 : if (opt_slab == 1) then
662 : 0 : msg = "2D Beigi method, the periodicity must be in the x-y plane. Modify vcutgeo and/or your geometry."
663 [ # # # # : 0 : ABI_CHECK(all(pdir == [1, 1, 0]), msg)
# # ]
664 [ # # ]: 0 : rcut = half*SQRT(DOT_PRODUCT(cryst%rprimd(:,3), cryst%rprimd(:,3)))
665 : : end if
666 : :
667 : 0 : end subroutine surface_setup
668 : : !!***
669 : :
670 : : !!****f* m_vcoul/integratefaux
671 : : !! NAME
672 : : !! integratefaux
673 : : !!
674 : : !! FUNCTION
675 : : !!
676 : : !! SOURCE
677 : :
678 : 420 : real(dp) function integratefaux(rcut, gprimd, ucvol, comm)
679 : :
680 : : real(dp),intent(in) :: rcut, gprimd(3,3), ucvol
681 : : integer,intent(in) :: comm
682 : :
683 : : !Local variables-------------------------------
684 : : integer,parameter :: nref = 3, nq = 50
685 : : integer :: ierr,iq,iqx1,iqy1,iqz1,iqx2,iqy2,iqz2,miniqy1,maxiqy1,nqhalf, nprocs, my_rank
686 : : real(dp) :: invnq,invnq3,qq,weightq,weightxy,weightxyz
687 : : real(dp) :: qq1(3),qq2(3),bb4sinpiqq_2(3,nq),sin2piqq(nq),bb4sinpiqq2_2(3,0:nq),sin2piqq2(3,0:nq)
688 : : real(dp) :: b1(3), b2(3), b3(3), bb(3)
689 : : real(dp) :: b1b1,b2b2,b3b3,b1b2,b2b3,b3b1
690 : : ! *************************************************************************
691 : :
692 : : ! nq is the number of sampling points along each axis for the numerical integration
693 : : ! nref is the area where the mesh is refined
694 : :
695 : 210 : integratefaux = zero
696 : 210 : invnq = one/DBLE(nq)
697 : 210 : invnq3 = invnq**3
698 : 210 : nqhalf = nq/2
699 : :
700 [ + + + + : 2100 : b1 = two_pi * gprimd(:,1); b2 = two_pi * gprimd(:,2); b3 = two_pi * gprimd(:,3)
+ + ]
701 [ + + + + : 2100 : b1b1 = dot_product(b1, b1); b2b2 = dot_product(b2, b2); b3b3 = dot_product(b3, b3)
+ + ]
702 : 210 : bb(1) = b1b1; bb(2) = b2b2; bb(3) = b3b3
703 : : b1b2 = dot_product(b1, b2); b2b3 = dot_product(b2, b3); b3b1 = dot_product(b3, b1)
704 : :
705 : 210 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
706 : :
707 : : ! In order to speed up the calculation, precompute the sines
708 [ + + ]: 10710 : do iq=1,nq
709 : 10500 : qq=DBLE(iq)*invnq-half
710 [ + + ]: 42210 : bb4sinpiqq_2(:,iq)=bb(:)*four*SIN(pi*qq)**2 ; sin2piqq(iq)=SIN(two_pi*qq)
711 : : end do
712 : :
713 [ + + ]: 10710 : do iqx1=1,nq
714 [ + + ]: 10500 : if (modulo(iqx1, nprocs) /= my_rank) cycle ! MPI parallelism
715 : 9300 : qq1(1)=DBLE(iqx1)*invnq-half
716 : : ! Here take advantage of the q <=> -q symmetry:
717 : : ! arrange the sampling of qx, qy space to avoid duplicating calculations. Need weights to do this ...
718 : : !do iqy1=1,nq
719 : 9300 : miniqy1 = nqhalf + 1; maxiqy1 = nq
720 [ + + ]: 9300 : if (iqx1 >= nqhalf) miniqy1 = nqhalf
721 [ + + ]: 9300 : if (iqx1 > nqhalf .and. iqx1 < nq) maxiqy1 = nq - 1
722 : :
723 [ + + ]: 242382 : do iqy1=miniqy1,maxiqy1
724 : 232872 : qq1(2) = DBLE(iqy1)*invnq - half
725 : : ! By default, the factor of two is for the q <=> -q symmetry
726 : 232872 : weightq = invnq3*two
727 : : ! But not all qx qy lines have a symmetric one ...
728 [ + + + + ]: 232872 : if( (iqx1 == nqhalf .or. iqx1 == nq) .and. (iqy1 == nqhalf .or. iqy1 == nq)) weightq = weightq*half
729 : :
730 [ + + ]: 11886972 : do iqz1=1,nq
731 : 11643600 : qq1(3) = DBLE(iqz1)*invnq - half
732 : :
733 : : ! Refine the mesh for the point close to the origin
734 [ + + + + : 11876472 : if( abs(iqx1-nqhalf) <= nref .and. abs(iqy1-nqhalf) <= nref .and. abs(iqz1-nqhalf) <= nref ) then
+ + ]
735 : : ! Note that the set of point is symmetric around the central point, while weights are taken into account
736 [ + + ]: 1692600 : do iq=0,nq
737 [ + + ]: 6640200 : qq2(:) = qq1(:)+ (DBLE(iq)*invnq-half)*invnq
738 [ + + + + ]: 11652900 : bb4sinpiqq2_2(:,iq) =bb(:)*four*SIN(pi*qq2(:))**2; sin2piqq2(:,iq)=SIN(two_pi*qq2(:))
739 : : end do
740 [ + + ]: 1692600 : do iqx2=0,nq
741 : 1660050 : qq2(1)=qq1(1) + (DBLE(iqx2)*invnq-half ) *invnq
742 [ + + ]: 86355150 : do iqy2=0,nq
743 : 84662550 : qq2(2)=qq1(2) + (DBLE(iqy2)*invnq-half ) *invnq
744 : 84662550 : weightxy=invnq3*weightq
745 [ + + ]: 84662550 : if (iqx2 == 0 .or. iqx2 == nq) weightxy = weightxy*half
746 [ + + ]: 84662550 : if (iqy2 == 0 .or. iqy2 == nq) weightxy = weightxy*half
747 [ + + ]: 4404112650 : do iqz2=0,nq
748 : 4317790050 : qq2(3) = qq1(3) + (DBLE(iqz2)*invnq - half) * invnq
749 : 4317790050 : weightxyz = weightxy
750 [ + + ]: 4317790050 : if (iqz2 == 0 .or. iqz2 == nq) weightxyz = weightxy*half
751 : : !
752 : : ! Treat the remaining divergence in the origin as if it would be a spherical integration of 1/q^2
753 : : if (iqx1/=nqhalf .or. iqy1/=nqhalf .or. iqz1/=nqhalf .or. &
754 [ + + + + : 4402452600 : iqx2/=nqhalf .or. iqy2/=nqhalf .or. iqz2/=nqhalf ) then
+ + ]
755 : : !integratefaux=integratefaux+ faux(qq2, rcut, b1, b2, b3) *invnq**6
756 : : integratefaux = integratefaux + &
757 : : faux_fast(qq2, bb4sinpiqq2_2(1,iqx2), bb4sinpiqq2_2(2,iqy2), bb4sinpiqq2_2(3,iqz2), &
758 : 4317789864 : sin2piqq2(1,iqx2), sin2piqq2(2,iqy2), sin2piqq2(3,iqz2), rcut, b1, b2, b3) * weightxyz
759 : : else
760 : 186 : integratefaux = integratefaux + 7.7955* ((two_pi)**3/ucvol*invnq3*invnq3 )**(-2./3.) *invnq3*invnq3
761 : : end if
762 : : end do
763 : : end do
764 : : end do
765 : : else
766 : : ! integratefaux=integratefaux+faux(qq1, rcut, b1, b2, b3)*invnq**3
767 : : integratefaux = integratefaux + &
768 : : faux_fast(qq1, bb4sinpiqq_2(1,iqx1), bb4sinpiqq_2(2,iqy1), bb4sinpiqq_2(3,iqz1), &
769 : 11611050 : sin2piqq(iqx1), sin2piqq(iqy1), sin2piqq(iqz1), rcut, b1, b2, b3) * weightq
770 : : end if
771 : : end do
772 : : end do
773 : : end do
774 : :
775 : 210 : call xmpi_sum(integratefaux, comm, ierr)
776 : :
777 : 210 : end function integratefaux
778 : : !!***
779 : :
780 : 9357 : real(dp) pure function faux(qq, rcut, b1, b2, b3)
781 : :
782 : : !Arguments ------------------------------------
783 : : real(dp),intent(in) :: qq(3)
784 : : real(dp),intent(in) :: rcut
785 : : real(dp),intent(in) :: b1(3), b2(3), b3(3)
786 : :
787 : : !Local variables-------------------------------
788 : : real(dp) :: bb4sinpiqq1_2, bb4sinpiqq2_2, bb4sinpiqq3_2, sin2piqq1, sin2piqq2, sin2piqq3, b1b1,b2b2,b3b3
789 : : ! *************************************************************************
790 : :
791 [ + + + + : 93570 : b1b1 = dot_product(b1, b1); b2b2 = dot_product(b2, b2); b3b3 = dot_product(b3, b3)
+ + ]
792 : :
793 : 9357 : bb4sinpiqq1_2 = b1b1 * four * SIN(pi*qq(1))**2
794 : 9357 : bb4sinpiqq2_2 = b2b2 * four * SIN(pi*qq(2))**2
795 : 9357 : bb4sinpiqq3_2 = b3b3 * four * SIN(pi*qq(3))**2
796 : 9357 : sin2piqq1 = SIN(two_pi*qq(1))
797 : 9357 : sin2piqq2 = SIN(two_pi*qq(2))
798 : 9357 : sin2piqq3 = SIN(two_pi*qq(3))
799 : :
800 : 9357 : faux = faux_fast(qq, bb4sinpiqq1_2, bb4sinpiqq2_2, bb4sinpiqq3_2, sin2piqq1, sin2piqq2, sin2piqq3, rcut, b1, b2, b3)
801 : :
802 : 9357 : end function faux
803 : :
804 : 4329410271 : real(dp) pure function faux_fast(qq, bb4sinpiqq1_2, bb4sinpiqq2_2, bb4sinpiqq3_2, sin2piqq1, sin2piqq2, sin2piqq3, &
805 : : rcut, b1, b2, b3)
806 : :
807 : : !Arguments ------------------------------------
808 : : real(dp),intent(in) :: qq(3)
809 : : real(dp),intent(in) :: bb4sinpiqq1_2, bb4sinpiqq2_2, bb4sinpiqq3_2, sin2piqq1, sin2piqq2, sin2piqq3, rcut
810 : : real(dp),intent(in) :: b1(3), b2(3), b3(3)
811 : :
812 : : !Local variables-------------------------------
813 : : real(dp) :: b1b2,b2b3,b3b1
814 : : ! *************************************************************************
815 : :
816 [ + + + + : 43294102710 : b1b2 = dot_product(b1, b2); b2b3 = dot_product(b2, b3); b3b1 = dot_product(b3, b1)
+ + ]
817 : :
818 : : faux_fast = bb4sinpiqq1_2 + bb4sinpiqq2_2 + bb4sinpiqq3_2 &
819 : : +two*( b1b2 * sin2piqq1*sin2piqq2 &
820 : : +b2b3 * sin2piqq2*sin2piqq3 &
821 : : +b3b1 * sin2piqq3*sin2piqq1 &
822 : 4329410271 : )
823 : :
824 [ + + ]: 4329410271 : if (rcut > tol6) then
825 [ + + ]: 93105424 : faux_fast = two_pi*two_pi/faux_fast * exp( -0.25d0*rcut**2* sum( ( qq(1)*b1(:)+qq(2)*b2(:)+qq(3)*b3(:) )**2 ) )
826 : : else
827 : 4306133915 : faux_fast = two_pi*two_pi/faux_fast
828 : : endif
829 : :
830 : 4329410271 : end function faux_fast
831 : : !!***
832 : :
833 : : !----------------------------------------------------------------------
834 : :
835 : : !!****f* m_vcoul/adapt_nmc
836 : : !! NAME
837 : : !! adapt_nmc
838 : : !!
839 : : !! FUNCTION
840 : : !! Empirical law to decrease the Monte Carlo sampling for large |q+G|, for which the accuracy is not an issue
841 : :
842 : 3393 : integer pure function adapt_nmc(nmc_max, qpg2) result(nmc)
843 : :
844 : : !Arguments ------------------------------------
845 : : integer,intent(in) :: nmc_max
846 : : real(dp),intent(in) :: qpg2
847 : : ! *************************************************************************
848 : :
849 : 3393 : nmc = NINT( nmc_max / ( 1.0_dp + 1.0_dp * qpg2**6 ) )
850 : 3393 : nmc = MIN(nmc_max, nmc)
851 : 3393 : nmc = MAX(1, nmc)
852 : :
853 : 3393 : end function adapt_nmc
854 : : !!***
855 : :
856 : : !----------------------------------------------------------------------
857 : :
858 : : !!****f* m_vcoul/vcoul_plot
859 : : !! NAME
860 : : !! vcoul_plot
861 : : !!
862 : : !! FUNCTION
863 : : !! Plot vccut(q,G) as a function of |q+G|. Calculate also vc in real space.
864 : : !!
865 : : !! INPUTS
866 : : !!
867 : : !! OUTPUT
868 : : !!
869 : : !! SOURCE
870 : :
871 : 0 : subroutine vcoul_plot(Vcp, Qmesh, Gsph, ng, vc, comm)
872 : :
873 : : !Arguments ------------------------------------
874 : : !scalars
875 : : class(vcoul_t),intent(in) :: Vcp
876 : : integer,intent(in) :: ng, comm
877 : : type(kmesh_t),intent(in) :: Qmesh
878 : : type(gsphere_t),intent(in) :: Gsph
879 : : !arrays
880 : : real(dp),intent(in) :: vc(ng,Qmesh%nibz)
881 : :
882 : : !Local variables-------------------------------
883 : : !scalars
884 : : integer,parameter :: master = 0
885 : : integer :: icount,idx_Sm1G,ierr,ig,igs,ii,iq_bz,iq_ibz,iqg,ir,isym,itim
886 : : integer :: my_start,my_stop,nqbz,nqibz,nr,ntasks,my_rank,unt
887 : : real(dp) :: arg,fact,l1,l2,l3,lmax,step,tmp,vcft,vc_bare
888 : : character(len=500) :: msg
889 : : character(len=fnlen) :: filnam
890 : : !arrays
891 : 0 : integer,allocatable :: insort(:)
892 : : real(dp) :: b1(3),b2(3),b3(3),gmet(3,3),gprimd(3,3),qbz(3),qpgc(3)
893 : 0 : real(dp),allocatable :: qpg_mod(:),rr(:,:,:),vcr(:,:),vcr_cut(:,:)
894 : : !************************************************************************
895 : :
896 [ # # ]: 0 : if (trim(Vcp%mode) /= 'CYLINDER') RETURN
897 : :
898 : 0 : my_rank = xmpi_comm_rank(comm)
899 : :
900 : 0 : nqibz=Vcp%nqibz; nqbz=Qmesh%nbz
901 [ # # # # : 0 : gmet=Gsph%gmet; gprimd=Gsph%gprimd
# # # # ]
902 : :
903 [ # # ]: 0 : b1(:)=two_pi*gprimd(:,1)
904 [ # # ]: 0 : b2(:)=two_pi*gprimd(:,2)
905 [ # # ]: 0 : b3(:)=two_pi*gprimd(:,3)
906 : :
907 : : ! Compare in Fourier space the true Coulomb with the cutted one.
908 [ # # ]: 0 : if (my_rank == master) then
909 [ # # # # : 0 : ABI_MALLOC(insort, (nqibz * ng))
# # # # #
# # # ]
910 [ # # # # : 0 : ABI_MALLOC(qpg_mod, (nqibz * ng))
# # # # #
# # # ]
911 : 0 : iqg = 1
912 [ # # ]: 0 : do iq_ibz=1,nqibz
913 [ # # ]: 0 : do ig=1,ng
914 [ # # # # ]: 0 : qpg_mod(iqg) = normv(Qmesh%ibz(:,iq_ibz) + Gsph%gvec(:,ig), gmet,'g')
915 : 0 : insort(iqg) = iqg; iqg = iqg + 1
916 : : end do
917 : : end do
918 : 0 : call sort_dp(nqibz * ng, qpg_mod, insort, tol14)
919 : :
920 : 0 : filnam='_VCoulFT_'
921 : 0 : call isfile(filnam, 'new')
922 [ # # ]: 0 : if (open_file(filnam, msg, newunit=unt, status='new', form='formatted') /= 0) then
923 : 0 : ABI_ERROR(msg)
924 : : end if
925 : : write(unt,'(a,i3,a,i6,a)')&
926 : 0 : '# |q+G| q-point (Tot no.',nqibz,') Gvec (',ng,') vc_bare(q,G) vc_cutoff(q,G) '
927 : :
928 [ # # ]: 0 : do iqg=1,nqibz*ng
929 : 0 : iq_ibz = (insort(iqg) - 1) / ng + 1
930 : 0 : ig = (insort(iqg)) - (iq_ibz-1) * ng
931 : 0 : vc_bare = zero
932 [ # # ]: 0 : if (qpg_mod(iqg) > tol16) vc_bare = four_pi / qpg_mod(iqg) ** 2
933 : : write(unt,'(f12.6,2x,3f8.4,2x,3i6,2x,2es14.6)')&
934 : 0 : qpg_mod(iqg), Qmesh%ibz(:,iq_ibz), Gsph%gvec(:,ig), vc_bare, vc(ig, iq_ibz)
935 : : end do
936 : :
937 : 0 : close(unt)
938 [ # # ]: 0 : ABI_FREE(insort)
939 [ # # ]: 0 : ABI_FREE(qpg_mod)
940 : : end if ! my_rank==master
941 : :
942 : : ! Fourier transform back to real space just to check cutoff implementation.
943 : 0 : ntasks= nqbz * ng
944 : 0 : call xmpi_split_work(ntasks, comm, my_start, my_stop)
945 : :
946 [ # # # # : 0 : l1 = NORM2(Vcp%rprimd(:,1))
# # ]
947 [ # # # # : 0 : l2 = NORM2(Vcp%rprimd(:,2))
# # ]
948 [ # # # # : 0 : l3 = NORM2(Vcp%rprimd(:,3))
# # ]
949 : :
950 : 0 : nr = 50
951 : 0 : lmax=MAX(l1,l2,l3) ; step=lmax/(nr-1)
952 : 0 : fact = one / (Vcp%ucvol * nqbz)
953 : :
954 : : ! numb coding
955 [ # # # # : 0 : ABI_CALLOC(rr, (3, nr, 3))
# # # # ]
956 [ # # ]: 0 : do ir=1,nr
957 [ # # # # ]: 0 : rr(1:3,ir,1:3)=(ir-1)*step
958 : : end do
959 : :
960 [ # # # # : 0 : ABI_CALLOC(vcr, (nr, 3))
# # ]
961 [ # # # # : 0 : ABI_CALLOC(vcr_cut, (nr, 3))
# # ]
962 : :
963 [ # # ]: 0 : do iq_bz=1,nqbz
964 : 0 : call Qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym, itim)
965 [ # # ]: 0 : if (ABS(qbz(1))<0.01) qbz(1)=zero
966 [ # # ]: 0 : if (ABS(qbz(2))<0.01) qbz(2)=zero
967 [ # # ]: 0 : if (ABS(qbz(3))<0.01) qbz(3)=zero
968 [ # # # # : 0 : igs=1; if (ALL(qbz(:)==zero)) igs=2
# # ]
969 [ # # ]: 0 : do ig=igs,ng
970 : 0 : icount=ig+(iq_bz-1)*ng
971 [ # # # # ]: 0 : if (icount < my_start .or. icount > my_stop) CYCLE
972 : 0 : idx_Sm1G = Gsph%rottbm1(ig,itim,isym) ! IS{^-1}G
973 : 0 : vcft=vc(idx_Sm1G,iq_ibz)
974 [ # # # # ]: 0 : qpgc(:)=qbz(:)+Gsph%gvec(:,ig) ; qpgc(:)=b1(:)*qpgc(1)+b2(:)*qpgc(2)+b3(:)*qpgc(3)
975 [ # # # # : 0 : tmp=NORM2(qpgc) ; tmp=tmp**2
# # ]
976 [ # # ]: 0 : do ii=1,3
977 [ # # ]: 0 : do ir=1,nr
978 [ # # ]: 0 : arg=DOT_PRODUCT(rr(:,ir,ii),qpgc)
979 : 0 : vcr_cut(ir,ii)=vcr_cut(ir,ii) + vcft*COS(arg)
980 : 0 : vcr (ir,ii)=vcr (ir,ii) + four_pi/tmp*COS(arg)
981 : : end do
982 : : end do
983 : : end do !ig
984 : : end do !iq_ibz
985 : :
986 : 0 : call xmpi_sum_master(vcr_cut,master,comm,ierr)
987 : 0 : call xmpi_sum_master(vcr ,master,comm,ierr)
988 : :
989 [ # # ]: 0 : if (my_rank == master) then
990 : 0 : filnam='_VCoulR_'
991 : 0 : call isfile(filnam, 'new')
992 [ # # ]: 0 : if (open_file(filnam,msg,newunit=unt,status='new',form='formatted') /= 0) then
993 : 0 : ABI_ERROR(msg)
994 : : end if
995 : 0 : write(unt,'(a)')'# length vc_bare(r) vc_cut(r) '
996 [ # # ]: 0 : do ir=1,nr
997 [ # # # # ]: 0 : write(unt,'(7es18.6)')(ir-1)*step,(fact*vcr(ir,ii),fact*vcr_cut(ir,ii),ii=1,3)
998 : : end do
999 : 0 : close(unt)
1000 : : end if
1001 : :
1002 [ # # ]: 0 : ABI_FREE(rr)
1003 [ # # ]: 0 : ABI_FREE(vcr)
1004 [ # # ]: 0 : ABI_FREE(vcr_cut)
1005 : :
1006 : 0 : end subroutine vcoul_plot
1007 : : !!***
1008 : :
1009 : : !----------------------------------------------------------------------
1010 : :
1011 : : !!****f* m_vcoul/vcoul_print
1012 : : !! NAME
1013 : : !! vcoul_print
1014 : : !!
1015 : : !! FUNCTION
1016 : : !! Print info on a Coulomb datatype.
1017 : : !!
1018 : : !! INPUTS
1019 : : !! units=Unit numbers for output
1020 : : !! [prtvol]=Verbosity level
1021 : : !!
1022 : : !! SOURCE
1023 : :
1024 [ - + ]: 407 : subroutine vcoul_print(Vcp, units, prtvol)
1025 : :
1026 : : !Arguments ------------------------------------
1027 : : !scalars
1028 : : class(vcoul_t),intent(in) :: Vcp
1029 : : integer,intent(in) :: units(:)
1030 : : integer,intent(in),optional :: prtvol
1031 : :
1032 : : !Local variables-------------------------------
1033 : : !scalars
1034 : : integer :: ii, my_prtvol, iqlwl
1035 : : character(len=500) :: msg
1036 : : ! *************************************************************************
1037 : :
1038 : 407 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
1039 : :
1040 : 408 : select case (Vcp%mode)
1041 : :
1042 : : case ('MINIBZ')
1043 : 1 : write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1044 : 1 : call wrtout(units, msg)
1045 : :
1046 : : case ('MINIBZ-ERF')
1047 : 1 : write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1048 : 1 : call wrtout(units, msg)
1049 : 1 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1050 : 1 : ' === Error function cutoff === ',ch10,ch10,&
1051 : 2 : ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1052 : 1 : call wrtout(units, msg)
1053 : :
1054 : : case ('MINIBZ-ERFC')
1055 : 1 : write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1056 : 1 : call wrtout(units, msg)
1057 : 1 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1058 : 1 : ' === Complement Error function cutoff === ',ch10,ch10,&
1059 : 2 : ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1060 : 1 : call wrtout(units, msg)
1061 : :
1062 : : case ('SPHERE')
1063 : 46 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1064 : 46 : ' === Spherical cutoff === ',ch10,ch10,&
1065 : 46 : ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10,&
1066 : 92 : ' Volume of the sphere .. ',four_pi/three*Vcp%rcut**3,' [Bohr^3] '
1067 : : !FB: This has no meaning here! & ' Sphere centered at .... ',Vcp%boxcenter,' (r.l.u) ',ch10
1068 : : !MG It might be useful if the system is not centered on the origin because in this case the
1069 : : ! matrix elements of the Coulomb have to be multiplied by a phase depending on boxcenter.
1070 : : ! I still have to decide if it is useful to code this possibility and which variable use to
1071 : : ! define the center (boxcenter is used in the tddft part).
1072 : 46 : call wrtout(units, msg)
1073 : :
1074 : : case ('CYLINDER')
1075 [ + + ]: 8 : ii=imin_loc(ABS(Vcp%pdir-1))
1076 : 2 : write(msg,'(5a,f10.4,3a,i2,2a,3f10.2,a)')ch10,&
1077 : 2 : ' === Cylindrical cutoff === ',ch10,ch10,&
1078 : 2 : ' Cutoff radius ............... ',Vcp%rcut,' [Bohr] ',ch10,&
1079 : 2 : ' Axis parallel to direction... ',ii,ch10,&
1080 : 4 : ' Passing through point ....... ',Vcp%boxcenter,' (r.l.u) '
1081 : 2 : call wrtout(units, msg)
1082 : :
1083 : 2 : write(msg,'(2a)')' Infinite length ....... ',ch10
1084 [ + + ]: 2 : if (Vcp%hcyl/=zero) write(msg,'(a,f8.5,2a)')' Finite length of ....... ',Vcp%hcyl,' [Bohr] ',ch10
1085 : 2 : call wrtout(units, msg)
1086 : :
1087 : : CASE ('SLAB')
1088 : 0 : write(msg,'(5a,f10.4,3a,3f10.2,2a)')ch10,&
1089 : 0 : ' === Surface cutoff === ',ch10,ch10,&
1090 : 0 : ' Cutoff radius .................... ',Vcp%rcut,' [Bohr] ',ch10,&
1091 : 0 : ' Central plane passing through .... ',Vcp%boxcenter,' (r.l.u) ',ch10
1092 : 0 : call wrtout(units, msg)
1093 : : !write(msg,'(a)')' Infinite length .......'
1094 : : !if (Vcp%hcyl/=zero) write(msg,'(a,f8.5,a)')' Finite length of .......',Vcp%hcyl,' [Bohr] '
1095 : : !call wrtout(units, msg)
1096 : :
1097 : : case ('AUXILIARY_FUNCTION')
1098 : 212 : write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1099 : 212 : call wrtout(units, msg)
1100 : :
1101 : : case ('AUX_GB')
1102 : 8 : write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1103 : 8 : call wrtout(units, msg)
1104 : :
1105 : : case ('CRYSTAL')
1106 : 117 : write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1107 : 117 : call wrtout(units, msg)
1108 : :
1109 : : case ('ERF')
1110 : 1 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1111 : 1 : ' === Error function cutoff === ',ch10,ch10,&
1112 : 2 : ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1113 : 1 : call wrtout(units, msg)
1114 : :
1115 : : case ('ERFC')
1116 : 18 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1117 : 18 : ' === Complement Error function cutoff === ',ch10,ch10,&
1118 : 36 : ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1119 : 18 : call wrtout(units, msg)
1120 : :
1121 : : case default
1122 [ + + + + : 407 : ABI_BUG(sjoin('Unknown cutoff mode: ', Vcp%mode))
+ - + + +
+ + - ]
1123 : : end select
1124 : :
1125 [ + - ]: 407 : if (Vcp%nqlwl > 0) then
1126 : 407 : write(msg,'(a,i3)')" q-points for optical limit: ",Vcp%nqlwl
1127 : 407 : call wrtout(units, msg)
1128 [ + + ]: 814 : do iqlwl=1,Vcp%nqlwl
1129 : 407 : write(msg,'(1x,i5,a,2x,3f12.6)') iqlwl,')',Vcp%qlwl(:,iqlwl)
1130 : 814 : call wrtout(units, msg)
1131 : : end do
1132 : : end if
1133 : :
1134 : 407 : end subroutine vcoul_print
1135 : : !!***
1136 : :
1137 : : !----------------------------------------------------------------------
1138 : :
1139 : : !!****f* m_vcoul/vcoul_free
1140 : : !! NAME
1141 : : !! vcoul_free
1142 : : !!
1143 : : !! FUNCTION
1144 : : !! Free dynamic memory
1145 : : !!
1146 : : !! SOURCE
1147 : :
1148 : 404 : subroutine vcoul_free(Vcp)
1149 : :
1150 : : !Arguments ------------------------------------
1151 : : class(vcoul_t),intent(inout) :: Vcp
1152 : : ! *************************************************************************
1153 : :
1154 [ + - ]: 404 : ABI_SFREE(Vcp%qibz)
1155 [ + - ]: 404 : ABI_SFREE(Vcp%qlwl)
1156 [ + - ]: 404 : ABI_SFREE(Vcp%vc_sqrt)
1157 [ + - ]: 404 : ABI_SFREE(Vcp%vc_sqrt_resid)
1158 [ + - ]: 404 : ABI_SFREE(Vcp%vcqlwl_sqrt)
1159 : :
1160 : 404 : end subroutine vcoul_free
1161 : : !!***
1162 : :
1163 : : !----------------------------------------------------------------------
1164 : :
1165 : : !!****f* m_vcoul/mc_init
1166 : : !! NAME
1167 : : !! mc_init
1168 : : !!
1169 : : !! FUNCTION
1170 : : !!
1171 : : !! SOURCE
1172 : :
1173 : 3 : subroutine mc_init(mc, rprimd, ucvol, gprimd, gmet, kptrlatt)
1174 : :
1175 : : !Arguments ------------------------------------
1176 : : class(mc_t),intent(out) :: mc
1177 : : real(dp),intent(in) :: rprimd(3,3), ucvol, gprimd(3,3), gmet(3,3)
1178 : : integer,intent(in) :: kptrlatt(3,3)
1179 : :
1180 : : !Local variables-------------------------------
1181 : : integer,parameter :: ncell=3
1182 : : integer :: nseed, i1,i2,i3,imc
1183 : : real(dp) :: lmin,vlength, ucvol_sc
1184 : : real(dp) :: rprimd_sc(3,3),gprimd_sc(3,3),gmet_sc(3,3),rmet_sc(3,3), qcart2red(3,3), qtmp(3),qmin(3),qmin_cart(3)
1185 [ + - ]: 3 : integer, allocatable :: seed(:)
1186 : : ! *************************************************************************
1187 : :
1188 [ + + + + ]: 39 : mc%gmet = gmet; mc%ucvol = ucvol
1189 : :
1190 : : ! Supercell defined by the k-mesh
1191 [ + + + + : 156 : rprimd_sc(:,:) = MATMUL(rprimd, kptrlatt)
+ + + + +
+ ]
1192 : 3 : call metric(gmet_sc, gprimd_sc, -1, rmet_sc, rprimd_sc, ucvol_sc)
1193 : :
1194 [ + + + + ]: 39 : qcart2red(:,:) = two_pi * gprimd(:,:)
1195 : 3 : call matrginv(qcart2red, 3, 3)
1196 : :
1197 : : ! Find the largest sphere inside the miniBZ in order to integrate the divergence analytically
1198 : 3 : mc%q0sph = HUGE(one)
1199 [ + + ]: 21 : do i1 = -ncell+1, ncell
1200 : 18 : qtmp(1) = dble(i1) * 0.5_dp
1201 [ + + ]: 129 : do i2 = -ncell+1, ncell
1202 : 108 : qtmp(2) = dble(i2) * 0.5_dp
1203 [ + + ]: 774 : do i3 = -ncell+1, ncell
1204 : 648 : qtmp(3) = dble(i3) * 0.5_dp
1205 [ + + ]: 648 : if (i1 == 0 .AND. i2 == 0 .AND. i3 == 0) cycle
1206 : 645 : vlength = normv(qtmp, gmet_sc, 'G')
1207 [ + + ]: 753 : if (vlength < mc%q0sph) mc%q0sph = vlength
1208 : : enddo
1209 : : enddo
1210 : : enddo
1211 : :
1212 : : ! Setup the random vectors for the Monte Carlo sampling of the miniBZ at q = 0
1213 : 3 : mc%nmc_max = 2500000
1214 [ - + - + ]: 3 : ABI_MALLOC(mc%qran,(3, mc%nmc_max))
1215 : 3 : call random_seed(size=nseed)
1216 [ - + + - : 9 : ABI_MALLOC(seed, (nseed))
+ - + - -
+ - + -
+ ]
1217 [ + + ]: 27 : do i1=1,nseed
1218 : 27 : seed(i1) = NINT(SQRT(DBLE(i1) * 103731))
1219 : : end do
1220 : 3 : call random_seed(put=seed)
1221 : 3 : call random_number(mc%qran)
1222 [ - + ]: 3 : ABI_FREE(seed)
1223 : :
1224 : : ! Overide the first "random vector" with 0
1225 [ + + ]: 12 : mc%qran(:,1) = zero
1226 : :
1227 : : ! Fold qran into the Wignez-Seitz cell around q = 0
1228 [ + + ]: 7500000 : do imc=2,mc%nmc_max
1229 : : lmin = HUGE(one)
1230 [ + + ]: 52499979 : do i1 = -ncell+1, ncell
1231 : 44999982 : qtmp(1) = mc%qran(1,imc) + dble(i1)
1232 [ + + ]: 322499871 : do i2 = -ncell+1, ncell
1233 : 269999892 : qtmp(2) = mc%qran(2,imc) + dble(i2)
1234 [ + + ]: 1934999226 : do i3 = -ncell+1, ncell
1235 : 1619999352 : qtmp(3) = mc%qran(3,imc) + dble(i3)
1236 : 1619999352 : vlength = normv(qtmp, gmet_sc, 'G')
1237 [ + + ]: 1889999244 : if (vlength < lmin) then
1238 : 42272448 : lmin = vlength
1239 : : ! Get the q-vector in cartesian coordinates
1240 [ + + + + : 676359168 : qmin_cart(:) = two_pi * MATMUL( gprimd_sc(:,:) , qtmp )
+ + ]
1241 : : ! Transform it back to the reciprocal space
1242 [ + + + + ]: 549541824 : qmin(:) = MATMUL(qcart2red , qmin_cart)
1243 : : end if
1244 : : enddo
1245 : : enddo
1246 : : enddo
1247 : :
1248 [ + + ]: 29999991 : mc%qran(:,imc) = qmin(:)
1249 : : enddo
1250 : :
1251 [ - + ]: 3 : end subroutine mc_init
1252 : : !!***
1253 : :
1254 : : !----------------------------------------------------------------------
1255 : :
1256 : : !!****f* m_vcoul/mc_integrate
1257 : : !! NAME
1258 : : !! mc_integrate
1259 : : !!
1260 : : !! FUNCTION
1261 : : !!
1262 : : !! SOURCE
1263 : :
1264 : 12 : subroutine mc_integrate(mc, mode, qibz, ng, gvec, rcut2, nkbz, vcoul, comm)
1265 : :
1266 : : !Arguments ------------------------------------
1267 : : class(mc_t),intent(in) :: mc
1268 : : real(dp),intent(in) :: rcut2
1269 : : integer,intent(in) :: nkbz, ng, comm
1270 : : character(len=*),intent(in) :: mode
1271 : : real(dp),intent(in) :: qibz(3)
1272 : : integer,intent(in) :: gvec(3, ng)
1273 : : real(dp),intent(out) :: vcoul(ng)
1274 : :
1275 : : !Local variables-------------------------------
1276 : : integer,parameter :: master = 0
1277 : : integer :: ig, ig0, imc, nmc, my_rank, nprocs, ierr
1278 : : logical :: q_is_gamma
1279 : : real(dp) :: qpg2, qpg(3)
1280 : : ! *************************************************************************
1281 : :
1282 : 12 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1283 [ + + + + ]: 21 : q_is_gamma = all(abs(qibz) < tol16)
1284 : :
1285 : : ! Find index of G=0 in gvec.
1286 : 12 : ig0 = -1
1287 [ + - ]: 12 : do ig=1,ng
1288 [ + + + - : 48 : if (all(gvec(:, ig) == 0)) then
- + ]
1289 : : ig0 = ig; exit
1290 : : end if
1291 : : end do
1292 [ - + ]: 12 : ABI_CHECK(ig0 /= -1, "Cannot find G=0 in gvec!")
1293 : :
1294 [ + + ]: 3408 : vcoul = zero
1295 : :
1296 : 12 : select case (trim(mode))
1297 : :
1298 : : case ('MINIBZ')
1299 [ + + ]: 1136 : do ig=1,ng
1300 [ + - ]: 1132 : if (mod(ig, nprocs) /= my_rank) cycle ! MPI parallelism.
1301 [ + + ]: 1132 : if (q_is_gamma .and. ig == ig0) cycle
1302 [ + + ]: 4524 : qpg(:) = qibz(:) + gvec(:,ig)
1303 : 1131 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1304 : 1131 : nmc = adapt_nmc(mc%nmc_max, qpg2)
1305 [ + + ]: 43850951 : do imc=1,nmc
1306 [ + + ]: 175399264 : qpg(:) = qibz(:) + gvec(:,ig) + mc%qran(:,imc)
1307 : 43849816 : qpg2 = normv(qpg, mc%gmet, 'G') ** 2
1308 : 43850948 : vcoul(ig) = vcoul(ig) + four_pi / qpg2 / REAL(nmc, dp)
1309 : : end do
1310 : : end do ! ig
1311 : :
1312 [ + + ]: 4 : if (q_is_gamma .and. my_rank == master) then
1313 : : ! Compute ig0 component
1314 : 1 : vcoul(ig0) = four_pi**2 * nkbz * mc%ucvol / ( 8.0_dp * pi**3 ) * mc%q0sph
1315 [ + + ]: 2500001 : do imc=1,mc%nmc_max
1316 [ + + ]: 10000000 : qpg(:) = qibz(:) + gvec(:,ig0) + mc%qran(:,imc)
1317 : 2500000 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1318 [ + + ]: 2500001 : if (qpg2 > mc%q0sph ** 2) vcoul(ig0) = vcoul(ig0) + four_pi / qpg2 / REAL(mc%nmc_max, dp)
1319 : : end do
1320 : : end if
1321 : :
1322 : : case ('MINIBZ-ERFC')
1323 [ + + ]: 1136 : do ig=1,ng
1324 [ + - ]: 1132 : if (mod(ig, nprocs) /= my_rank) cycle ! MPI parallelism.
1325 [ + + ]: 1132 : if (q_is_gamma .and. ig == ig0) cycle
1326 [ + + ]: 4524 : qpg(:) = qibz(:) + gvec(:,ig)
1327 : 1131 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1328 : 1131 : nmc = adapt_nmc(mc%nmc_max, qpg2)
1329 [ + + ]: 43850951 : do imc=1,nmc
1330 [ + + ]: 175399264 : qpg(:) = qibz(:) + gvec(:,ig) + mc%qran(:,imc)
1331 : 43849816 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1332 : 43850948 : vcoul(ig) = vcoul(ig) + four_pi / qpg2 / REAL(nmc,dp) * ( one - EXP( -0.25d0 * rcut2 * qpg2 ) )
1333 : : end do
1334 : : end do ! ig
1335 : :
1336 [ + + ]: 4 : if (q_is_gamma .and. my_rank == master) then
1337 : : ! Override ig0 component
1338 : : vcoul(ig0) = four_pi**2 * nkbz * mc%ucvol / ( 8.0_dp * pi**3 ) &
1339 : 1 : * ( mc%q0sph - SQRT(pi/rcut2) * abi_derf(0.5_dp*SQRT(rcut2)*mc%q0sph) )
1340 [ + + ]: 2500001 : do imc=1,mc%nmc_max
1341 [ + + ]: 10000000 : qpg(:) = qibz(:) + gvec(:,ig0) + mc%qran(:,imc)
1342 : 2500000 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1343 [ + + ]: 2500001 : if (qpg2 > mc%q0sph**2) then
1344 : 800113 : vcoul(ig0) = vcoul(ig0) + four_pi / qpg2 / REAL(mc%nmc_max,dp) * (one - EXP( -0.25d0 * rcut2 * qpg2))
1345 : : end if
1346 : : end do
1347 : : end if
1348 : :
1349 : : case ('MINIBZ-ERF')
1350 [ + + ]: 1136 : do ig=1,ng
1351 [ + - ]: 1132 : if (mod(ig, nprocs) /= my_rank) cycle ! MPI parallelism.
1352 [ + + ]: 1132 : if (q_is_gamma .and. ig == ig0) cycle
1353 [ + + ]: 4524 : qpg(:) = qibz(:) + gvec(:,ig)
1354 : 1131 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1355 : 1131 : nmc = adapt_nmc(mc%nmc_max, qpg2)
1356 [ + + ]: 43850951 : do imc=1,nmc
1357 [ + + ]: 175399264 : qpg(:) = qibz(:) + gvec(:,ig) + mc%qran(:,imc)
1358 : 43849816 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1359 : 43850948 : vcoul(ig) = vcoul(ig) + four_pi / qpg2 / REAL(nmc,dp) * EXP( -0.25d0 * rcut2 * qpg2 )
1360 : : end do
1361 : : end do ! ig
1362 : :
1363 [ + + ]: 4 : if (q_is_gamma .and. my_rank == master) then
1364 : : ! Override ig=ig0 component
1365 : 1 : vcoul(ig0) = four_pi**2 * nkbz * mc%ucvol / ( 8.0_dp * pi**3 ) * SQRT(pi/rcut2) * abi_derf(0.5_dp*SQRT(rcut2)*mc%q0sph)
1366 : :
1367 [ + + ]: 2500001 : do imc=1,mc%nmc_max
1368 [ + + ]: 10000000 : qpg(:) = qibz(:) + gvec(:,ig0) + mc%qran(:,imc)
1369 : 2500000 : qpg2 = normv(qpg, mc%gmet, 'G')**2
1370 [ + + ]: 2500001 : if (qpg2 > mc%q0sph**2) then
1371 : 800113 : vcoul(ig0) = vcoul(ig0) + four_pi / qpg2 / REAL(mc%nmc_max,dp) * EXP( -0.25d0 * rcut2 * qpg2 )
1372 : : end if
1373 : : end do
1374 : : end if
1375 : :
1376 : : case default
1377 [ + + + - ]: 12 : ABI_ERROR(sjoin("Invalid mode:", mode))
1378 : : end select
1379 : :
1380 : : ! Collect result on each MPI proc.
1381 : 12 : call xmpi_sum(vcoul, comm, ierr)
1382 : :
1383 : 12 : end subroutine mc_integrate
1384 : : !!***
1385 : :
1386 : : !----------------------------------------------------------------------
1387 : :
1388 : : !!****f* m_vcoul/mc_free
1389 : : !! NAME
1390 : : !! mc_free
1391 : : !!
1392 : : !! FUNCTION
1393 : : !! Free dynamic memory
1394 : : !!
1395 : : !! SOURCE
1396 : :
1397 : 3 : subroutine mc_free(mc)
1398 : :
1399 : : !Arguments ------------------------------------
1400 : : class(mc_t),intent(inout) :: mc
1401 : : ! *************************************************************************
1402 : :
1403 [ + - ]: 3 : ABI_SFREE(mc%qran)
1404 : :
1405 : 3 : end subroutine mc_free
1406 : : !!***
1407 : :
1408 : : !----------------------------------------------------------------------
1409 : :
1410 : : !!****f* m_vcoul/beigi_cylinder_limit
1411 : : !! NAME
1412 : : !! beigi_cylinder_limit
1413 : : !!
1414 : : !! FUNCTION
1415 : : !!
1416 : : !! SOURCE
1417 : :
1418 : 1 : subroutine beigi_cylinder_limit(opt_cylinder, cryst, nqibz, nkbz, rcut, hcyl, boxcenter, pdir, i_sz)
1419 : :
1420 : : !Arguments ------------------------------------
1421 : : integer,intent(in) :: opt_cylinder, nqibz, nkbz, pdir(3)
1422 : : type(crystal_t),intent(in) :: cryst
1423 : : real(dp),intent(in) :: rcut, hcyl, boxcenter(3)
1424 : : real(dp),intent(out) :: i_sz
1425 : :
1426 : : !Local variables-------------------------------
1427 : : integer :: ii, iq, npar, npt, gamma_pt(3,1)
1428 : : real(dp) :: step, bz_plane, dx, integ, q0_vol, q0_volsph, b1(3),b2(3),b3(3)
1429 : 1 : real(dp),allocatable :: cov(:,:),par(:),qfit(:,:),sigma(:),var(:), vcfit(:,:),xx(:),yy(:)
1430 : : ! *************************************************************************
1431 : :
1432 [ + + + + ]: 7 : b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
1433 : :
1434 : 1 : npt =100
1435 : 1 : npar=8; gamma_pt = RESHAPE(([0, 0, 0]), [3, 1])
1436 [ - + ]: 1 : ABI_MALLOC(qfit, (3, npt))
1437 [ - + ]: 1 : ABI_MALLOC(vcfit, (1, npt))
1438 [ - + ]: 1 : if (nqibz == 1) then
1439 : 0 : ABI_ERROR("nqibz == 1 not supported when Beigi's method is used")
1440 : : endif
1441 [ + + + + ]: 401 : qfit(:,:)=zero
1442 [ + + ]: 101 : step=half/(npt * (nqibz-1)) ; qfit(3,:)=arth(tol6,step,npt)
1443 : : !step=(half/(nqibz-1)/tol6)**(one/npt) ; qfit(3,:)=geop(tol6,step,npt)
1444 : :
1445 [ + + ]: 101 : do iq=1,npt
1446 : : call cutoff_cylinder(qfit(:,iq),1,gamma_pt,rcut,hcyl,pdir,boxcenter,&
1447 : 101 : Cryst%rprimd,vcfit(:,iq),opt_cylinder, xmpi_comm_self)
1448 : : end do
1449 : :
1450 [ - + ]: 1 : ABI_MALLOC(xx, (npt))
1451 [ - + ]: 1 : ABI_MALLOC(yy, (npt))
1452 [ - + ]: 1 : ABI_MALLOC(sigma, (npt))
1453 [ - + ]: 1 : ABI_MALLOC(par, (npar))
1454 [ - + ]: 1 : ABI_MALLOC(var, (npar))
1455 [ - + ]: 1 : ABI_MALLOC(cov, (npar, npar))
1456 : :
1457 [ + + ]: 101 : do ii=1,npt
1458 : 101 : xx(ii) = normv(qfit(:,ii), cryst%gmet, 'G')
1459 : : end do
1460 [ - + ]: 1 : ABI_FREE(qfit)
1461 [ + + + + ]: 201 : sigma=one ; yy(:)=vcfit(1,:)
1462 [ - + ]: 1 : ABI_FREE(vcfit)
1463 : : !call llsfit_svd(xx,yy,sigma,npar,K0fit,chisq,par,var,cov,info)
1464 : : !do ii=1,npt
1465 : : ! write(99,*)xx(ii),yy(ii),DOT_PRODUCT(par,K0fit(xx(ii),npar))
1466 : : !end do
1467 : 1 : bz_plane=l2norm(b1.x.b2)
1468 : : !integ=K0fit_int(xx(npt),par,npar)
1469 : : !write(std_out,*)' SVD fit : chi-square',chisq
1470 : : !write(std_out,*)' fit-parameters : ',par
1471 : : !write(std_out,*)' variance ',var
1472 : : !write(std_out,*)' bz_plane ',bz_plane
1473 : : !write(std_out,*)' SCD integ ',integ
1474 : : ! Here Im assuming homogeneous mesh
1475 : 1 : dx=(xx(2)-xx(1))
1476 : 1 : integ=yy(2)*dx*3.0/2.0
1477 [ + + ]: 97 : integ=integ+SUM(yy(3:npt-2))*dx
1478 : 1 : integ=integ+yy(npt-1)*dx*3.0/2.0
1479 : : !write(std_out,*)' simple integral',integ
1480 : 1 : q0_volsph = (two_pi)**3 / (nkbz * cryst%ucvol)
1481 : 1 : q0_vol=bz_plane*two*xx(npt)
1482 : : !write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
1483 : 1 : i_sz = bz_plane * two * integ / q0_vol
1484 : : !write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
1485 : : !write(std_out,*)' Cylindrical cutoff value ',i_sz
1486 : : !i_sz=four_pi*7.44*q0_vol**(-two_thirds)
1487 : :
1488 : 1 : ABI_FREE(xx)
1489 : 1 : ABI_FREE(yy)
1490 [ - + ]: 1 : ABI_FREE(sigma)
1491 [ - + ]: 1 : ABI_FREE(par)
1492 [ - + ]: 1 : ABI_FREE(var)
1493 [ - + ]: 1 : ABI_FREE(cov)
1494 : :
1495 : 1 : end subroutine beigi_cylinder_limit
1496 : : !!***
1497 : :
1498 : : !----------------------------------------------------------------------
1499 : :
1500 : : !!****f* m_vcoul/beigi_surface_limit
1501 : : !! NAME
1502 : : !! beigi_surface_limit
1503 : : !!
1504 : : !! FUNCTION
1505 : : !!
1506 : : !! SOURCE
1507 : :
1508 : 0 : subroutine beigi_surface_limit(opt_slab, cryst, nqibz, nkbz, rcut, alpha, boxcenter, pdir, i_sz)
1509 : :
1510 : : !Arguments ------------------------------------
1511 : : integer,intent(in) :: opt_slab, nqibz, nkbz, pdir(3)
1512 : : type(crystal_t),intent(in) :: cryst
1513 : : real(dp),intent(in) :: rcut, alpha(3), boxcenter(3)
1514 : : real(dp),intent(out) :: i_sz
1515 : :
1516 : : !Local variables-------------------------------
1517 : : integer :: ii, npt, gamma_pt(3,1)
1518 : : real(dp) :: step, bz_plane, dx, integ, q0_vol, q0_volsph, b1(3),b2(3),b3(3)
1519 : 0 : real(dp),allocatable :: qfit(:,:),sigma(:),vcfit(:,:),xx(:),yy(:), qcart(:,:)
1520 : : ! *************************************************************************
1521 : :
1522 [ # # ]: 0 : b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
1523 : :
1524 : 0 : gamma_pt=RESHAPE([0, 0, 0], [3, 1]) ! Gamma point
1525 : 0 : npt=100 ! Number of points in 1D
1526 [ # # ]: 0 : ABI_MALLOC(qfit, (3, npt))
1527 [ # # ]: 0 : ABI_MALLOC(qcart, (3, npt))
1528 [ # # ]: 0 : ABI_MALLOC(vcfit, (1, npt))
1529 [ # # ]: 0 : if (nqibz == 1) then
1530 : 0 : ABI_ERROR("nqibz == 1 not supported when Beigi's method is used")
1531 : : endif
1532 [ # # # # ]: 0 : qfit(:,:)=zero
1533 [ # # # # ]: 0 : qcart(:,:)=zero
1534 : : ! Size of the third vector
1535 : 0 : bz_plane=l2norm(b3)
1536 : 0 : q0_volsph=(two_pi)**3 / (nkbz * cryst%ucvol)
1537 : : ! radius that gives the same volume as q0_volsph
1538 : : ! Let's assume that c is perpendicular to the plane
1539 : : ! We also assume isotropic BZ around gamma
1540 : 0 : step=sqrt((q0_volsph/bz_plane)/pi)/npt
1541 : :
1542 : : !step=half/(npt*(nqibz-1))
1543 : : ! Let's take qpoints along 1 line, the vcut does depend only on the norm
1544 [ # # ]: 0 : qcart(1,:) = arth(tol6,step,npt)
1545 : :
1546 [ # # ]: 0 : do ii=1,npt
1547 [ # # ]: 0 : qfit(:,ii) = MATMUL(TRANSPOSE(Cryst%rprimd),qcart(:,ii)) / (2*pi)
1548 : : call cutoff_slab(qfit(:,ii), 1, gamma_pt, cryst%gprimd, rcut, &
1549 : 0 : boxcenter, pdir, alpha, vcfit(:,ii), opt_slab)
1550 : : end do
1551 : :
1552 [ # # ]: 0 : ABI_MALLOC(xx, (npt))
1553 [ # # ]: 0 : ABI_MALLOC(yy, (npt))
1554 [ # # ]: 0 : ABI_MALLOC(sigma, (npt))
1555 [ # # ]: 0 : do ii=1,npt
1556 : : !xx(ii)=qfit(1,:)
1557 : 0 : xx(ii) = normv(qfit(:,ii), cryst%gmet, 'G')
1558 : : end do
1559 [ # # ]: 0 : ABI_FREE(qfit)
1560 [ # # ]: 0 : sigma=one
1561 [ # # ]: 0 : yy(:)=vcfit(1,:)
1562 : : !yy(:)=one
1563 [ # # ]: 0 : ABI_FREE(vcfit)
1564 : :
1565 : 0 : dx=(xx(2)-xx(1))
1566 : : ! integ = \int dr r f(r)
1567 : 0 : integ=xx(2)*yy(2)*dx*3.0/2.0
1568 [ # # ]: 0 : integ=integ+DOT_PRODUCT(xx(3:npt-2),yy(3:npt-2))*dx
1569 : 0 : integ=integ+xx(npt-1)*yy(npt-1)*dx*3.0/2.0
1570 : : !write(std_out,*)' simple integral',integ
1571 : 0 : q0_vol=bz_plane*pi*xx(npt)**2
1572 : : !write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
1573 : 0 : i_sz=bz_plane*2*pi*integ/q0_vol
1574 : : !write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
1575 : : !write(std_out,*)' Cylindrical cutoff value ',i_sz
1576 : : !i_sz=four_pi*7.44*q0_vol**(-two_thirds)
1577 : 0 : ABI_FREE(xx)
1578 : 0 : ABI_FREE(yy)
1579 : :
1580 [ # # # # ]: 0 : end subroutine beigi_surface_limit
1581 : : !!***
1582 : :
1583 : : !!****f* m_vcoul/carrier_isz
1584 : : !! NAME
1585 : : !! carrier_isz
1586 : : !!
1587 : : !! FUNCTION
1588 : : !!
1589 : : !! SOURCE
1590 : :
1591 : 210 : real(dp) function carrier_isz(cryst, nqbz, qbz, rcut, comm) result(i_sz)
1592 : :
1593 : : !Arguments ------------------------------------
1594 : : type(crystal_t),intent(in) :: cryst
1595 : : integer,intent(in) :: nqbz, comm
1596 : : real(dp), intent(in) :: qbz(3, nqbz), rcut
1597 : :
1598 : : !Local variables-------------------------------
1599 : : integer :: iq_bz
1600 : : real(dp) :: qbz_norm, bz_geometry_factor, qbz_cart(3), b1(3), b2(3), b3(3)
1601 : : !************************************************************************
1602 : :
1603 [ + + + + : 2100 : b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
+ + ]
1604 : :
1605 : : bz_geometry_factor = zero
1606 [ + + ]: 9777 : do iq_bz=1,nqbz
1607 [ + + ]: 38268 : qbz_cart(:) = qbz(1,iq_bz)*b1(:) + qbz(2,iq_bz)*b2(:) + qbz(3,iq_bz)*b3(:)
1608 [ + + + + : 38268 : qbz_norm = NORM2(qbz_cart(:))
+ + ]
1609 [ + + ]: 9777 : if (qbz_norm > TOLQ0) bz_geometry_factor = bz_geometry_factor - faux(qbz(:,iq_bz), rcut, b1, b2, b3)
1610 : : end do
1611 : :
1612 : 210 : bz_geometry_factor = bz_geometry_factor + integratefaux(rcut, cryst%gprimd, cryst%ucvol, comm) * nqbz
1613 : 210 : i_sz = four_pi * bz_geometry_factor ! Final result stored here
1614 : :
1615 : 210 : end function carrier_isz
1616 : : !!***
1617 : :
1618 : : !!****f* m_vcoul/gygi_baldereschi_isz
1619 : : !! NAME
1620 : : !! gygi_baldereschi_isz
1621 : : !!
1622 : : !! FUNCTION
1623 : : !!
1624 : : !! SOURCE
1625 : :
1626 : 8 : real(dp) function gygi_baldereschi_isz(cryst, nqbz, qbz, vc_ecut, ng, gvec) result(i_sz)
1627 : :
1628 : : !Arguments ------------------------------------
1629 : : type(crystal_t),intent(in) :: cryst
1630 : : integer,intent(in) :: nqbz, ng
1631 : : real(dp), intent(in) :: qbz(3, nqbz), vc_ecut
1632 : : integer,intent(in) :: gvec(3,ng)
1633 : :
1634 : : !Local variables-------------------------------
1635 : : integer :: iq_bz, ig
1636 : : real(dp) :: bz_geometry_factor, intfauxgb, alfa, qpg2, qpg(3)
1637 : : !************************************************************************
1638 : :
1639 : : ! the choice of alfa (the width of the gaussian) is somehow empirical
1640 : 8 : alfa = 150.0 / vc_ecut
1641 : :
1642 : 8 : bz_geometry_factor=zero
1643 [ + + ]: 224 : do iq_bz=1,nqbz
1644 [ + + ]: 59192 : do ig = 1,ng
1645 [ + + ]: 235872 : qpg(:) = qbz(:,iq_bz) + gvec(:,ig)
1646 : 58968 : qpg2 = normv(qpg, cryst%gmet, 'G')**2
1647 [ + + ]: 59184 : if (qpg2 > TOLQ0) bz_geometry_factor = bz_geometry_factor - EXP(-alfa*qpg2)/qpg2
1648 : : end do
1649 : : end do
1650 : :
1651 : 8 : intfauxgb = cryst%ucvol/four_pi/SQRT(0.5*two_pi*alfa)
1652 : 8 : bz_geometry_factor = bz_geometry_factor + intfauxgb * nqbz
1653 : :
1654 : 8 : i_sz = four_pi*bz_geometry_factor
1655 : :
1656 : 8 : end function gygi_baldereschi_isz
1657 : : !!***
1658 : :
1659 : : !----------------------------------------------------------------------
1660 : :
1661 : : !!****f* m_vcoul/vcgen_init
1662 : : !! NAME
1663 : : !! vcgen_init
1664 : : !!
1665 : : !! FUNCTION
1666 : : !!
1667 : : !! SOURCE
1668 : :
1669 [ # # ]: 0 : subroutine vcgen_init(vcgen, cryst, kptrlatt, nkbz, nqibz, nqbz, qbz, rcut, gw_icutcoul, vcutgeo, vc_ecut, comm)
1670 : :
1671 : : !Arguments ------------------------------------
1672 : : class(vcgen_t),intent(out) :: vcgen
1673 : : type(crystal_t),intent(in) :: cryst
1674 : : integer,intent(in) :: kptrlatt(3,3), nkbz, nqibz, nqbz, gw_icutcoul
1675 : : real(dp),intent(in) :: qbz(3,nqbz), rcut, vc_ecut, vcutgeo(3)
1676 : : integer,intent(in) :: comm
1677 : :
1678 : : !Local variables-------------------------------
1679 : : integer,parameter :: istwfk1 = 1
1680 : : integer :: gvec0(3) ! npw_,
1681 : : real(dp) :: q0_vol, bz_geometry_factor, rcut2
1682 : : character(len=500) :: msg
1683 : : real(dp) :: vcoul0(1), q_gamma(3)
1684 : : !integer,allocatable :: gvec_(:,:)
1685 : : ! *************************************************************************
1686 : :
1687 : : ABI_UNUSED([vc_ecut])
1688 : :
1689 : : ! Save dimension and other useful quantities in Vcp
1690 : 0 : vcgen%rcut = rcut ! Cutoff radius for cylinder.
1691 : 0 : vcgen%hcyl = zero ! Length of finite cylinder (Rozzi's method, default is Beigi).
1692 [ # # ]: 0 : vcgen%boxcenter = zero ! Boxcenter at the moment is supposed to be at the origin.
1693 [ # # ]: 0 : vcgen%vcutgeo = vcutgeo(:) ! Info on the orientation and extension of the cutoff region.
1694 : 0 : vcgen%nkbz = nkbz
1695 : :
1696 : 0 : call gw_icutcoul_to_mode(gw_icutcoul, vcgen%mode)
1697 : 0 : q_gamma = zero
1698 : 0 : gvec0 = 0
1699 : :
1700 : 0 : select case (trim(vcgen%mode))
1701 : : case ('MINIBZ', 'MINIBZ-ERFC', 'MINIBZ-ERF')
1702 : 0 : call vcgen%mc%init(cryst%rprimd, cryst%ucvol, cryst%gprimd, cryst%gmet, kptrlatt)
1703 : 0 : rcut2 = vcgen%rcut**2
1704 : 0 : call vcgen%mc%integrate(vcgen%mode, q_gamma, 1, gvec0, rcut2, nkbz, vcoul0, xmpi_comm_self)
1705 : : ! Treat the limit q --> 0.
1706 : 0 : vcgen%i_sz = vcoul0(1)
1707 : :
1708 : : case ('SPHERE')
1709 : : ! A non-positive value of rcut activates the recipe of Spencer & Alavi, PRB 77, 193110 (2008) [[cite:Spencer2008]].
1710 [ # # ]: 0 : if (vcgen%rcut < tol12) then
1711 : 0 : vcgen%rcut = (cryst%ucvol * nkbz * 3.d0 / four_pi) ** third
1712 : 0 : write(msg,'(2a,2x,f8.4,a)')ch10,' Using calculated rcut: ',vcgen%rcut,' to have same volume as the BvK crystal'
1713 : 0 : call wrtout(std_out, msg)
1714 : : end if
1715 [ # # ]: 0 : vcgen%vcutgeo = zero
1716 : :
1717 : : ! Treat the limit q --> 0
1718 : : ! The small cube is approximated by a sphere, while vc(q=0) = 2piR**2.
1719 : : ! if a single q-point is used, the expression for the volume is exact.
1720 : 0 : vcgen%i_sz = two_pi * vcgen%rcut**2
1721 : :
1722 : : case ('CYLINDER')
1723 : 0 : call cylinder_setup(cryst, vcgen%vcutgeo, vcgen%hcyl, vcgen%pdir, vcgen%opt_cylinder)
1724 : :
1725 : : ! If Beigi, treat the limit q --> 0.
1726 [ # # ]: 0 : if (vcgen%opt_cylinder == 1) then
1727 : : call beigi_cylinder_limit(vcgen%opt_cylinder, cryst, nqibz, nkbz, &
1728 : 0 : vcgen%rcut, vcgen%hcyl, vcgen%boxcenter, vcgen%pdir, vcgen%i_sz)
1729 : : else
1730 : : ! In Rozzi's method the lim q+G --> 0 is finite.
1731 : : call cutoff_cylinder(q_gamma, 1, gvec0, vcgen%rcut, vcgen%hcyl, vcgen%pdir,&
1732 : 0 : vcgen%boxcenter, cryst%rprimd, vcoul0, vcgen%opt_cylinder, xmpi_comm_self)
1733 : 0 : vcgen%i_sz = vcoul0(1)
1734 : : end if
1735 : :
1736 : : case ('SLAB')
1737 : 0 : call surface_setup(cryst, vcgen%vcutgeo, vcgen%alpha, vcgen%rcut, vcgen%pdir, vcgen%opt_slab)
1738 : :
1739 : : ! If Beigi, treat the limit q --> 0.
1740 [ # # ]: 0 : if (vcgen%opt_slab == 1) then
1741 : : ! Integrate numerically in the plane close to 0
1742 : : call beigi_surface_limit(vcgen%opt_slab, cryst, nqibz, nkbz, vcgen%rcut, vcgen%alpha, &
1743 : 0 : vcgen%boxcenter, vcgen%pdir, vcgen%i_sz)
1744 : : else
1745 : : ! In Rozzi's method the lim q+G --> 0 is finite.
1746 : : call cutoff_slab(q_gamma, 1, gvec0, cryst%gprimd, vcgen%rcut, &
1747 : 0 : vcgen%boxcenter, vcgen%pdir, vcgen%alpha, vcoul0, vcgen%opt_slab)
1748 : 0 : vcgen%i_sz = vcoul0(1)
1749 : : end if
1750 : :
1751 : : case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB")
1752 : :
1753 [ # # ]: 0 : if (vcgen%mode == "CRYSTAL") then
1754 : : ! Analytic integration of 4pi/q^2 over the volume element:
1755 : : ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
1756 : : ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
1757 : : ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (see gwa.pdf, appendix A.4)
1758 : 0 : q0_vol = (two_pi) **3 / (nkbz*cryst%ucvol); bz_geometry_factor=zero
1759 : 0 : vcgen%i_sz = four_pi*7.44*q0_vol**(-two_thirds)
1760 : :
1761 [ # # ]: 0 : else if (vcgen%mode == "AUXILIARY_FUNCTION") then
1762 : : ! Numerical integration of the exact-exchange divergence through the
1763 : : ! auxiliary function of Carrier et al. PRB 75, 205126 (2007) [[cite:Carrier2007]].
1764 : 0 : vcgen%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, comm)
1765 : :
1766 [ # # ]: 0 : else if (vcgen%mode == "AUX_GB") then
1767 : : ! We use the auxiliary function of a Gygi-Baldereschi variant [[cite:Gigy1986]]
1768 : : ! TODO:
1769 : 0 : ABI_ERROR("AUX_GB not implemented in vcgen_init")
1770 : : !call get_kg(kk_bz, istwfk1, vc_ecut, cryst%gmet, npw_, gvec_)
1771 : : !vcgen%i_sz = gygi_baldereschi_isz(cryst, nqbz, qbz, vc_ecut, ng, gvec_)
1772 : : !ABI_FREE(gvec_)
1773 : :
1774 : : else
1775 : 0 : ABI_ERROR(sjoin("Need treatment of 1/q^2 singularity! for mode", vcgen%mode))
1776 : : end if
1777 : :
1778 : : case ('ERF')
1779 : 0 : vcgen%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, xmpi_comm_self)
1780 : :
1781 : : case ('ERFC')
1782 : : ! === Treat 1/q^2 singularity ===
1783 : : ! * There is NO singularity in this case.
1784 : 0 : vcgen%i_sz = pi * vcgen%rcut**2 ! Final result stored here
1785 : :
1786 : : case default
1787 [ # # # # : 0 : ABI_BUG(sjoin('Unsupported cutoff mode:', vcgen%mode))
# # # # ]
1788 : : end select
1789 : :
1790 : : !call wrtout(std_out, sjoin(" vcgen%i_sz", ftoa(vcgen%i_sz)))
1791 : :
1792 : 0 : end subroutine vcgen_init
1793 : : !!***
1794 : :
1795 : : !----------------------------------------------------------------------
1796 : :
1797 : : !!****f* m_vcoul/vcgen_get_vc_sqrt
1798 : : !! NAME
1799 : : !! vcgen_get_vc_sqrt
1800 : : !!
1801 : : !! FUNCTION
1802 : : !! Compute sqrt(vc(q, g))
1803 : : !!
1804 : : !! SOURCE
1805 : :
1806 : 0 : subroutine vcgen_get_vc_sqrt(vcgen, qpt, npw, gvec, q0, cryst, vc_sqrt, comm, &
1807 : 0 : vc) ! optional
1808 : :
1809 : : !Arguments ------------------------------------
1810 : : class(vcgen_t),intent(in) :: vcgen
1811 : : real(dp),intent(in) :: qpt(3), q0(3)
1812 : : integer,intent(in) :: npw, gvec(3,npw), comm
1813 : : type(crystal_t),intent(in) :: cryst
1814 : : complex(gwp),intent(out) :: vc_sqrt(npw)
1815 : : real(dp),optional,intent(out) :: vc(npw)
1816 : :
1817 : : !Local variables-------------------------------
1818 : : integer :: ig, ig0
1819 : : real(dp) :: rcut2
1820 : : logical :: q_is_gamma
1821 : 0 : real(dp),allocatable :: vcoul(:)
1822 : : ! *************************************************************************
1823 : :
1824 : 0 : q_is_gamma = normv(qpt, cryst%gmet, "G") < GW_TOLQ0
1825 : :
1826 : : ! Find the index of G=0 in gvec.
1827 : 0 : ig0 = -1
1828 [ # # ]: 0 : do ig=1,npw
1829 [ # # # # : 0 : if (all(gvec(:,ig) == 0)) then
# # ]
1830 : : ig0 = ig; exit
1831 : : end if
1832 : : end do
1833 [ # # ]: 0 : ABI_CHECK(ig0 /= -1, "Cannot find G=0 in gvec!")
1834 : :
1835 [ # # # # : 0 : ABI_MALLOC(vcoul, (npw))
# # # # #
# # # ]
1836 : :
1837 : 0 : select case (trim(vcgen%mode))
1838 : : case ('MINIBZ', 'MINIBZ-ERFC', 'MINIBZ-ERF')
1839 : 0 : rcut2 = vcgen%rcut**2
1840 : 0 : call vcgen%mc%integrate(vcgen%mode, qpt, npw, gvec, rcut2, vcgen%nkbz, vcoul, comm)
1841 : :
1842 : : case ('SPHERE')
1843 : 0 : call cutoff_sphere(qpt, npw, gvec, cryst%gmet, vcgen%rcut, vcoul)
1844 : :
1845 : : case ('CYLINDER')
1846 : : call cutoff_cylinder(qpt, npw, gvec, vcgen%rcut, vcgen%hcyl, vcgen%pdir, &
1847 : 0 : vcgen%boxcenter, cryst%rprimd, vcoul, vcgen%opt_cylinder, comm)
1848 : :
1849 : : case ('SLAB')
1850 : : call cutoff_slab(qpt, npw, gvec, cryst%gprimd, vcgen%rcut, &
1851 : 0 : vcgen%boxcenter, vcgen%pdir, vcgen%alpha, vcoul, vcgen%opt_slab)
1852 : :
1853 : : case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB", "ERF", "ERFC")
1854 : : ! Compute |q+G| with special treatment of (q=0, g=0).
1855 [ # # ]: 0 : do ig=1,npw
1856 [ # # ]: 0 : if (q_is_gamma .and. ig == ig0) then
1857 [ # # ]: 0 : vcoul(ig) = normv(q0 + gvec(:,ig), cryst%gmet, "G")
1858 : : !print *, "q_is_gamma with ", q0, "and vcoul:", vcoul(ig); stop
1859 : : else
1860 [ # # ]: 0 : vcoul(ig) = normv(qpt + gvec(:,ig), cryst%gmet, "G")
1861 : : end if
1862 : : end do
1863 : :
1864 [ # # ]: 0 : if (vcgen%mode == "ERF") then
1865 [ # # ]: 0 : vcoul(:) = four_pi/(vcoul(:)**2) * EXP( -0.25d0 * (vcgen%rcut*vcoul(:))**2 )
1866 [ # # ]: 0 : else if (vcgen%mode == "ERFC") then
1867 [ # # ]: 0 : vcoul(:) = four_pi/(vcoul(:)**2) * ( one - EXP( -0.25d0 * (vcgen%rcut*vcoul(:))**2 ) )
1868 : : else
1869 [ # # ]: 0 : vcoul = four_pi/vcoul**2
1870 : : end if
1871 : :
1872 : : case default
1873 [ # # # # : 0 : ABI_BUG(sjoin('Unsupported cutoff mode:', vcgen%mode))
# # ]
1874 : : end select
1875 : :
1876 : : ! Store final results in complex array as Rozzi's cutoff can give real negative values
1877 [ # # ]: 0 : vc_sqrt = sqrt(cmplx(vcoul, zero))
1878 : :
1879 [ # # # # ]: 0 : if (present(vc)) vc = vcoul
1880 [ # # ]: 0 : ABI_FREE(vcoul)
1881 : :
1882 [ # # ]: 0 : end subroutine vcgen_get_vc_sqrt
1883 : : !!***
1884 : :
1885 : : !----------------------------------------------------------------------
1886 : :
1887 : : !!****f* m_vcoul/vcgen_free
1888 : : !! NAME
1889 : : !! vcgen_free
1890 : : !!
1891 : : !! FUNCTION
1892 : : !! Free dynamic memory.
1893 : : !!
1894 : : !! SOURCE
1895 : :
1896 : 0 : subroutine vcgen_free(vcgen)
1897 : :
1898 : : !Arguments ------------------------------------
1899 : : class(vcgen_t),intent(inout) :: vcgen
1900 : : ! *************************************************************************
1901 : :
1902 : 0 : call vcgen%mc%free()
1903 : :
1904 : 0 : end subroutine vcgen_free
1905 : : !!***
1906 : :
1907 : : !----------------------------------------------------------------------
1908 : :
1909 : : !!****f* m_vcoul/vcgen_print
1910 : : !! NAME
1911 : : !! vcgen_print
1912 : : !!
1913 : : !! FUNCTION
1914 : : !! Print info on the object.
1915 : : !!
1916 : : !! SOURCE
1917 : :
1918 : 0 : subroutine vcgen_print(vcgen, units, header, prtvol)
1919 : :
1920 : : !Arguments ------------------------------------
1921 : : class(vcgen_t),intent(in) :: vcgen
1922 : : integer,intent(in) :: units(:), prtvol
1923 : : character(len=*),intent(in) :: header
1924 : :
1925 : : !Local variables-------------------------------
1926 : : !scalars
1927 : : integer :: ii
1928 : : character(len=500) :: msg
1929 : : ! *************************************************************************
1930 : :
1931 : : ABI_UNUSED(prtvol)
1932 : :
1933 [ # # # # ]: 0 : msg = ch10//' === '//trim(adjustl(header))//' === '
1934 : 0 : call wrtout(units, msg)
1935 : :
1936 : 0 : select case (vcgen%mode)
1937 : :
1938 : : case ('MINIBZ')
1939 : 0 : write(msg,'(3a)')ch10,' cutoff-mode = ',trim(vcgen%mode)
1940 : 0 : call wrtout(units, msg)
1941 : :
1942 : : case ('MINIBZ-ERF')
1943 : 0 : write(msg,'(3a)')ch10,' cutoff-mode = ',trim(vcgen%mode)
1944 : 0 : call wrtout(units, msg)
1945 : 0 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1946 : 0 : ' === Error function cutoff === ',ch10,ch10,&
1947 : 0 : ' Cutoff radius ......... ',vcgen%rcut,' [Bohr] ',ch10
1948 : 0 : call wrtout(units, msg)
1949 : :
1950 : : case ('MINIBZ-ERFC')
1951 : 0 : write(msg,'(3a)')ch10,' cutoff-mode = ',trim(vcgen%mode)
1952 : 0 : call wrtout(units, msg)
1953 : 0 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1954 : 0 : ' === Complement Error function cutoff === ',ch10,ch10,&
1955 : 0 : ' Cutoff radius ......... ',vcgen%rcut,' [Bohr] ',ch10
1956 : 0 : call wrtout(units, msg)
1957 : :
1958 : : case ('SPHERE')
1959 : 0 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1960 : 0 : ' === Spherical cutoff === ',ch10,ch10,&
1961 : 0 : ' Cutoff radius ......... ',vcgen%rcut,' [Bohr] ',ch10,&
1962 : 0 : ' Volume of the sphere .. ',four_pi/three*vcgen%rcut**3,' [Bohr^3] '
1963 : : !FB: This has no meaning here! & ' Sphere centered at .... ',vcgen%boxcenter,' (r.l.u) ',ch10
1964 : : !MG It might be useful if the system is not centered on the origin because in this case the
1965 : : ! matrix elements of the Coulomb have to be multiplied by a phase depending on boxcenter.
1966 : : ! I still have to decide if it is useful to code this possibility and which variable use to
1967 : : ! define the center (boxcenter is used in the tddft part).
1968 : 0 : call wrtout(units, msg)
1969 : :
1970 : : case ('CYLINDER')
1971 [ # # ]: 0 : ii=imin_loc(ABS(vcgen%pdir-1))
1972 : 0 : write(msg,'(5a,f10.4,3a,i2,2a,3f10.2,a)')ch10,&
1973 : 0 : ' === Cylindrical cutoff === ',ch10,ch10,&
1974 : 0 : ' Cutoff radius ............... ',vcgen%rcut,' [Bohr] ',ch10,&
1975 : 0 : ' Axis parallel to direction... ',ii,ch10,&
1976 : 0 : ' Passing through point ....... ',vcgen%boxcenter,' (r.l.u) '
1977 : 0 : call wrtout(units, msg)
1978 : :
1979 : 0 : write(msg,'(2a)')' Infinite length ....... ',ch10
1980 [ # # ]: 0 : if (vcgen%hcyl/=zero) write(msg,'(a,f8.5,2a)')' Finite length of ....... ',vcgen%hcyl,' [Bohr] ',ch10
1981 : 0 : call wrtout(units, msg)
1982 : :
1983 : : CASE ('SLAB')
1984 : 0 : write(msg,'(5a,f10.4,3a,3f10.2,2a)')ch10,&
1985 : 0 : ' === Surface cutoff === ',ch10,ch10,&
1986 : 0 : ' Cutoff radius .................... ',vcgen%rcut,' [Bohr] ',ch10,&
1987 : 0 : ' Central plane passing through .... ',vcgen%boxcenter,' (r.l.u) ',ch10
1988 : 0 : call wrtout(units, msg)
1989 : : !write(msg,'(a)')' Infinite length .......'
1990 : : !if (vcgen%hcyl/=zero) write(msg,'(a,f8.5,a)')' Finite length of .......',vcgen%hcyl,' [Bohr] '
1991 : : !call wrtout(units, msg)
1992 : :
1993 : : case ('AUXILIARY_FUNCTION')
1994 : 0 : write(msg,'(3a)')ch10,' cutoff-mode = ',trim(vcgen%mode)
1995 : 0 : call wrtout(units, msg)
1996 : :
1997 : : case ('AUX_GB')
1998 : 0 : write(msg,'(3a)')ch10,' cutoff-mode = ',trim(vcgen%mode)
1999 : 0 : call wrtout(units, msg)
2000 : :
2001 : : case ('CRYSTAL')
2002 : 0 : write(msg,'(3a)')ch10,' cutoff-mode = ',trim(vcgen%mode)
2003 : 0 : call wrtout(units, msg)
2004 : :
2005 : : case ('ERF')
2006 : 0 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
2007 : 0 : ' === Error function cutoff === ',ch10,ch10,&
2008 : 0 : ' Cutoff radius ......... ',vcgen%rcut,' [Bohr] ',ch10
2009 : 0 : call wrtout(units, msg)
2010 : :
2011 : : case ('ERFC')
2012 : 0 : write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
2013 : 0 : ' === Complement Error function cutoff === ',ch10,ch10,&
2014 : 0 : ' Cutoff radius ......... ',vcgen%rcut,' [Bohr] ',ch10
2015 : 0 : call wrtout(units, msg)
2016 : :
2017 : : case default
2018 [ # # # # : 0 : ABI_BUG(sjoin('Unknown cutoff mode: ', vcgen%mode))
# # # # #
# # # ]
2019 : : end select
2020 : :
2021 : : !if (vcgen%nqlwl > 0) then
2022 : : ! write(msg,'(a,i3)')" q-points for optical limit: ",vcgen%nqlwl
2023 : : ! call wrtout(units, msg)
2024 : : ! do iqlwl=1,vcgen%nqlwl
2025 : : ! write(msg,'(1x,i5,a,2x,3f12.6)') iqlwl,')',vcgen%qlwl(:,iqlwl)
2026 : : ! call wrtout(units, msg)
2027 : : ! end do
2028 : : !end if
2029 : :
2030 [ # # ]: 0 : end subroutine vcgen_print
2031 : : !!***
2032 : :
2033 [ - + - - : 1218 : end module m_vcoul
+ + - + +
- - + - +
+ - - + -
+ - + ]
2034 : : !!***
|
