Branch data Line data Source code
1 : : !!****m* ABINIT/m_double_grid
2 : : !! NAME
3 : : !! m_double_grid
4 : : !!
5 : : !! FUNCTION
6 : : !! This module defines the double grid object. This object contains the coarse mesh
7 : : !! and the dense mesh used for the interpolation of the BSE Hamiltonian,
8 : : !! and contains the mapping between the two meshes.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2008-2026 ABINIT group (YG, SP, MJV)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : MODULE m_double_grid
25 : :
26 : : use defs_basis
27 : : use m_errors
28 : : use m_abicore
29 : : use m_hide_blas
30 : : use m_bz_mesh
31 : : use m_krank
32 : :
33 : : use m_numeric_tools, only : wrap2_zero_one, interpol3d_indices
34 : : use m_matrix, only : matr3inv
35 : :
36 : : implicit none
37 : :
38 : : private
39 : : !!***
40 : :
41 : : !!****t* m_double_grid/double_grid_t
42 : : !! NAME
43 : : !! double_grid_t
44 : : !!
45 : : !! FUNCTION
46 : : !! The double grid contains a coarse mesh and a dense mesh
47 : : !! It also contains the mapping between the two meshes
48 : : !!
49 : : !! SOURCE
50 : :
51 : : type,public :: double_grid_t
52 : :
53 : : integer :: kmult(3)
54 : : ! Number of subdivisions in the coarse box in each direction
55 : :
56 : : integer :: ndiv
57 : : ! Total number of small sub-boxes in the large box associated to the coarse k-mesh.
58 : : ! i.e. product(kmult)
59 : :
60 : : integer :: maxcomp_coarse(3)
61 : : ! Dimensions of the box containing the coarse points in integer coord
62 : :
63 : : integer :: nbz_coarse
64 : : ! Number of k-points in the coarse BZ (open mesh)
65 : :
66 : : integer :: nbz_closedcoarse
67 : : ! Number of k-points inside the coarse BZ (closed mesh)
68 : : ! = PROD(maxcomp_coarse+1)
69 : :
70 : : integer :: nbz_dense
71 : : ! Number of k-point inside the dense BZ (open mesh)
72 : : ! = PROD(maxcomp_coarse.*kmult)
73 : :
74 : : integer, allocatable :: inttoik_coarse(:)
75 : : ! inttoik_coarse(nbz_closedcoarse)
76 : : ! Index of the kpoint in the coarse BZ.
77 : :
78 : : integer, allocatable :: iktoint_coarse(:)
79 : : ! iktoint_coarse(nbz_coarse)
80 : : ! Index of int by kpoint
81 : :
82 : : integer, allocatable :: inttoik_dense(:)
83 : :
84 : : integer, allocatable :: iktoint_dense(:)
85 : :
86 : : integer,allocatable :: indices_coarse(:,:)
87 : : ! indices_coarse(3,nbz_closedcoarse)
88 : : ! Indices (i1,i2,i3) for each point, ordinated by the integer coord
89 : :
90 : : integer,allocatable :: indices_dense(:,:)
91 : : ! indices_dense(6,nbz_dense)
92 : : ! Indices (i1,i2,i3);(j1,j2,j3) for each point, ordinated by integer coord
93 : :
94 : : integer :: kptrlatt_dense(3,3)
95 : : ! kptrlatt of the dense mesh
96 : :
97 : : real(dp) :: klatt_dense(3,3)
98 : :
99 : : integer :: nshiftk_dense
100 : : ! Number of shifts in the dense mesh.
101 : :
102 : : real(dp),allocatable :: shiftk_dense(:,:)
103 : : ! shift_dense(3,nshiftk_dense)
104 : : ! Shifts of the dense mesh.
105 : :
106 : : ! Dense lattice
107 : : integer :: kptrlatt_coarse(3,3)
108 : : ! kptrlatt of the coarse mesh
109 : :
110 : : real(dp) :: klatt_coarse(3,3)
111 : :
112 : : integer :: nshiftk_coarse
113 : : ! Number of shifts in the coarse mesh.
114 : :
115 : : real(dp),allocatable :: shiftk_coarse(:,:)
116 : : ! shift_coarse(3,nshiftk_coarse)
117 : : ! Shifts of the coarse mesh.
118 : :
119 : : ! Coarse lattice
120 : : integer, allocatable :: g0_coarse(:,:)
121 : : integer, allocatable :: g0_dense(:,:)
122 : : ! g0_dense/coarse(3,nkpt_closedcoarse/dense)
123 : : ! G0 vector between the kpt obtained with indices
124 : : ! and the kpt obtained insize bz
125 : :
126 : : integer, allocatable :: dense_to_coarse(:)
127 : : ! dense_to_coarse(nbz_dense)
128 : : ! Give the ibz_coarse corresponding to the dense mesh (the (0,0,0) point)
129 : :
130 : : integer, allocatable :: coarse_to_dense(:,:)
131 : : ! coarse_to_dense(nbz_coarse,ndiv)
132 : : ! Give all the ibz_dense corresponding to the (0,0,0) coarse point
133 : :
134 : : end type double_grid_t
135 : :
136 : : public :: double_grid_init ! Initializes the double grid with coarse mesh and dense mesh read from file
137 : : public :: double_grid_free ! Deallocate all memory
138 : : public :: get_kpt_from_indices_coarse ! Returns the k-point index and g0 vector associated to the set of indices
139 : : public :: compute_corresp ! Compute correspondance data between k-dense and k-coarse
140 : : !public :: get_kpt_from_indices_dense
141 : :
142 : : public :: kptfine_av ! Find the k-points of a fine grid that are around a k-point of a coarse mesh.
143 : : public :: k_neighbors ! Find 8 neighbors of given k-point on a coarse grid, and return
144 : : !!***
145 : :
146 : : !----------------------------------------------------------------------
147 : :
148 : : CONTAINS !=============================================================================
149 : : !!***
150 : :
151 : : !!****f* m_double_grid/double_grid_init
152 : : !! NAME
153 : : !! double_grid_init
154 : : !!
155 : : !! FUNCTION
156 : : !! Initialize the double_grid datatype "grid" from coarse and dense mesh
157 : : !!
158 : : !! INPUTS
159 : : !! Kmesh_coarse = descriptor of the coarse BZ sampling
160 : : !! Kmesh_dense = descriptor of the dense BZ sampling
161 : : !! kptrlatt_coarse(3,3) = vectors in R space that defines the reciprocal cell
162 : : !! kmult(3) = multiplication factors from coarse to dense
163 : : !!
164 : : !! OUTPUT
165 : : !! grid = double_grid to be created
166 : : !!
167 : : !! SOURCE
168 : :
169 [ - + - + : 4 : subroutine double_grid_init(Kmesh_coarse, Kmesh_dense, kptrlatt_coarse, kmult, grid)
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
170 : :
171 : : !Argument ------------------------------------
172 : : !scalars
173 : : type(double_grid_t),intent(out) :: grid
174 : : type(kmesh_t),intent(in) :: Kmesh_coarse,Kmesh_dense
175 : : !arrays
176 : : integer,intent(in) :: kptrlatt_coarse(3,3),kmult(3)
177 : :
178 : : !Local variables -----------------------------
179 : : !scalars
180 : : integer :: ii, info
181 : : !arrays
182 : : integer :: ipiv(3)
183 : : real(dp) :: rlatt_coarse(3,3),klatt_coarse(3,3),curmat(3,3)
184 : : !*********************************************
185 : :
186 [ - + ]: 4 : ABI_CHECK(Kmesh_coarse%nshift == 1, "Coarse mesh works only with nshiftk=1")
187 [ - + ]: 4 : ABI_CHECK(Kmesh_dense%nshift == 1, "Dense mesh : Works only with nshiftk=1")
188 : :
189 : 4 : grid%nshiftk_coarse = Kmesh_coarse%nshift
190 : 4 : grid%nshiftk_dense = Kmesh_dense%nshift
191 : :
192 : 4 : grid%nbz_coarse = Kmesh_coarse%nbz
193 : :
194 [ - + + - : 12 : ABI_MALLOC(grid%shiftk_coarse,(3,grid%nshiftk_coarse))
+ - + - -
+ - + -
+ ]
195 [ - + + - : 12 : ABI_MALLOC(grid%shiftk_dense,(3,grid%nshiftk_dense))
+ - + - -
+ - + -
+ ]
196 : :
197 [ + + + + ]: 20 : grid%shiftk_coarse(:,:) = Kmesh_coarse%shift(:,:)
198 [ + + + + ]: 20 : grid%shiftk_dense(:,:) = Kmesh_dense%shift(:,:)
199 : :
200 [ + + + + ]: 52 : grid%kptrlatt_coarse(:,:) = kptrlatt_coarse(:,:)
201 [ + + + + ]: 52 : rlatt_coarse(:,:) = kptrlatt_coarse(:,:)
202 : 4 : call matr3inv(rlatt_coarse,klatt_coarse)
203 [ + + + + ]: 52 : grid%klatt_coarse(:,:) = klatt_coarse(:,:)
204 : :
205 : 4 : grid%nbz_dense = Kmesh_dense%nbz
206 : 4 : grid%nbz_coarse = Kmesh_coarse%nbz
207 : :
208 [ + + ]: 16 : grid%kmult(:) = kmult(:)
209 : 4 : grid%ndiv = kmult(1)*kmult(2)*kmult(3)
210 : :
211 [ - + + - : 12 : ABI_MALLOC(grid%indices_dense,(6,Kmesh_dense%nbz))
+ - + - -
+ - + -
+ ]
212 [ - + + - : 12 : ABI_MALLOC(grid%g0_dense,(3,Kmesh_dense%nbz))
+ - + - -
+ - + -
+ ]
213 [ - + + - : 12 : ABI_MALLOC(grid%iktoint_dense,(Kmesh_dense%nbz))
+ - + - -
+ - + -
+ ]
214 [ - + + - : 8 : ABI_MALLOC(grid%inttoik_dense,(Kmesh_dense%nbz))
- + - + -
+ ]
215 : :
216 [ + + ]: 16 : grid%maxcomp_coarse(:) = -1
217 : :
218 [ + + + + ]: 52 : curmat(:,:) = grid%kptrlatt_coarse(:,:)
219 : :
220 : : ! Gaussian elimination
221 : 4 : call dgetrf(3,3,curmat,3,ipiv,info)
222 : :
223 : 4 : grid%nbz_closedcoarse = 1
224 : :
225 [ + + ]: 16 : do ii = 1,3
226 : 12 : grid%maxcomp_coarse(ii) = ABS(NINT(curmat(ipiv(ii),ipiv(ii))))
227 : 16 : grid%nbz_closedcoarse = grid%nbz_closedcoarse*(grid%maxcomp_coarse(ii)+1)
228 : : end do
229 : :
230 [ - + + - : 12 : ABI_MALLOC(grid%indices_coarse,(3,grid%nbz_closedcoarse))
+ - + - -
+ - + -
+ ]
231 [ - + + - : 8 : ABI_MALLOC(grid%g0_coarse,(3,grid%nbz_closedcoarse))
- + - + -
+ ]
232 [ - + + - : 12 : ABI_MALLOC(grid%iktoint_coarse,(Kmesh_coarse%nbz))
+ - + - -
+ - + -
+ ]
233 [ - + + - : 12 : ABI_MALLOC(grid%inttoik_coarse,(grid%nbz_closedcoarse))
+ - + - -
+ - + -
+ ]
234 : :
235 : : ! We should pass 'grid' at this stage !
236 : :
237 : : call create_indices_coarse(Kmesh_coarse%bz, Kmesh_coarse%nbz, grid%klatt_coarse, &
238 : : & grid%nshiftk_coarse, grid%shiftk_coarse, grid%maxcomp_coarse, grid%nbz_closedcoarse, grid%indices_coarse, &
239 : 4 : & grid%g0_coarse,grid%iktoint_coarse,grid%inttoik_coarse)
240 : :
241 : : call create_indices_dense(grid%klatt_coarse, grid%maxcomp_coarse, Kmesh_dense%bz, Kmesh_dense%nbz, &
242 : : & grid%nshiftk_dense, grid%shiftk_dense, grid%kmult, grid%indices_dense, grid%g0_dense, grid%iktoint_dense, &
243 : 4 : & grid%inttoik_dense)
244 : :
245 [ - + + - : 8 : ABI_MALLOC(grid%dense_to_coarse,(Kmesh_dense%nbz))
- + - + -
+ ]
246 [ - + + - : 16 : ABI_MALLOC(grid%coarse_to_dense,(Kmesh_coarse%nbz,grid%ndiv))
- + + - +
- + - - +
- + - + ]
247 : :
248 : : call compute_neighbours(grid%nbz_dense, grid%iktoint_dense, grid%indices_dense, &
249 : : & grid%maxcomp_coarse, grid%inttoik_coarse, grid%g0_coarse, grid%nbz_closedcoarse, grid%nbz_coarse,&
250 : 4 : & grid%ndiv, grid%dense_to_coarse, grid%coarse_to_dense)
251 : :
252 : 4 : end subroutine double_grid_init
253 : : !!***
254 : :
255 : : !----------------------------------------------------------------------------
256 : :
257 : : !!****f* m_double_grid/create_indices_coarse
258 : : !! NAME
259 : : !! create_indices_coarse
260 : : !!
261 : : !! FUNCTION
262 : : !! Create mapping between kpoints and integer indexing
263 : : !!
264 : : !! INPUTS
265 : : !! bz(3,nbz) = k-points in the Brillouin Zone
266 : : !! nbz = number of k-points
267 : : !! klatt(3,3) = reciprocal space vectors defining the reciprocal cell
268 : : !! nshiftk = Number of shifts
269 : : !! shiftk(3,nshiftk) = Shiftks of the Brillouin Zone
270 : : !! maxcomp(3) = Maximum int along each direction
271 : : !! nbz_closed = Number of k-points inside the closed Brillouin Zone (adding periodic images)
272 : : !!
273 : : !! OUTPUT
274 : : !! indices(3,nbz_closed) = indices for each k-point in the closed BZ
275 : : !! g0(3,nbz_closed) = g vectors between k-point inside bz and k-point given by indices
276 : : !! iktoint(nbz) = mapping between k-points in the bz and int indices
277 : : !! inttoik(nbz_closed) = mapping between int indices and k-points in the bz
278 : : !!
279 : : !! SOURCE
280 : :
281 : 4 : subroutine create_indices_coarse(bz, nbz, klatt, nshiftk, shiftk, maxcomp, nbz_closed, indices, g0, iktoint, inttoik)
282 : :
283 : : !Argument ------------------------------------
284 : : !scalars
285 : : integer,intent(in) :: nbz,nshiftk,nbz_closed
286 : : !arrays
287 : : integer,intent(in) :: maxcomp(3)
288 : : integer,intent(out) :: indices(3,nbz_closed)
289 : : integer,intent(out) :: g0(3,nbz_closed)
290 : : integer,intent(out) :: iktoint(nbz), inttoik(nbz_closed)
291 : : real(dp),intent(in) :: bz(3,nbz),klatt(3,3),shiftk(3,nshiftk)
292 : :
293 : : !Local variables -----------------------------
294 : : !scalars
295 : : integer :: ik,ii,i1,i2, i3
296 : : logical :: found
297 : : !arrays
298 : : integer :: curg0(3)
299 : : real(dp) :: curk1(3),ktoget(3)
300 : : !*********************************************
301 : :
302 [ - + ]: 4 : ABI_CHECK(nshiftk==1,"nshiftk != 1 not supported")
303 : :
304 [ + + ]: 16 : do i1 = 0,maxcomp(1)
305 [ + + ]: 52 : do i2 = 0,maxcomp(2)
306 [ + + ]: 156 : do i3 = 0,maxcomp(3)
307 : 108 : ii = (i1*(maxcomp(2)+1)+i2)*(maxcomp(3)+1)+i3+1
308 [ + + ]: 432 : ktoget(:) = shiftk(:,1)+(/i1,i2,i3/)
309 [ + + + + ]: 1404 : curk1(:) = MATMUL(klatt(:,:),ktoget(:))
310 : 360 : found = .FALSE.
311 [ + - ]: 360 : do ik = 1,nbz
312 [ + + ]: 360 : if(isamek(curk1(:),bz(:,ik),curg0)) then
313 [ + + ]: 432 : indices(:,ii) = (/i1,i2,i3/)
314 [ + + ]: 432 : g0(:,ii) = curg0
315 [ + + + + : 108 : if (i1 /= maxcomp(1) .and. i2 /= maxcomp(2) .and. i3 /= maxcomp(3)) then
+ + ]
316 : 32 : iktoint(ik) = ii
317 : : end if
318 : 108 : inttoik(ii) = ik
319 : 108 : found = .TRUE.
320 : : exit
321 : : end if
322 : : end do
323 : 36 : if (.not. found) then
324 : 0 : write(std_out,*) "curk1 = ",curk1
325 : 0 : write(std_out,*) bz
326 : 0 : ABI_ERROR("A k-point generated from kptrlatt cannot be found in the BZ")
327 : : end if
328 : : end do
329 : : end do
330 : : end do
331 : :
332 : 4 : end subroutine create_indices_coarse
333 : : !!***
334 : :
335 : : !----------------------------------------------------------------------
336 : :
337 : : !!****f* m_double_grid/get_kpt_from_indices_coarse
338 : : !! NAME
339 : : !! get_kpt_from_indices_coarse
340 : : !!
341 : : !! FUNCTION
342 : : !! Returns the k-point index and g0 vector associated to the set of indices
343 : : !!
344 : : !! INPUTS
345 : : !! indices(3) = index of the searched k-point
346 : : !! maxcomp(3) = Maximum int along each direction
347 : : !! inttoik(nkpt) = mapping between int indices and k-points in the bz
348 : : !! allg0(3,nkpt) = g vectors between k-point inside bz and k-point given by indices
349 : : !! nkpt = number of k-points
350 : : !!
351 : : !! OUTPUT
352 : : !! ikpt = index of k-point we search
353 : : !! g0(3) = g-vector obtained
354 : : !!
355 : : !! SOURCE
356 : :
357 : 21056 : subroutine get_kpt_from_indices_coarse(indices,maxcomp,inttoik,allg0,nkpt,ikpt,g0)
358 : :
359 : : !Argument ------------------------------------
360 : : !scalars
361 : : integer,intent(in) :: nkpt
362 : : integer,intent(out) :: ikpt
363 : : !arrays
364 : : integer,intent(in) :: indices(3),maxcomp(3)
365 : : integer,intent(in) :: inttoik(nkpt),allg0(3,nkpt)
366 : : integer,intent(out) :: g0(3)
367 : :
368 : : !Local variables -----------------------------
369 : : !scalars
370 : : integer :: curicoord
371 : : !*********************************************
372 : :
373 : 21056 : curicoord = (indices(1)*(maxcomp(2)+1)+indices(2))*(maxcomp(3)+1)+indices(3)+1
374 : 21056 : ikpt = inttoik(curicoord)
375 [ + + ]: 84224 : g0 = allg0(:,curicoord)
376 : :
377 : 21056 : end subroutine get_kpt_from_indices_coarse
378 : : !!***
379 : :
380 : : !----------------------------------------------------------------------
381 : :
382 : : !!****f* m_double_grid/create_indices_dense
383 : : !! NAME
384 : : !! create_indices_dense
385 : : !!
386 : : !! FUNCTION
387 : : !! Create mapping between kpoints and integer indexing
388 : : !!
389 : : !! INPUTS
390 : : !! klatt_coarse(3,3) = reciprocal space vectors defining the reciprocal cell of coarse BZ
391 : : !! maxcomp(3) = Maximum int along each direction
392 : : !! bz_dense(3,nbz_dense) = k-points in the dense BZ
393 : : !! nbz_dense = number of k-points in the dense BZ
394 : : !! nshiftk = Number of shifts
395 : : !! shiftk(3,nshiftk) = Shiftks of the Brillouin Zone
396 : : !! kmult(3) = multiplication factors
397 : : !! nbz_coarse = number of k-points in the coarse BZ
398 : : !! kptrlatt_coarse(3,3) = real space vectors defining the reciprocal cell of coarse BZ
399 : : !!
400 : : !! OUTPUT
401 : : !! indices(6,nbz_dense) = indices for each k-point in the closed BZ
402 : : !! g0(3,nbz_dense) = g vectors between k-point inside bz and k-point given by indices
403 : : !! iktoint(nbz_dense) = mapping between k-points in the bz and int indices
404 : : !! inttoik(nbz_dense) = mapping between int indices and k-points in the bz
405 : : !!
406 : : !! SOURCE
407 : :
408 : 4 : subroutine create_indices_dense(klatt_coarse, maxcomp, &
409 : 4 : & bz_dense, nbz_dense, nshiftk, shiftk, kmult, indices, g0, inttoik, iktoint)
410 : :
411 : : !Argument ------------------------------------
412 : : !scalars
413 : : integer,intent(in) :: nbz_dense, nshiftk
414 : : !arrays
415 : : integer,intent(in) :: kmult(3),maxcomp(3)
416 : : integer,intent(out) :: indices(6,nbz_dense),g0(3,nbz_dense)
417 : : integer,intent(out) :: inttoik(nbz_dense),iktoint(nbz_dense)
418 : : real(dp),intent(in) :: bz_dense(3,nbz_dense),klatt_coarse(3,3)
419 : : real(dp),intent(in) :: shiftk(3,nshiftk)
420 : :
421 : : !Local variables -----------------------------
422 : : integer :: ik,ii,ii_coarse
423 : : integer :: i1,i2,i3,j1,j2,j3
424 : : logical :: found
425 : : !arrays
426 : : integer :: curg0(3)
427 : : real(dp) :: curk1(3),ktoget(3)
428 : : !*********************************************
429 : :
430 : 4 : call wrtout(std_out, "Create Indices Dense", "COLL")
431 : :
432 [ - + ]: 4 : ABI_CHECK(nshiftk==1,"nshiftk != 1 not supported")
433 : :
434 [ + + ]: 12 : do i1 = 0,maxcomp(1)-1
435 [ + + ]: 28 : do i2 = 0,maxcomp(2)-1
436 [ + + ]: 56 : do i3 = 0,maxcomp(3)-1
437 : 32 : ii_coarse = (i1*(maxcomp(2)+1)+i2)*(maxcomp(3)+1)+i3+1
438 : :
439 [ + + ]: 112 : do j1 = 0,kmult(1)-1
440 [ + + ]: 224 : do j2 = 0,kmult(2)-1
441 [ + + ]: 448 : do j3 = 0,kmult(3)-1
442 : : ii = ((i1*kmult(1)+j1)*(maxcomp(2)*kmult(2)) +&
443 : : & (i2*kmult(2)+j2))*(maxcomp(3)*kmult(3))+&
444 : 256 : & (i3*kmult(3)+j3)+1
445 : :
446 : 256 : ktoget(1) = i1+((REAL(j1)+shiftk(1,1))/kmult(1))
447 : 256 : ktoget(2) = i2+((REAL(j2)+shiftk(2,1))/kmult(2))
448 : 256 : ktoget(3) = i3+((REAL(j3)+shiftk(3,1))/kmult(3))
449 : :
450 [ + + + + ]: 3328 : curk1(:) = MATMUL(klatt_coarse(:,:),ktoget(:))
451 : 8320 : found = .FALSE.
452 [ + - ]: 8320 : do ik = 1,nbz_dense
453 [ + + ]: 8320 : if(isamek(curk1(:),bz_dense(:,ik),curg0)) then
454 [ + + ]: 1792 : indices(:,ii) = (/i1,i2,i3,j1,j2,j3/)
455 [ + + ]: 1024 : g0(:,ii) = curg0
456 : 256 : inttoik(ii) = ik
457 : 256 : iktoint(ik) = ii
458 : 256 : found = .TRUE.
459 : : exit
460 : : end if
461 : : end do
462 : 128 : if(.not. found) then
463 : 0 : write(std_out,*) "curk1 = ",curk1
464 : 0 : write(std_out,*) bz_dense
465 : 0 : ABI_ERROR("Problem when creating indices")
466 : : end if
467 : : end do
468 : : end do
469 : : end do
470 : :
471 : : end do
472 : : end do
473 : : end do
474 : :
475 : 4 : end subroutine create_indices_dense
476 : : !!***
477 : :
478 : : !----------------------------------------------------------------------
479 : :
480 : : !!****f* m_double_grid/get_kpt_from_indices_dense
481 : : !! NAME
482 : : !! get_kpt_from_indices_coarse
483 : : !!
484 : : !! FUNCTION
485 : : !! Returns the k-point index and g0 vector associated to the set of indices
486 : : !!
487 : : !! INPUTS
488 : : !! indices(6) = index of the searched k-point
489 : : !! maxcomp(3) = Maximum int along each direction
490 : : !! kmult(3) = multiplication factors
491 : : !! inttoik(nkpt) = mapping between int indices and k-points in the bz
492 : : !! allg0(3,nkpt) = g vectors between k-point inside bz and k-point given by indices
493 : : !! nkpt = number of k-points
494 : : !!
495 : : !! OUTPUT
496 : : !! ikpt = index of k-point we search
497 : : !! g0(3) = g-vector obtained
498 : : !!
499 : : !! SOURCE
500 : :
501 : : subroutine get_kpt_from_indices_dense(indices,maxcomp,kmult,inttoik,allg0,nkpt,ikpt,g0)
502 : :
503 : : !Argument ------------------------------------
504 : : !scalars
505 : : integer, intent(in) :: nkpt
506 : : integer, intent(out) :: ikpt
507 : : !arrays
508 : : integer, intent(in) :: indices(6),maxcomp(3),inttoik(nkpt)
509 : : integer, intent(in) :: allg0(3,nkpt),kmult(3)
510 : : integer, intent(out) :: g0(3)
511 : :
512 : : !Local variables -----------------------------
513 : : !scalars
514 : : integer :: curicoord
515 : : !*********************************************
516 : :
517 : : curicoord = ((indices(1)*kmult(1)+indices(4))*(maxcomp(2)*kmult(2))+&
518 : : (indices(2)*kmult(2)+indices(5)))*(maxcomp(3)*kmult(3))+&
519 : : (indices(3)*kmult(3)+indices(6))+1
520 : :
521 : : ikpt = inttoik(curicoord)
522 : : g0 = allg0(:,curicoord)
523 : :
524 : : end subroutine get_kpt_from_indices_dense
525 : : !!***
526 : :
527 : : !----------------------------------------------------------------------
528 : :
529 : : !!****f* m_double_grid/compute_neighbours
530 : : !! NAME
531 : : !! compute_neighbours
532 : : !!
533 : : !! FUNCTION
534 : : !! Compute correspondance between points in the dense BZ and in the coarse BZ
535 : : !!
536 : : !! INPUTS
537 : : !! nbz_dense, nbz_closedcoarse, nbz_coarse, ndiv
538 : : !! iktoint_dense(nbz_dense)
539 : : !! indices_dense(6,nbz_dense)
540 : : !! maxcomp_coarse(3)
541 : : !! inttoik_coarse(nbz_closedcoarse)
542 : : !! g0_coarse(3,nbz_closedcoarse)
543 : : !!
544 : : !! OUTPUT
545 : : !! dense_to_coarse(nbz_dense)
546 : : !! coarse_to_dense(nbz_coarse,ndiv)
547 : : !!
548 : : !! SOURCE
549 : :
550 : 4 : subroutine compute_neighbours(nbz_dense, iktoint_dense, indices_dense, maxcomp_coarse, &
551 : 4 : & inttoik_coarse, g0_coarse, nbz_closedcoarse, nbz_coarse, ndiv, dense_to_coarse, coarse_to_dense)
552 : :
553 : : !Argument ------------------------------------
554 : : !scalars
555 : : integer,intent(in) :: nbz_dense, nbz_closedcoarse, nbz_coarse, ndiv
556 : : !arrays
557 : : integer,intent(in) :: iktoint_dense(nbz_dense)
558 : : integer,intent(in) :: indices_dense(6,nbz_dense)
559 : : integer,intent(in) :: maxcomp_coarse(3)
560 : : integer,intent(in) :: inttoik_coarse(nbz_closedcoarse)
561 : : integer,intent(in) :: g0_coarse(3,nbz_closedcoarse)
562 : : integer,intent(out) :: dense_to_coarse(nbz_dense)
563 : : integer,intent(out) :: coarse_to_dense(nbz_coarse,ndiv)
564 : :
565 : : !Local variables -----------------------------
566 : : !scalars
567 : : integer :: ik_dense, iorder, ik_coarse
568 : : !arrays
569 : 4 : integer :: curindex(nbz_coarse)
570 : : integer :: curindices_dense(6), curindices_coarse(3)
571 : : integer :: g0(3)
572 : : !*********************************************
573 : :
574 : : DBG_ENTER("COLL")
575 : :
576 [ + + + + ]: 292 : coarse_to_dense = 1
577 [ + + ]: 260 : dense_to_coarse = 1
578 : :
579 [ + + ]: 36 : curindex = 1
580 [ + + ]: 260 : do ik_dense = 1, nbz_dense
581 : : ! From ik_ibz in the dense mesh -> indices_dense
582 : 256 : iorder = iktoint_dense(ik_dense)
583 : :
584 : : ! From indices_dense -> indices_coarse
585 [ + + ]: 1792 : curindices_dense = indices_dense(:,iorder)
586 [ + + ]: 1024 : curindices_coarse = curindices_dense(1:3)
587 : : ! From indices_coarse -> ik_ibz in the coarse mesh
588 : : call get_kpt_from_indices_coarse(curindices_coarse,maxcomp_coarse,&
589 : 256 : & inttoik_coarse,g0_coarse,nbz_closedcoarse,ik_coarse,g0)
590 : :
591 : 256 : dense_to_coarse(ik_dense) = ik_coarse
592 : 256 : coarse_to_dense(ik_coarse, curindex(ik_coarse)) = ik_dense
593 : :
594 : 260 : curindex(ik_coarse) = curindex(ik_coarse) + 1
595 : : end do
596 : :
597 : : DBG_EXIT("COLL")
598 : :
599 : 4 : end subroutine compute_neighbours
600 : : !!***
601 : :
602 : : !---------------------------------------------------------------------
603 : :
604 : : !!****f* m_double_grid/compute_corresp
605 : : !! NAME
606 : : !! compute_corresp
607 : : !!
608 : : !! FUNCTION
609 : : !! Pre-process tables with mapping between divisions and coarse k-points
610 : : !!
611 : : !! INPUTS
612 : : !! double_grid
613 : : !!
614 : : !! OUTPUT
615 : : !! div2kdense(double_grid%nbz_coarse,double_grid%ndiv)
616 : : !! (k_coarse,idiv) -> k_dense
617 : : !! kdense2div(double_grid%nbz_dense)
618 : : !! k_dense -> idiv
619 : : !!
620 : : !! SOURCE
621 : :
622 : 4 : subroutine compute_corresp(double_grid, div2kdense, kdense2div)
623 : :
624 : : !Argument ------------------------------------
625 : : !scalars
626 : : class(double_grid_t),intent(in) :: double_grid
627 : : !arrays
628 : : integer,intent(out) :: div2kdense(double_grid%nbz_coarse,double_grid%ndiv)
629 : : integer,intent(out) :: kdense2div(double_grid%nbz_dense)
630 : :
631 : : !Local variables -----------------------------
632 : : !scalars
633 : : integer :: iorder,ik_dense,ik_coarse
634 : : !arrays
635 : : integer :: curindices_dense(6)
636 : 4 : integer,allocatable :: curindex(:)
637 : : !*********************************************
638 : :
639 [ - + + - : 12 : ABI_MALLOC(curindex,(double_grid%nbz_coarse))
+ - + - -
+ + - ]
640 [ + + ]: 36 : curindex = 1
641 : :
642 [ + + ]: 260 : do ik_dense = 1,double_grid%nbz_dense
643 : : ! From ik_ibz in the dense mesh -> indices_dense
644 : 256 : iorder = double_grid%iktoint_dense(ik_dense)
645 : : !g01 = double_grid%g0_dense(:,iorder)
646 : :
647 : : ! From indices_dense -> indices_coarse
648 : : curindices_dense = double_grid%indices_dense(:,iorder)
649 : :
650 : 256 : ik_coarse = double_grid%dense_to_coarse(ik_dense)
651 : 256 : div2kdense(ik_coarse,curindex(ik_coarse)) = ik_dense
652 : 256 : kdense2div(ik_dense) = curindex(ik_coarse)
653 : :
654 : 260 : curindex(ik_coarse) = curindex(ik_coarse) + 1
655 : : end do
656 : :
657 [ - + ]: 4 : ABI_FREE(curindex)
658 : :
659 : 4 : end subroutine compute_corresp
660 : : !!***
661 : :
662 : : !----------------------------------------------------------------------
663 : :
664 : : !!****f* m_double_grid/double_grid_free
665 : : !! NAME
666 : : !! double_grid_free
667 : : !!
668 : : !! FUNCTION
669 : : !! Deallocate all dynamics entities present in a double_grid structure.
670 : : !!
671 : : !! SOURCE
672 : :
673 : 4 : subroutine double_grid_free(grid)
674 : :
675 : : !Arguments ------------------------------------
676 : : class(double_grid_t),intent(inout) :: grid
677 : : ! *********************************************************************
678 : :
679 : : !integer
680 [ + - ]: 4 : ABI_SFREE(grid%inttoik_coarse)
681 [ + - ]: 4 : ABI_SFREE(grid%inttoik_dense)
682 [ + - ]: 4 : ABI_SFREE(grid%iktoint_coarse)
683 [ + - ]: 4 : ABI_SFREE(grid%iktoint_dense)
684 [ + - ]: 4 : ABI_SFREE(grid%indices_coarse)
685 [ + - ]: 4 : ABI_SFREE(grid%indices_dense)
686 [ + - ]: 4 : ABI_SFREE(grid%g0_coarse)
687 [ + - ]: 4 : ABI_SFREE(grid%g0_dense)
688 [ + - ]: 4 : ABI_SFREE(grid%dense_to_coarse)
689 [ + - ]: 4 : ABI_SFREE(grid%coarse_to_dense)
690 : :
691 : : !real
692 [ + - ]: 4 : ABI_SFREE(grid%shiftk_dense)
693 [ + - ]: 4 : ABI_SFREE(grid%shiftk_coarse)
694 : :
695 : 4 : end subroutine double_grid_free
696 : : !!***
697 : :
698 : : !----------------------------------------------------------------------
699 : :
700 : : !!****f* m_double_grid/kptfine_av
701 : : !! NAME
702 : : !! kptfine_av
703 : : !!
704 : : !! FUNCTION
705 : : !! Find the k-points of a fine grid that are around a k-point of a coarse mesh.
706 : : !!
707 : : !! INPUTS
708 : : !! center(3) = the point of the coarse mesh around which you want know which
709 : : !! k-points of the fine mesh belong to.
710 : : !! qptrlatt(3,3) = qptrlatt of the considered calculation (this is obtained
711 : : !! from the input variable ngqpt and shiftq.
712 : : !! kpt_fine(3,nkpt_fine) = this table contain all the k-points of the fine grid
713 : : !! in the full BZ (no sym op. allowed) and is read from the header
714 : : !! of the dense WF file.
715 : : !! nkpt_fine = number of k-points of the fine grid read from the header of the
716 : : !! dense WF file.
717 : : !!
718 : : !! OUTPUT
719 : : !! kpt_fine_sub(nkpt_sub) = k-points of the fine grid that are around center(3)
720 : : !! nkpt_sub = number of k-points of the fine grid that are around center(3)
721 : : !! wgt_sub(nkpt_sub) = weight of the k-points of the fine grid that are around center(3).
722 : : !!
723 : : !! SOURCE
724 : :
725 : 24 : subroutine kptfine_av(center,qptrlatt,kpt_fine,nkpt_fine,kpt_fine_sub,nkpt_sub,wgt_sub)
726 : :
727 : : !Arguments ------------------------------------
728 : : !scalars
729 : : integer,intent(in) :: nkpt_fine
730 : : integer,intent(out) :: nkpt_sub
731 : : !arrays
732 : : integer,intent(in) :: qptrlatt(3,3)
733 : : real(dp),intent(in) :: kpt_fine(3,nkpt_fine)
734 : : real(dp),intent(in) :: center(3)
735 : : integer,pointer :: kpt_fine_sub(:)
736 : : real(dp),pointer :: wgt_sub(:)
737 : :
738 : : !Local variables-------------------------------
739 : : !scalars
740 : : integer :: ikpt,aa,bb,cc
741 : : integer :: ii,jj
742 : : !arrays
743 : : real(dp) :: center_ref(3)
744 : : real(dp) :: kpt_fine_ref(3)
745 : : real(dp) :: kpt_tmp(3),kpt_tmp2(3)
746 : 24 : integer,allocatable :: kpt_fine_sub_tmp(:)
747 : 24 : real(dp),allocatable :: wgt_sub_tmp(:)
748 : : logical :: found(3)
749 : : ! *************************************************************************
750 : :
751 [ - + + - : 72 : ABI_MALLOC(kpt_fine_sub_tmp,(nkpt_fine))
+ - + - -
+ - + ]
752 [ - + + - : 72 : ABI_MALLOC(wgt_sub_tmp,(nkpt_fine))
+ - + - -
+ + - ]
753 : :
754 : : !It is easier to work in real space using the qptrlatt matrices because in this
755 : : !referential any k-points sampling will be cast into an orthorhombic shape.
756 : : !In that space we can simply take all k-points of the fine grid that between
757 : : !center_ref-0.5 and center_ref+0.5
758 : :
759 [ + + + + : 600 : center_ref = MATMUL(qptrlatt,center)
+ + + + ]
760 : :
761 : : !When considering points center(3) that lying close or on a BZ edge we need to
762 : : !take the k-points of the fine grid taking into account unklamp vectors. This
763 : : !is done with the aa, bb and cc loops.
764 : :
765 : : ii = 1
766 [ + + ]: 1560 : do ikpt=1,nkpt_fine
767 [ + + ]: 6144 : kpt_tmp = kpt_fine(:,ikpt)
768 [ + + ]: 6168 : do aa=-1,1
769 : 4608 : kpt_tmp2(1) = kpt_tmp(1)+aa
770 [ + + ]: 19968 : do bb=-1,1
771 : 13824 : kpt_tmp2(2) = kpt_tmp(2)+bb
772 [ + + ]: 59904 : do cc=-1,1
773 : 41472 : kpt_tmp2(3) = kpt_tmp(3)+cc
774 [ + + + + : 1036800 : kpt_fine_ref = MATMUL(qptrlatt,kpt_tmp2)
+ + + + ]
775 [ + + + + ]: 41472 : if((kpt_fine_ref(1)>=center_ref(1)-0.5-tol8).and.&
776 : 13824 : & (kpt_fine_ref(1)<=center_ref(1)+0.5+tol8)) then
777 [ + + + + ]: 10368 : if((kpt_fine_ref(2)>=center_ref(2)-0.5-tol8).and.&
778 : : & (kpt_fine_ref(2)<=center_ref(2)+0.5+tol8)) then
779 [ + + + + ]: 2592 : if((kpt_fine_ref(3)>=center_ref(3)-0.5-tol8).and.&
780 : : & (kpt_fine_ref(3)<=center_ref(3)+0.5+tol8)) then
781 : 648 : kpt_fine_sub_tmp(ii) = ikpt
782 : 648 : ii = ii +1
783 : : end if
784 : : end if
785 : : end if
786 : : end do
787 : : end do
788 : : end do
789 : : end do
790 : :
791 : 24 : nkpt_sub = ii-1
792 [ - + + - : 72 : ABI_MALLOC(kpt_fine_sub,(nkpt_sub))
+ - + - -
+ - + ]
793 [ - + + - : 72 : ABI_MALLOC(wgt_sub,(nkpt_sub))
+ - + - -
+ - + ]
794 : :
795 [ + + ]: 672 : do jj=1,nkpt_sub
796 : 672 : kpt_fine_sub(jj) = kpt_fine_sub_tmp(jj)
797 : : end do
798 : :
799 : : !We then compute a weight function. This weight function is simply a
800 : : !rectangular weight function that take the value 1 for k-points of the fine
801 : : !grid inside the cube, 0.5 for k-points that are lying on one face of the cube,
802 : : !0.25 for k-points that are lying on an edge of the cube and 0.125 for k-points
803 : : !that are lying on a peak of the cube.
804 : :
805 [ + + ]: 672 : wgt_sub(:) = 1.0
806 : :
807 [ + + ]: 672 : do ikpt=1,nkpt_sub
808 [ + + ]: 2592 : found(:) = .True.
809 [ + + ]: 2592 : kpt_tmp = kpt_fine(:,kpt_fine_sub(ikpt))
810 [ + + ]: 2616 : do aa=-1,1
811 : 1944 : kpt_tmp2(1) = kpt_tmp(1)+aa
812 [ + + ]: 8424 : do bb=-1,1
813 : 5832 : kpt_tmp2(2) = kpt_tmp(2)+bb
814 [ + + ]: 25272 : do cc=-1,1
815 : 17496 : kpt_tmp2(3) = kpt_tmp(3)+cc
816 [ + + + + : 437400 : kpt_fine_ref = MATMUL(qptrlatt,kpt_tmp2)
+ + + + ]
817 [ + + + + ]: 17496 : if((ABS(kpt_fine_ref(1)-center_ref(1)-0.5)< tol8) .or.&
818 : : (ABS(kpt_fine_ref(1)-center_ref(1)+0.5) < tol8)) then
819 [ + + ]: 3888 : if(found(1)) then
820 : 432 : wgt_sub(ikpt) = wgt_sub(ikpt)*0.5; found(1) = .False.
821 : : end if
822 : : end if
823 [ + + + + ]: 17496 : if((ABS(kpt_fine_ref(2)-center_ref(2)-0.5) < tol8) .or.&
824 : : (ABS(kpt_fine_ref(2)-center_ref(2)+0.5) < tol8)) then
825 [ + + ]: 3888 : if(found(2)) then
826 : 432 : wgt_sub(ikpt) = wgt_sub(ikpt)*0.5; found(2) = .False.
827 : : end if
828 : : end if
829 [ + + + + ]: 17496 : if((ABS(kpt_fine_ref(3)-center_ref(3)-0.5)< tol8) .or.&
830 : 5832 : (ABS(kpt_fine_ref(3)-center_ref(3)+0.5) < tol8)) then
831 [ + + ]: 3888 : if(found(3)) then
832 : 432 : wgt_sub(ikpt) = wgt_sub(ikpt)*0.5; found(3) = .False.
833 : : end if
834 : : end if
835 : : end do
836 : : end do
837 : : end do
838 : : end do
839 : :
840 [ - + ]: 24 : ABI_FREE(kpt_fine_sub_tmp)
841 [ - + ]: 24 : ABI_FREE(wgt_sub_tmp)
842 : :
843 : 24 : end subroutine kptfine_av
844 : : !!***
845 : :
846 : : !!****f* m_double_grid/k_neighbors
847 : : !!
848 : : !! NAME
849 : : !! k_neighbors
850 : : !!
851 : : !! FUNCTION
852 : : !! find 8 neighbors of given k-point on a coarse grid, and return
853 : : !! them along with relative k-shift within coarse grid cell
854 : : !!
855 : : !! INPUTS
856 : : !! kpt = k-point to be interpolated to, in full BZ
857 : : !! kptrlatt = lattice vectors for coarse k-grid
858 : : !! invrankkpt = rank list to find k-points
859 : : !!
860 : : !! OUTPUT
861 : : !! rel_kpt = k-point coordinates renormalized to coarse grid cell
862 : : !! kpt_phon_indices = indices of k-points on corners of cell
863 : : !!
864 : : !! TODO
865 : : !! This routine is not used anymore. Deprecate or Remove?
866 : : !!
867 : : !! SOURCE
868 : :
869 : 0 : subroutine k_neighbors(kpt, kptrlatt,krank, rel_kpt, kpt_phon_indices)
870 : :
871 : : ! inputs
872 : : real(dp), intent(in) :: kpt(3)
873 : : integer, intent(in) :: kptrlatt(3,3)
874 : : type(krank_t), intent(in) :: krank
875 : :
876 : : ! outputs
877 : : real(dp), intent(out) :: rel_kpt(3)
878 : : integer, intent(out) :: kpt_phon_indices(8)
879 : : ! local vars
880 : : integer :: symrankkpt
881 : : integer :: ir1,ir2,ir3, pr1,pr2,pr3
882 : : real(dp) :: redkpt(3), cornerkpt(3), res
883 : : ! *************************************************************************
884 : :
885 : : !wrap fine kpt to [0,1]
886 : 0 : call wrap2_zero_one(kpt(1),redkpt(1),res)
887 : 0 : call wrap2_zero_one(kpt(2),redkpt(2),res)
888 : 0 : call wrap2_zero_one(kpt(3),redkpt(3),res)
889 : : !find 8 indices of points neighboring ikpt_phon, for interpolation
890 : : call interpol3d_indices (redkpt,kptrlatt(1,1),kptrlatt(2,2),kptrlatt(3,3), &
891 : 0 : ir1,ir2,ir3, pr1,pr2,pr3)
892 : :
893 : : !transpose ir pr to ikpt_phon indices
894 : : !order of kpt_phons:
895 : : !ir1 ir2 ir3
896 [ # # ]: 0 : cornerkpt = (/real(ir1-1)/kptrlatt(1,1),real(ir2-1)/kptrlatt(2,2), real(ir3-1)/kptrlatt(3,3)/)
897 : 0 : symrankkpt = krank%get_rank(cornerkpt)
898 : 0 : kpt_phon_indices(1) = krank%invrank(symrankkpt)
899 : : !pr1 ir2 ir3
900 [ # # ]: 0 : cornerkpt = (/real(pr1-1)/kptrlatt(1,1),real(ir2-1)/kptrlatt(2,2), real(ir3-1)/kptrlatt(3,3)/)
901 : 0 : symrankkpt = krank%get_rank (cornerkpt)
902 : 0 : kpt_phon_indices(2) = krank%invrank(symrankkpt)
903 : : !ir1 pr2 ir3
904 [ # # ]: 0 : cornerkpt = (/real(ir1-1)/kptrlatt(1,1),real(pr2-1)/kptrlatt(2,2), real(ir3-1)/kptrlatt(3,3)/)
905 : 0 : symrankkpt = krank%get_rank (cornerkpt)
906 : 0 : kpt_phon_indices(3) = krank%invrank(symrankkpt)
907 : : !pr1 pr2 ir3
908 [ # # ]: 0 : cornerkpt = (/real(pr1-1)/kptrlatt(1,1),real(pr2-1)/kptrlatt(2,2), real(ir3-1)/kptrlatt(3,3)/)
909 : 0 : symrankkpt = krank%get_rank (cornerkpt)
910 : 0 : kpt_phon_indices(4) = krank%invrank(symrankkpt)
911 : : !ir1 ir2 pr3
912 [ # # ]: 0 : cornerkpt = (/real(ir1-1)/kptrlatt(1,1),real(ir2-1)/kptrlatt(2,2), real(pr3-1)/kptrlatt(3,3)/)
913 : 0 : symrankkpt = krank%get_rank (cornerkpt)
914 : 0 : kpt_phon_indices(5) = krank%invrank(symrankkpt)
915 : : !pr1 ir2 pr3
916 [ # # ]: 0 : cornerkpt = (/real(pr1-1)/kptrlatt(1,1),real(ir2-1)/kptrlatt(2,2), real(pr3-1)/kptrlatt(3,3)/)
917 : 0 : symrankkpt = krank%get_rank (cornerkpt)
918 : 0 : kpt_phon_indices(6) = krank%invrank(symrankkpt)
919 : : !ir1 pr2 pr3
920 [ # # ]: 0 : cornerkpt = (/real(ir1-1)/kptrlatt(1,1),real(pr2-1)/kptrlatt(2,2), real(pr3-1)/kptrlatt(3,3)/)
921 : 0 : symrankkpt = krank%get_rank (cornerkpt)
922 : 0 : kpt_phon_indices(7) = krank%invrank(symrankkpt)
923 : : !pr1 pr2 pr3
924 [ # # ]: 0 : cornerkpt = (/real(pr1-1)/kptrlatt(1,1),real(pr2-1)/kptrlatt(2,2), real(pr3-1)/kptrlatt(3,3)/)
925 : 0 : symrankkpt = krank%get_rank (cornerkpt)
926 : 0 : kpt_phon_indices(8) = krank%invrank(symrankkpt)
927 : :
928 : : !retrieve the gkq matrix for all q, at the neighbor k vectors
929 : 0 : rel_kpt(1) = redkpt(1)*kptrlatt(1,1)-real(ir1-1)
930 : 0 : rel_kpt(2) = redkpt(2)*kptrlatt(2,2)-real(ir2-1)
931 : 0 : rel_kpt(3) = redkpt(3)*kptrlatt(3,3)-real(ir3-1)
932 : :
933 : 0 : end subroutine k_neighbors
934 : : !!***
935 : :
936 [ # # # # : 0 : END MODULE m_double_grid
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
937 : : !!***
|
