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