Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice
10 : !> summation in real space. These integrals are e.g. needed in periodic RI for RPA, GW
11 : !> \par History
12 : !> 2018.05 created [Jan Wilhelm]
13 : !> \author Jan Wilhelm
14 : ! **************************************************************************************************
15 : MODULE kpoint_coulomb_2c
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE basis_set_types, ONLY: gto_basis_set_type
19 : USE cell_types, ONLY: cell_type,&
20 : get_cell,&
21 : pbc
22 : USE constants_operator, ONLY: operator_coulomb
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_create, dbcsr_init_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
25 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
26 : dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
27 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
28 : USE generic_shg_integrals, ONLY: int_operators_r12_ab_shg
29 : USE generic_shg_integrals_init, ONLY: contraction_matrix_shg
30 : USE kinds, ONLY: dp
31 : USE kpoint_types, ONLY: get_kpoint_info,&
32 : kpoint_type
33 : USE mathconstants, ONLY: gaussi,&
34 : twopi,&
35 : z_one
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE particle_types, ONLY: particle_type
38 : USE qs_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : USE qs_kind_types, ONLY: get_qs_kind,&
41 : qs_kind_type
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_coulomb_2c'
49 :
50 : PUBLIC :: build_2c_coulomb_matrix_kp, build_2c_coulomb_matrix_kp_small_cell
51 :
52 : ! **************************************************************************************************
53 :
54 : TYPE two_d_util_type
55 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: block
56 : END TYPE
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief ...
62 : !> \param matrix_v_kp ...
63 : !> \param kpoints ...
64 : !> \param basis_type ...
65 : !> \param cell ...
66 : !> \param particle_set ...
67 : !> \param qs_kind_set ...
68 : !> \param atomic_kind_set ...
69 : !> \param size_lattice_sum ...
70 : !> \param operator_type ...
71 : !> \param ikp_start ...
72 : !> \param ikp_end ...
73 : ! **************************************************************************************************
74 100 : SUBROUTINE build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, &
75 : atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
76 :
77 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
78 : TYPE(kpoint_type), POINTER :: kpoints
79 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
80 : TYPE(cell_type), POINTER :: cell
81 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
82 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
83 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
84 : INTEGER :: size_lattice_sum, operator_type, &
85 : ikp_start, ikp_end
86 :
87 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_2c_coulomb_matrix_kp'
88 :
89 : INTEGER :: handle
90 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
91 :
92 100 : CALL timeset(routineN, handle)
93 :
94 100 : CALL allocate_tmp(matrix_v_L_tmp, matrix_v_kp, ikp_start)
95 :
96 : CALL lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, &
97 : qs_kind_set, atomic_kind_set, size_lattice_sum, matrix_v_L_tmp, &
98 100 : operator_type, ikp_start, ikp_end)
99 :
100 100 : CALL deallocate_tmp(matrix_v_L_tmp)
101 :
102 100 : CALL timestop(handle)
103 :
104 100 : END SUBROUTINE build_2c_coulomb_matrix_kp
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param matrix_v_kp ...
109 : !> \param kpoints ...
110 : !> \param basis_type ...
111 : !> \param cell ...
112 : !> \param particle_set ...
113 : !> \param qs_kind_set ...
114 : !> \param atomic_kind_set ...
115 : !> \param size_lattice_sum ...
116 : !> \param matrix_v_L_tmp ...
117 : !> \param operator_type ...
118 : !> \param ikp_start ...
119 : !> \param ikp_end ...
120 : ! **************************************************************************************************
121 100 : SUBROUTINE lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, &
122 : qs_kind_set, atomic_kind_set, size_lattice_sum, matrix_v_L_tmp, &
123 : operator_type, ikp_start, ikp_end)
124 :
125 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
126 : TYPE(kpoint_type), POINTER :: kpoints
127 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
128 : TYPE(cell_type), POINTER :: cell
129 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
130 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
131 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
132 : INTEGER :: size_lattice_sum
133 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
134 : INTEGER :: operator_type, ikp_start, ikp_end
135 :
136 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lattice_sum'
137 :
138 : INTEGER :: factor, handle, handle2, i_block, i_x, i_x_inner, i_x_outer, ik, j_y, j_y_inner, &
139 : j_y_outer, k_z, k_z_inner, k_z_outer, x_max, x_min, y_max, y_min, z_max, z_min
140 : INTEGER, DIMENSION(3) :: nkp_grid
141 : REAL(KIND=dp) :: coskl, sinkl
142 : REAL(KIND=dp), DIMENSION(3) :: vec_L, vec_s
143 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
144 100 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L, blocks_v_L_store
145 : TYPE(two_d_util_type), ALLOCATABLE, &
146 100 : DIMENSION(:, :, :) :: blocks_v_kp
147 :
148 100 : CALL timeset(routineN, handle)
149 :
150 : CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
151 100 : x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
152 :
153 100 : CALL allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
154 100 : CALL allocate_blocks_v_L(blocks_v_L, matrix_v_L_tmp)
155 100 : CALL allocate_blocks_v_L(blocks_v_L_store, matrix_v_L_tmp)
156 :
157 488 : DO i_x_inner = 0, 2*nkp_grid(1) - 1
158 3784 : DO j_y_inner = 0, 2*nkp_grid(2) - 1
159 24100 : DO k_z_inner = 0, 2*nkp_grid(3) - 1
160 :
161 50944 : DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
162 160000 : DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
163 547520 : DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
164 :
165 407936 : i_x = i_x_inner + i_x_outer
166 407936 : j_y = j_y_inner + j_y_outer
167 407936 : k_z = k_z_inner + k_z_outer
168 :
169 : IF (i_x > x_max .OR. i_x < x_min .OR. &
170 : j_y > y_max .OR. j_y < y_min .OR. &
171 407936 : k_z > z_max .OR. k_z < z_min) CYCLE
172 :
173 625408 : vec_s = [REAL(i_x, dp), REAL(j_y, dp), REAL(k_z, dp)]
174 :
175 2032576 : vec_L = MATMUL(hmat, vec_s)
176 :
177 : ! Compute (P 0 | Q vec_L) and store it in matrix_v_L_tmp
178 : CALL compute_v_transl(matrix_v_L_tmp, blocks_v_L, vec_L, particle_set, &
179 : qs_kind_set, atomic_kind_set, basis_type, cell, &
180 156352 : operator_type)
181 :
182 845092 : DO i_block = 1, SIZE(blocks_v_L)
183 : blocks_v_L_store(i_block)%block(:, :) = blocks_v_L_store(i_block)%block(:, :) &
184 139785176 : + blocks_v_L(i_block)%block(:, :)
185 : END DO
186 :
187 : END DO
188 : END DO
189 : END DO
190 :
191 20416 : CALL timeset(routineN//"_R_to_k", handle2)
192 :
193 : ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
194 169728 : DO ik = ikp_start, ikp_end
195 :
196 : ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
197 597248 : coskl = COS(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
198 597248 : sinkl = SIN(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
199 :
200 748352 : DO i_block = 1, SIZE(blocks_v_L)
201 :
202 : blocks_v_kp(ik, 1, i_block)%block(:, :) = blocks_v_kp(ik, 1, i_block)%block(:, :) &
203 282880576 : + coskl*blocks_v_L_store(i_block)%block(:, :)
204 : blocks_v_kp(ik, 2, i_block)%block(:, :) = blocks_v_kp(ik, 2, i_block)%block(:, :) &
205 283029888 : + sinkl*blocks_v_L_store(i_block)%block(:, :)
206 :
207 : END DO
208 :
209 : END DO
210 :
211 93088 : DO i_block = 1, SIZE(blocks_v_L)
212 :
213 19736480 : blocks_v_L_store(i_block)%block(:, :) = 0.0_dp
214 :
215 : END DO
216 :
217 44128 : CALL timestop(handle2)
218 :
219 : END DO
220 : END DO
221 : END DO
222 :
223 100 : CALL set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
224 :
225 100 : CALL deallocate_blocks_v_kp(blocks_v_kp)
226 100 : CALL deallocate_blocks_v_L(blocks_v_L)
227 100 : CALL deallocate_blocks_v_L(blocks_v_L_store)
228 :
229 100 : CALL timestop(handle)
230 :
231 100 : END SUBROUTINE lattice_sum
232 :
233 : ! **************************************************************************************************
234 : !> \brief ...
235 : !> \param matrix_v_kp ...
236 : !> \param blocks_v_kp ...
237 : !> \param ikp_start ...
238 : !> \param ikp_end ...
239 : ! **************************************************************************************************
240 100 : SUBROUTINE set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
241 :
242 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
243 : TYPE(two_d_util_type), ALLOCATABLE, &
244 : DIMENSION(:, :, :) :: blocks_v_kp
245 : INTEGER :: ikp_start, ikp_end
246 :
247 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_blocks_to_matrix_v_kp'
248 :
249 : INTEGER :: col, handle, i_block, i_real_im, ik, row
250 100 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
251 : TYPE(dbcsr_iterator_type) :: iter
252 :
253 100 : CALL timeset(routineN, handle)
254 :
255 740 : DO ik = ikp_start, ikp_end
256 :
257 2020 : DO i_real_im = 1, 2
258 :
259 1280 : i_block = 1
260 :
261 1280 : CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_im)%matrix)
262 :
263 6220 : DO WHILE (dbcsr_iterator_blocks_left(iter))
264 :
265 4940 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
266 :
267 2423116 : data_block(:, :) = blocks_v_kp(ik, i_real_im, i_block)%block(:, :)
268 :
269 4940 : i_block = i_block + 1
270 :
271 : END DO
272 :
273 3200 : CALL dbcsr_iterator_stop(iter)
274 :
275 : END DO
276 :
277 : END DO
278 :
279 100 : CALL timestop(handle)
280 :
281 100 : END SUBROUTINE set_blocks_to_matrix_v_kp
282 :
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param matrix_v_L_tmp ...
286 : !> \param blocks_v_L ...
287 : !> \param vec_L ...
288 : !> \param particle_set ...
289 : !> \param qs_kind_set ...
290 : !> \param atomic_kind_set ...
291 : !> \param basis_type ...
292 : !> \param cell ...
293 : !> \param operator_type ...
294 : ! **************************************************************************************************
295 156352 : SUBROUTINE compute_v_transl(matrix_v_L_tmp, blocks_v_L, vec_L, particle_set, &
296 : qs_kind_set, atomic_kind_set, basis_type, cell, operator_type)
297 :
298 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
299 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L
300 : REAL(KIND=dp), DIMENSION(3) :: vec_L
301 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
302 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
303 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
304 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
305 : TYPE(cell_type), POINTER :: cell
306 : INTEGER :: operator_type
307 :
308 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_transl'
309 :
310 : INTEGER :: col, handle, i_block, kind_a, kind_b, row
311 156352 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
312 : REAL(dp), DIMENSION(3) :: ra, rab_L, rb
313 156352 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
314 156352 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: contr_a, contr_b
315 : TYPE(dbcsr_iterator_type) :: iter
316 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
317 :
318 156352 : CALL timeset(routineN, handle)
319 :
320 156352 : NULLIFY (basis_set_a, basis_set_b, data_block)
321 :
322 156352 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
323 :
324 156352 : CALL dbcsr_set(matrix_v_L_tmp, 0.0_dp)
325 :
326 156352 : i_block = 1
327 :
328 156352 : CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
329 :
330 736036 : DO WHILE (dbcsr_iterator_blocks_left(iter))
331 :
332 579684 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
333 :
334 579684 : kind_a = kind_of(row)
335 579684 : kind_b = kind_of(col)
336 :
337 579684 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, basis_type=basis_type)
338 579684 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, basis_type=basis_type)
339 :
340 579684 : ra(1:3) = pbc(particle_set(row)%r(1:3), cell)
341 579684 : rb(1:3) = pbc(particle_set(col)%r(1:3), cell)
342 :
343 2318736 : rab_L(1:3) = rb(1:3) - ra(1:3) + vec_L(1:3)
344 :
345 579684 : CALL contraction_matrix_shg(basis_set_a, contr_a)
346 579684 : CALL contraction_matrix_shg(basis_set_b, contr_b)
347 :
348 139377240 : blocks_v_L(i_block)%block = 0.0_dp
349 :
350 : CALL int_operators_r12_ab_shg(operator_type, blocks_v_L(i_block)%block, rab=rab_L, &
351 : fba=basis_set_a, fbb=basis_set_b, scona_shg=contr_a, sconb_shg=contr_b, &
352 579684 : calculate_forces=.FALSE.)
353 :
354 579684 : i_block = i_block + 1
355 :
356 579684 : DEALLOCATE (contr_a, contr_b)
357 :
358 : END DO
359 :
360 156352 : CALL dbcsr_iterator_stop(iter)
361 :
362 156352 : DEALLOCATE (kind_of)
363 :
364 156352 : CALL timestop(handle)
365 :
366 156352 : END SUBROUTINE compute_v_transl
367 :
368 : ! **************************************************************************************************
369 : !> \brief ...
370 : !> \param blocks_v_kp ...
371 : ! **************************************************************************************************
372 100 : SUBROUTINE deallocate_blocks_v_kp(blocks_v_kp)
373 : TYPE(two_d_util_type), ALLOCATABLE, &
374 : DIMENSION(:, :, :) :: blocks_v_kp
375 :
376 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_blocks_v_kp'
377 :
378 : INTEGER :: handle, i_block, i_real_img, ik
379 :
380 100 : CALL timeset(routineN, handle)
381 :
382 940 : DO ik = LBOUND(blocks_v_kp, 1), UBOUND(blocks_v_kp, 1)
383 2020 : DO i_real_img = 1, SIZE(blocks_v_kp, 2)
384 6860 : DO i_block = 1, SIZE(blocks_v_kp, 3)
385 6220 : DEALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block)
386 : END DO
387 : END DO
388 : END DO
389 :
390 5040 : DEALLOCATE (blocks_v_kp)
391 :
392 100 : CALL timestop(handle)
393 :
394 100 : END SUBROUTINE
395 :
396 : ! **************************************************************************************************
397 : !> \brief ...
398 : !> \param blocks_v_L ...
399 : ! **************************************************************************************************
400 200 : SUBROUTINE deallocate_blocks_v_L(blocks_v_L)
401 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L
402 :
403 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_blocks_v_L'
404 :
405 : INTEGER :: handle, i_block
406 :
407 200 : CALL timeset(routineN, handle)
408 :
409 900 : DO i_block = 1, SIZE(blocks_v_L, 1)
410 900 : DEALLOCATE (blocks_v_L(i_block)%block)
411 : END DO
412 :
413 900 : DEALLOCATE (blocks_v_L)
414 :
415 200 : CALL timestop(handle)
416 :
417 200 : END SUBROUTINE
418 :
419 : ! **************************************************************************************************
420 : !> \brief ...
421 : !> \param blocks_v_L ...
422 : !> \param matrix_v_L_tmp ...
423 : ! **************************************************************************************************
424 200 : SUBROUTINE allocate_blocks_v_L(blocks_v_L, matrix_v_L_tmp)
425 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L
426 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
427 :
428 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_blocks_v_L'
429 :
430 : INTEGER :: col, handle, i_block, nblocks, row
431 200 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
432 : TYPE(dbcsr_iterator_type) :: iter
433 :
434 200 : CALL timeset(routineN, handle)
435 :
436 200 : nblocks = 0
437 :
438 200 : CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
439 :
440 900 : DO WHILE (dbcsr_iterator_blocks_left(iter))
441 :
442 700 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
443 :
444 700 : nblocks = nblocks + 1
445 :
446 : END DO
447 :
448 200 : CALL dbcsr_iterator_stop(iter)
449 :
450 1300 : ALLOCATE (blocks_v_L(nblocks))
451 :
452 200 : i_block = 1
453 :
454 200 : CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
455 :
456 900 : DO WHILE (dbcsr_iterator_blocks_left(iter))
457 :
458 700 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
459 :
460 2800 : ALLOCATE (blocks_v_L(i_block)%block(SIZE(data_block, 1), SIZE(data_block, 2)))
461 217232 : blocks_v_L(i_block)%block = 0.0_dp
462 :
463 1600 : i_block = i_block + 1
464 :
465 : END DO
466 :
467 200 : CALL dbcsr_iterator_stop(iter)
468 :
469 200 : CALL timestop(handle)
470 :
471 400 : END SUBROUTINE
472 :
473 : ! **************************************************************************************************
474 : !> \brief ...
475 : !> \param blocks_v_kp ...
476 : !> \param matrix_v_kp ...
477 : !> \param ikp_start ...
478 : !> \param ikp_end ...
479 : ! **************************************************************************************************
480 100 : SUBROUTINE allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
481 : TYPE(two_d_util_type), ALLOCATABLE, &
482 : DIMENSION(:, :, :) :: blocks_v_kp
483 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
484 : INTEGER :: ikp_start, ikp_end
485 :
486 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_blocks_v_kp'
487 :
488 : INTEGER :: col, handle, i_block, i_real_img, ik, &
489 : nblocks, row
490 100 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
491 : TYPE(dbcsr_iterator_type) :: iter
492 :
493 100 : CALL timeset(routineN, handle)
494 :
495 100 : nblocks = 0
496 :
497 100 : CALL dbcsr_iterator_start(iter, matrix_v_kp(ikp_start, 1)%matrix)
498 :
499 450 : DO WHILE (dbcsr_iterator_blocks_left(iter))
500 :
501 350 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
502 :
503 350 : nblocks = nblocks + 1
504 :
505 : END DO
506 :
507 100 : CALL dbcsr_iterator_stop(iter)
508 :
509 6490 : ALLOCATE (blocks_v_kp(ikp_start:ikp_end, SIZE(matrix_v_kp, 2), nblocks))
510 :
511 740 : DO ik = ikp_start, ikp_end
512 :
513 2020 : DO i_real_img = 1, SIZE(matrix_v_kp, 2)
514 :
515 1280 : i_block = 1
516 :
517 1280 : CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_img)%matrix)
518 :
519 6220 : DO WHILE (dbcsr_iterator_blocks_left(iter))
520 :
521 4940 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
522 :
523 0 : ALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block(SIZE(data_block, 1), &
524 19760 : SIZE(data_block, 2)))
525 2423116 : blocks_v_kp(ik, i_real_img, i_block)%block = 0.0_dp
526 :
527 11160 : i_block = i_block + 1
528 :
529 : END DO
530 :
531 3200 : CALL dbcsr_iterator_stop(iter)
532 :
533 : END DO
534 :
535 : END DO
536 :
537 100 : CALL timestop(handle)
538 :
539 100 : END SUBROUTINE
540 :
541 : ! **************************************************************************************************
542 : !> \brief ...
543 : !> \param cell ...
544 : !> \param kpoints ...
545 : !> \param size_lattice_sum ...
546 : !> \param factor ...
547 : !> \param hmat ...
548 : !> \param x_min ...
549 : !> \param x_max ...
550 : !> \param y_min ...
551 : !> \param y_max ...
552 : !> \param z_min ...
553 : !> \param z_max ...
554 : !> \param nkp_grid ...
555 : ! **************************************************************************************************
556 112 : SUBROUTINE get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
557 : x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
558 :
559 : TYPE(cell_type), POINTER :: cell
560 : TYPE(kpoint_type), POINTER :: kpoints
561 : INTEGER :: size_lattice_sum, factor
562 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
563 : INTEGER :: x_min, x_max, y_min, y_max, z_min, z_max
564 : INTEGER, DIMENSION(3) :: nkp_grid
565 :
566 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_factor_and_xyz_min_max'
567 :
568 : INTEGER :: handle, nkp
569 : INTEGER, DIMENSION(3) :: periodic
570 :
571 112 : CALL timeset(routineN, handle)
572 :
573 112 : CALL get_kpoint_info(kpoints, nkp_grid=nkp_grid, nkp=nkp)
574 112 : CALL get_cell(cell=cell, h=hmat, periodic=periodic)
575 :
576 112 : IF (periodic(1) == 0) THEN
577 86 : CPASSERT(nkp_grid(1) == 1)
578 : END IF
579 112 : IF (periodic(2) == 0) THEN
580 20 : CPASSERT(nkp_grid(2) == 1)
581 : END IF
582 112 : IF (periodic(3) == 0) THEN
583 26 : CPASSERT(nkp_grid(3) == 1)
584 : END IF
585 :
586 112 : IF (MODULO(nkp_grid(1), 2) == 1) THEN
587 86 : factor = 3**(size_lattice_sum - 1)
588 26 : ELSE IF (MODULO(nkp_grid(1), 2) == 0) THEN
589 26 : factor = 2**(size_lattice_sum - 1)
590 : END IF
591 :
592 112 : IF (MODULO(nkp_grid(1), 2) == 1) THEN
593 86 : x_min = -(factor*nkp_grid(1) - 1)/2
594 86 : x_max = (factor*nkp_grid(1) - 1)/2
595 26 : ELSE IF (MODULO(nkp_grid(1), 2) == 0) THEN
596 26 : x_min = -factor*nkp_grid(1)/2
597 26 : x_max = factor*nkp_grid(1)/2 - 1
598 : END IF
599 112 : IF (periodic(1) == 0) THEN
600 86 : x_min = 0
601 86 : x_max = 0
602 : END IF
603 :
604 112 : IF (MODULO(nkp_grid(2), 2) == 1) THEN
605 20 : y_min = -(factor*nkp_grid(2) - 1)/2
606 20 : y_max = (factor*nkp_grid(2) - 1)/2
607 92 : ELSE IF (MODULO(nkp_grid(2), 2) == 0) THEN
608 92 : y_min = -factor*nkp_grid(2)/2
609 92 : y_max = factor*nkp_grid(2)/2 - 1
610 : END IF
611 112 : IF (periodic(2) == 0) THEN
612 20 : y_min = 0
613 20 : y_max = 0
614 : END IF
615 :
616 112 : IF (MODULO(nkp_grid(3), 2) == 1) THEN
617 26 : z_min = -(factor*nkp_grid(3) - 1)/2
618 26 : z_max = (factor*nkp_grid(3) - 1)/2
619 86 : ELSE IF (MODULO(nkp_grid(3), 2) == 0) THEN
620 86 : z_min = -factor*nkp_grid(3)/2
621 86 : z_max = factor*nkp_grid(3)/2 - 1
622 : END IF
623 112 : IF (periodic(3) == 0) THEN
624 26 : z_min = 0
625 26 : z_max = 0
626 : END IF
627 :
628 112 : CALL timestop(handle)
629 :
630 112 : END SUBROUTINE
631 :
632 : ! **************************************************************************************************
633 : !> \brief ...
634 : !> \param matrix_v_L_tmp ...
635 : !> \param matrix_v_kp ...
636 : !> \param ikp_start ...
637 : ! **************************************************************************************************
638 100 : SUBROUTINE allocate_tmp(matrix_v_L_tmp, matrix_v_kp, ikp_start)
639 :
640 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
641 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
642 : INTEGER :: ikp_start
643 :
644 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_tmp'
645 :
646 : INTEGER :: handle
647 :
648 100 : CALL timeset(routineN, handle)
649 :
650 100 : NULLIFY (matrix_v_L_tmp)
651 100 : CALL dbcsr_init_p(matrix_v_L_tmp)
652 : CALL dbcsr_create(matrix=matrix_v_L_tmp, &
653 : template=matrix_v_kp(ikp_start, 1)%matrix, &
654 100 : matrix_type=dbcsr_type_no_symmetry)
655 100 : CALL dbcsr_reserve_all_blocks(matrix_v_L_tmp)
656 100 : CALL dbcsr_set(matrix_v_L_tmp, 0.0_dp)
657 :
658 100 : CALL timestop(handle)
659 :
660 100 : END SUBROUTINE
661 :
662 : ! **************************************************************************************************
663 : !> \brief ...
664 : !> \param matrix_v_L_tmp ...
665 : ! **************************************************************************************************
666 100 : SUBROUTINE deallocate_tmp(matrix_v_L_tmp)
667 :
668 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
669 :
670 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_tmp'
671 :
672 : INTEGER :: handle
673 :
674 100 : CALL timeset(routineN, handle)
675 :
676 100 : CALL dbcsr_release_p(matrix_v_L_tmp)
677 :
678 100 : CALL timestop(handle)
679 :
680 100 : END SUBROUTINE
681 :
682 : ! **************************************************************************************************
683 : !> \brief ...
684 : !> \param V_k ...
685 : !> \param qs_env ...
686 : !> \param kpoints ...
687 : !> \param size_lattice_sum ...
688 : !> \param basis_type ...
689 : !> \param ikp_start ...
690 : !> \param ikp_end ...
691 : ! **************************************************************************************************
692 12 : SUBROUTINE build_2c_coulomb_matrix_kp_small_cell(V_k, qs_env, kpoints, size_lattice_sum, &
693 : basis_type, ikp_start, ikp_end)
694 : COMPLEX(KIND=dp), DIMENSION(:, :, :) :: V_k
695 : TYPE(qs_environment_type), POINTER :: qs_env
696 : TYPE(kpoint_type), POINTER :: kpoints
697 : INTEGER :: size_lattice_sum
698 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
699 : INTEGER :: ikp_start, ikp_end
700 :
701 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_2c_coulomb_matrix_kp_small_cell'
702 :
703 : INTEGER :: factor, handle, handle2, i_cell, i_x, i_x_inner, i_x_outer, ik, ikp_local, j_y, &
704 : j_y_inner, j_y_outer, k_z, k_z_inner, k_z_outer, n_atom, n_bf, x_max, x_min, y_max, &
705 : y_min, z_max, z_min
706 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bf_end_from_atom, bf_start_from_atom
707 : INTEGER, DIMENSION(3) :: nkp_grid
708 : REAL(KIND=dp) :: coskl, sinkl
709 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: V_L
710 : REAL(KIND=dp), DIMENSION(3) :: vec_L, vec_s
711 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
712 12 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
713 : TYPE(cell_type), POINTER :: cell
714 : TYPE(mp_para_env_type), POINTER :: para_env
715 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
716 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
717 :
718 12 : CALL timeset(routineN, handle)
719 :
720 : CALL get_qs_env(qs_env=qs_env, &
721 : para_env=para_env, &
722 : particle_set=particle_set, &
723 : cell=cell, &
724 : qs_kind_set=qs_kind_set, &
725 12 : atomic_kind_set=atomic_kind_set)
726 :
727 : CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
728 12 : x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
729 :
730 12 : CALL get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
731 :
732 48 : ALLOCATE (V_L(n_bf, n_bf))
733 :
734 220 : DO i_x_inner = 0, 2*nkp_grid(1) - 1
735 11228 : DO j_y_inner = 0, 2*nkp_grid(2) - 1
736 72656 : DO k_z_inner = 0, 2*nkp_grid(3) - 1
737 :
738 2396160 : V_L(:, :) = 0.0_dp
739 61440 : i_cell = 0
740 :
741 163840 : DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
742 552960 : DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
743 1699840 : DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
744 :
745 1208320 : i_x = i_x_inner + i_x_outer
746 1208320 : j_y = j_y_inner + j_y_outer
747 1208320 : k_z = k_z_inner + k_z_outer
748 :
749 : IF (i_x > x_max .OR. i_x < x_min .OR. &
750 : j_y > y_max .OR. j_y < y_min .OR. &
751 1208320 : k_z > z_max .OR. k_z < z_min) CYCLE
752 :
753 455680 : i_cell = i_cell + 1
754 :
755 1822720 : vec_s = [REAL(i_x, dp), REAL(j_y, dp), REAL(k_z, dp)]
756 :
757 455680 : IF (MODULO(i_cell, para_env%num_pe) .NE. para_env%mepos) CYCLE
758 :
759 2961920 : vec_L = MATMUL(hmat, vec_s)
760 :
761 : ! Compute (P 0 | Q vec_L) and add it to V_R
762 : CALL add_V_L(V_L, vec_L, n_atom, bf_start_from_atom, bf_end_from_atom, &
763 1597440 : particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
764 :
765 : END DO
766 : END DO
767 : END DO
768 :
769 61440 : CALL para_env%sync()
770 61440 : CALL para_env%sum(V_L)
771 :
772 61440 : CALL timeset(routineN//"_R_to_k", handle2)
773 :
774 61440 : ikp_local = 0
775 :
776 : ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
777 33091584 : DO ik = 1, ikp_end
778 :
779 33030144 : IF (MODULO(ik, para_env%num_pe) .NE. para_env%mepos) CYCLE
780 :
781 16515072 : ikp_local = ikp_local + 1
782 :
783 16515072 : IF (ik < ikp_start) CYCLE
784 :
785 : ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
786 53477376 : coskl = COS(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
787 53477376 : sinkl = SIN(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
788 :
789 : V_k(:, :, ikp_local) = V_k(:, :, ikp_local) + z_one*coskl*V_L(:, :) + &
790 524611584 : gaussi*sinkl*V_L(:, :)
791 :
792 : END DO
793 :
794 133888 : CALL timestop(handle2)
795 :
796 : END DO
797 : END DO
798 : END DO
799 :
800 12 : CALL timestop(handle)
801 :
802 24 : END SUBROUTINE build_2c_coulomb_matrix_kp_small_cell
803 :
804 : ! **************************************************************************************************
805 : !> \brief ...
806 : !> \param qs_env ...
807 : !> \param n_atom ...
808 : !> \param basis_type ...
809 : !> \param bf_start_from_atom ...
810 : !> \param bf_end_from_atom ...
811 : !> \param n_bf ...
812 : ! **************************************************************************************************
813 12 : SUBROUTINE get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
814 :
815 : TYPE(qs_environment_type), POINTER :: qs_env
816 : INTEGER :: n_atom
817 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
818 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bf_start_from_atom, bf_end_from_atom
819 : INTEGER :: n_bf
820 :
821 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_basis_sizes'
822 :
823 : INTEGER :: handle, iatom, ikind, n_kind, nsgf
824 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
825 12 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
826 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
827 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
828 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
829 :
830 12 : CALL timeset(routineN, handle)
831 :
832 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
833 12 : qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
834 12 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
835 :
836 12 : n_atom = SIZE(particle_set)
837 12 : n_kind = SIZE(qs_kind_set)
838 :
839 36 : DO ikind = 1, n_kind
840 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
841 24 : basis_type=basis_type)
842 36 : CPASSERT(ASSOCIATED(basis_set_a))
843 : END DO
844 :
845 60 : ALLOCATE (bf_start_from_atom(n_atom), bf_end_from_atom(n_atom))
846 :
847 12 : n_bf = 0
848 44 : DO iatom = 1, n_atom
849 32 : bf_start_from_atom(iatom) = n_bf + 1
850 32 : ikind = kind_of(iatom)
851 32 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type=basis_type)
852 32 : n_bf = n_bf + nsgf
853 44 : bf_end_from_atom(iatom) = n_bf
854 : END DO
855 :
856 12 : CALL timestop(handle)
857 :
858 24 : END SUBROUTINE get_basis_sizes
859 :
860 : ! **************************************************************************************************
861 : !> \brief ...
862 : !> \param V_L ...
863 : !> \param vec_L ...
864 : !> \param n_atom ...
865 : !> \param bf_start_from_atom ...
866 : !> \param bf_end_from_atom ...
867 : !> \param particle_set ...
868 : !> \param qs_kind_set ...
869 : !> \param atomic_kind_set ...
870 : !> \param basis_type ...
871 : !> \param cell ...
872 : ! **************************************************************************************************
873 227840 : SUBROUTINE add_V_L(V_L, vec_L, n_atom, bf_start_from_atom, bf_end_from_atom, &
874 : particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
875 :
876 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: V_L
877 : REAL(KIND=dp), DIMENSION(3) :: vec_L
878 : INTEGER :: n_atom
879 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bf_start_from_atom, bf_end_from_atom
880 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
881 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
882 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
883 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
884 : TYPE(cell_type), POINTER :: cell
885 :
886 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_V_L'
887 :
888 : INTEGER :: a_1, a_2, atom_a, atom_b, b_1, b_2, &
889 : handle, kind_a, kind_b
890 227840 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
891 : REAL(dp), DIMENSION(3) :: ra, rab_L, rb
892 227840 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: V_L_ab
893 227840 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: contr_a, contr_b
894 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
895 :
896 227840 : CALL timeset(routineN, handle)
897 :
898 227840 : NULLIFY (basis_set_a, basis_set_b)
899 :
900 227840 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
901 :
902 890880 : DO atom_a = 1, n_atom
903 :
904 2839040 : DO atom_b = 1, n_atom
905 :
906 1948160 : kind_a = kind_of(atom_a)
907 1948160 : kind_b = kind_of(atom_b)
908 :
909 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, &
910 1948160 : basis_type=basis_type)
911 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, &
912 1948160 : basis_type=basis_type)
913 :
914 1948160 : ra(1:3) = pbc(particle_set(atom_a)%r(1:3), cell)
915 1948160 : rb(1:3) = pbc(particle_set(atom_b)%r(1:3), cell)
916 :
917 7792640 : rab_L(1:3) = rb(1:3) - ra(1:3) + vec_L(1:3)
918 :
919 1948160 : CALL contraction_matrix_shg(basis_set_a, contr_a)
920 1948160 : CALL contraction_matrix_shg(basis_set_b, contr_b)
921 :
922 1948160 : a_1 = bf_start_from_atom(atom_a)
923 1948160 : a_2 = bf_end_from_atom(atom_a)
924 1948160 : b_1 = bf_start_from_atom(atom_b)
925 1948160 : b_2 = bf_end_from_atom(atom_b)
926 :
927 7792640 : ALLOCATE (V_L_ab(a_2 - a_1 + 1, b_2 - b_1 + 1))
928 :
929 : CALL int_operators_r12_ab_shg(operator_coulomb, V_L_ab, rab=rab_L, &
930 : fba=basis_set_a, fbb=basis_set_b, &
931 : scona_shg=contr_a, sconb_shg=contr_b, &
932 1948160 : calculate_forces=.FALSE.)
933 :
934 13862400 : V_L(a_1:a_2, b_1:b_2) = V_L(a_1:a_2, b_1:b_2) + V_L_ab(:, :)
935 :
936 2611200 : DEALLOCATE (contr_a, contr_b, V_L_ab)
937 :
938 : END DO
939 :
940 : END DO
941 :
942 227840 : DEALLOCATE (kind_of)
943 :
944 227840 : CALL timestop(handle)
945 :
946 227840 : END SUBROUTINE add_V_L
947 :
948 0 : END MODULE kpoint_coulomb_2c
|
