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