Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 two_d_util_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 104 : 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 104 : CALL timeset(routineN, handle)
93 :
94 104 : 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 104 : operator_type, ikp_start, ikp_end)
99 :
100 104 : CALL deallocate_tmp(matrix_v_L_tmp)
101 :
102 104 : CALL timestop(handle)
103 :
104 104 : 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 104 : 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 104 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L, blocks_v_L_store
145 : TYPE(two_d_util_type), ALLOCATABLE, &
146 104 : DIMENSION(:, :, :) :: blocks_v_kp
147 :
148 104 : CALL timeset(routineN, handle)
149 :
150 : CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
151 104 : x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
152 :
153 104 : CALL allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
154 104 : CALL allocate_blocks_v_L(blocks_v_L, matrix_v_L_tmp)
155 104 : CALL allocate_blocks_v_L(blocks_v_L_store, matrix_v_L_tmp)
156 :
157 500 : DO i_x_inner = 0, 2*nkp_grid(1) - 1
158 3812 : DO j_y_inner = 0, 2*nkp_grid(2) - 1
159 24204 : DO k_z_inner = 0, 2*nkp_grid(3) - 1
160 :
161 51104 : DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
162 160240 : DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
163 548016 : DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
164 :
165 408272 : i_x = i_x_inner + i_x_outer
166 408272 : j_y = j_y_inner + j_y_outer
167 408272 : 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 408272 : k_z > z_max .OR. k_z < z_min) CYCLE
172 :
173 625704 : vec_s = [REAL(i_x, dp), REAL(j_y, dp), REAL(k_z, dp)]
174 :
175 2033538 : 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 156426 : operator_type)
181 :
182 : !$OMP PARALLEL DO DEFAULT(NONE) &
183 : !$OMP SHARED(blocks_v_L, blocks_v_L_store) &
184 517408 : !$OMP PRIVATE(i_block)
185 : DO i_block = 1, SIZE(blocks_v_L)
186 : blocks_v_L_store(i_block)%block(:, :) = blocks_v_L_store(i_block)%block(:, :) &
187 : + blocks_v_L(i_block)%block(:, :)
188 : END DO
189 : !$OMP END PARALLEL DO
190 :
191 : END DO
192 : END DO
193 : END DO
194 :
195 20496 : CALL timeset(routineN//"_R_to_k", handle2)
196 :
197 : ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
198 169952 : DO ik = ikp_start, ikp_end
199 :
200 : ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
201 597824 : coskl = COS(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
202 597824 : sinkl = SIN(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
203 :
204 749184 : DO i_block = 1, SIZE(blocks_v_L)
205 :
206 : blocks_v_kp(ik, 1, i_block)%block(:, :) = blocks_v_kp(ik, 1, i_block)%block(:, :) &
207 283342240 : + coskl*blocks_v_L_store(i_block)%block(:, :)
208 : blocks_v_kp(ik, 2, i_block)%block(:, :) = blocks_v_kp(ik, 2, i_block)%block(:, :) &
209 283491696 : + sinkl*blocks_v_L_store(i_block)%block(:, :)
210 :
211 : END DO
212 :
213 : END DO
214 :
215 93488 : DO i_block = 1, SIZE(blocks_v_L)
216 :
217 19969520 : blocks_v_L_store(i_block)%block(:, :) = 0.0_dp
218 :
219 : END DO
220 :
221 44304 : CALL timestop(handle2)
222 :
223 : END DO
224 : END DO
225 : END DO
226 :
227 104 : CALL set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
228 :
229 104 : CALL deallocate_blocks_v_kp(blocks_v_kp)
230 104 : CALL deallocate_blocks_v_L(blocks_v_L)
231 104 : CALL deallocate_blocks_v_L(blocks_v_L_store)
232 :
233 104 : CALL timestop(handle)
234 :
235 104 : END SUBROUTINE lattice_sum
236 :
237 : ! **************************************************************************************************
238 : !> \brief ...
239 : !> \param matrix_v_kp ...
240 : !> \param blocks_v_kp ...
241 : !> \param ikp_start ...
242 : !> \param ikp_end ...
243 : ! **************************************************************************************************
244 104 : SUBROUTINE set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
245 :
246 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
247 : TYPE(two_d_util_type), ALLOCATABLE, &
248 : DIMENSION(:, :, :) :: blocks_v_kp
249 : INTEGER :: ikp_start, ikp_end
250 :
251 : CHARACTER(LEN=*), PARAMETER :: routineN = 'set_blocks_to_matrix_v_kp'
252 :
253 : INTEGER :: col, handle, i_block, i_real_im, ik, row
254 104 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
255 : TYPE(dbcsr_iterator_type) :: iter
256 :
257 104 : CALL timeset(routineN, handle)
258 :
259 750 : DO ik = ikp_start, ikp_end
260 :
261 2042 : DO i_real_im = 1, 2
262 :
263 1292 : i_block = 1
264 :
265 1292 : CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_im)%matrix)
266 :
267 6276 : DO WHILE (dbcsr_iterator_blocks_left(iter))
268 :
269 4984 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
270 :
271 2452768 : data_block(:, :) = blocks_v_kp(ik, i_real_im, i_block)%block(:, :)
272 :
273 4984 : i_block = i_block + 1
274 :
275 : END DO
276 :
277 3230 : CALL dbcsr_iterator_stop(iter)
278 :
279 : END DO
280 :
281 : END DO
282 :
283 104 : CALL timestop(handle)
284 :
285 104 : END SUBROUTINE set_blocks_to_matrix_v_kp
286 :
287 : ! **************************************************************************************************
288 : !> \brief ...
289 : !> \param matrix_v_L_tmp ...
290 : !> \param blocks_v_L ...
291 : !> \param vec_L ...
292 : !> \param particle_set ...
293 : !> \param qs_kind_set ...
294 : !> \param atomic_kind_set ...
295 : !> \param basis_type ...
296 : !> \param cell ...
297 : !> \param operator_type ...
298 : ! **************************************************************************************************
299 156426 : SUBROUTINE compute_v_transl(matrix_v_L_tmp, blocks_v_L, vec_L, particle_set, &
300 : qs_kind_set, atomic_kind_set, basis_type, cell, operator_type)
301 :
302 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
303 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L
304 : REAL(KIND=dp), DIMENSION(3) :: vec_L
305 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
306 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
307 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
308 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
309 : TYPE(cell_type), POINTER :: cell
310 : INTEGER :: operator_type
311 :
312 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_transl'
313 :
314 : INTEGER :: col, handle, i_block, kind_a, kind_b, row
315 156426 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
316 : REAL(dp), DIMENSION(3) :: ra, rab_L, rb
317 156426 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
318 156426 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: contr_a, contr_b
319 : TYPE(dbcsr_iterator_type) :: iter
320 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
321 :
322 156426 : CALL timeset(routineN, handle)
323 :
324 156426 : NULLIFY (basis_set_a, basis_set_b, data_block)
325 :
326 156426 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
327 :
328 156426 : CALL dbcsr_set(matrix_v_L_tmp, 0.0_dp)
329 :
330 156426 : i_block = 1
331 :
332 156426 : CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
333 :
334 736438 : DO WHILE (dbcsr_iterator_blocks_left(iter))
335 :
336 580012 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
337 :
338 580012 : kind_a = kind_of(row)
339 580012 : kind_b = kind_of(col)
340 :
341 580012 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, basis_type=basis_type)
342 580012 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, basis_type=basis_type)
343 :
344 580012 : ra(1:3) = pbc(particle_set(row)%r(1:3), cell)
345 580012 : rb(1:3) = pbc(particle_set(col)%r(1:3), cell)
346 :
347 2320048 : rab_L(1:3) = rb(1:3) - ra(1:3) + vec_L(1:3)
348 :
349 580012 : CALL contraction_matrix_shg(basis_set_a, contr_a)
350 580012 : CALL contraction_matrix_shg(basis_set_b, contr_b)
351 :
352 139635064 : blocks_v_L(i_block)%block = 0.0_dp
353 :
354 : CALL int_operators_r12_ab_shg(operator_type, blocks_v_L(i_block)%block, rab=rab_L, &
355 : fba=basis_set_a, fbb=basis_set_b, scona_shg=contr_a, sconb_shg=contr_b, &
356 580012 : calculate_forces=.FALSE.)
357 :
358 580012 : i_block = i_block + 1
359 :
360 580012 : DEALLOCATE (contr_a, contr_b)
361 :
362 : END DO
363 :
364 156426 : CALL dbcsr_iterator_stop(iter)
365 :
366 156426 : DEALLOCATE (kind_of)
367 :
368 156426 : CALL timestop(handle)
369 :
370 156426 : END SUBROUTINE compute_v_transl
371 :
372 : ! **************************************************************************************************
373 : !> \brief ...
374 : !> \param blocks_v_kp ...
375 : ! **************************************************************************************************
376 104 : SUBROUTINE deallocate_blocks_v_kp(blocks_v_kp)
377 : TYPE(two_d_util_type), ALLOCATABLE, &
378 : DIMENSION(:, :, :) :: blocks_v_kp
379 :
380 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_blocks_v_kp'
381 :
382 : INTEGER :: handle, i_block, i_real_img, ik
383 :
384 104 : CALL timeset(routineN, handle)
385 :
386 958 : DO ik = LBOUND(blocks_v_kp, 1), UBOUND(blocks_v_kp, 1)
387 2042 : DO i_real_img = 1, SIZE(blocks_v_kp, 2)
388 6922 : DO i_block = 1, SIZE(blocks_v_kp, 3)
389 6276 : DEALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block)
390 : END DO
391 : END DO
392 : END DO
393 :
394 5088 : DEALLOCATE (blocks_v_kp)
395 :
396 104 : CALL timestop(handle)
397 :
398 104 : END SUBROUTINE deallocate_blocks_v_kp
399 :
400 : ! **************************************************************************************************
401 : !> \brief ...
402 : !> \param blocks_v_L ...
403 : ! **************************************************************************************************
404 208 : SUBROUTINE deallocate_blocks_v_L(blocks_v_L)
405 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L
406 :
407 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_blocks_v_L'
408 :
409 : INTEGER :: handle, i_block
410 :
411 208 : CALL timeset(routineN, handle)
412 :
413 934 : DO i_block = 1, SIZE(blocks_v_L, 1)
414 934 : DEALLOCATE (blocks_v_L(i_block)%block)
415 : END DO
416 :
417 934 : DEALLOCATE (blocks_v_L)
418 :
419 208 : CALL timestop(handle)
420 :
421 208 : END SUBROUTINE deallocate_blocks_v_L
422 :
423 : ! **************************************************************************************************
424 : !> \brief ...
425 : !> \param blocks_v_L ...
426 : !> \param matrix_v_L_tmp ...
427 : ! **************************************************************************************************
428 208 : SUBROUTINE allocate_blocks_v_L(blocks_v_L, matrix_v_L_tmp)
429 : TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_L
430 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
431 :
432 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_blocks_v_L'
433 :
434 : INTEGER :: col, handle, i_block, nblocks, row
435 208 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
436 : TYPE(dbcsr_iterator_type) :: iter
437 :
438 208 : CALL timeset(routineN, handle)
439 :
440 208 : nblocks = 0
441 :
442 208 : CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
443 :
444 934 : DO WHILE (dbcsr_iterator_blocks_left(iter))
445 :
446 726 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
447 :
448 726 : nblocks = nblocks + 1
449 :
450 : END DO
451 :
452 208 : CALL dbcsr_iterator_stop(iter)
453 :
454 1350 : ALLOCATE (blocks_v_L(nblocks))
455 :
456 208 : i_block = 1
457 :
458 208 : CALL dbcsr_iterator_start(iter, matrix_v_L_tmp)
459 :
460 934 : DO WHILE (dbcsr_iterator_blocks_left(iter))
461 :
462 726 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
463 :
464 2904 : ALLOCATE (blocks_v_L(i_block)%block(SIZE(data_block, 1), SIZE(data_block, 2)))
465 232590 : blocks_v_L(i_block)%block = 0.0_dp
466 :
467 1660 : i_block = i_block + 1
468 :
469 : END DO
470 :
471 208 : CALL dbcsr_iterator_stop(iter)
472 :
473 208 : CALL timestop(handle)
474 :
475 416 : END SUBROUTINE allocate_blocks_v_L
476 :
477 : ! **************************************************************************************************
478 : !> \brief ...
479 : !> \param blocks_v_kp ...
480 : !> \param matrix_v_kp ...
481 : !> \param ikp_start ...
482 : !> \param ikp_end ...
483 : ! **************************************************************************************************
484 104 : SUBROUTINE allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
485 : TYPE(two_d_util_type), ALLOCATABLE, &
486 : DIMENSION(:, :, :) :: blocks_v_kp
487 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
488 : INTEGER :: ikp_start, ikp_end
489 :
490 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_blocks_v_kp'
491 :
492 : INTEGER :: col, handle, i_block, i_real_img, ik, &
493 : nblocks, row
494 104 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block
495 : TYPE(dbcsr_iterator_type) :: iter
496 :
497 104 : CALL timeset(routineN, handle)
498 :
499 104 : nblocks = 0
500 :
501 104 : CALL dbcsr_iterator_start(iter, matrix_v_kp(ikp_start, 1)%matrix)
502 :
503 467 : DO WHILE (dbcsr_iterator_blocks_left(iter))
504 :
505 363 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
506 :
507 363 : nblocks = nblocks + 1
508 :
509 : END DO
510 :
511 104 : CALL dbcsr_iterator_stop(iter)
512 :
513 6593 : ALLOCATE (blocks_v_kp(ikp_start:ikp_end, SIZE(matrix_v_kp, 2), nblocks))
514 :
515 750 : DO ik = ikp_start, ikp_end
516 :
517 2042 : DO i_real_img = 1, SIZE(matrix_v_kp, 2)
518 :
519 1292 : i_block = 1
520 :
521 1292 : CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_img)%matrix)
522 :
523 6276 : DO WHILE (dbcsr_iterator_blocks_left(iter))
524 :
525 4984 : CALL dbcsr_iterator_next_block(iter, row, col, data_block)
526 :
527 0 : ALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block(SIZE(data_block, 1), &
528 19936 : SIZE(data_block, 2)))
529 2452768 : blocks_v_kp(ik, i_real_img, i_block)%block = 0.0_dp
530 :
531 11260 : i_block = i_block + 1
532 :
533 : END DO
534 :
535 3230 : CALL dbcsr_iterator_stop(iter)
536 :
537 : END DO
538 :
539 : END DO
540 :
541 104 : CALL timestop(handle)
542 :
543 104 : END SUBROUTINE allocate_blocks_v_kp
544 :
545 : ! **************************************************************************************************
546 : !> \brief ...
547 : !> \param cell ...
548 : !> \param kpoints ...
549 : !> \param size_lattice_sum ...
550 : !> \param factor ...
551 : !> \param hmat ...
552 : !> \param x_min ...
553 : !> \param x_max ...
554 : !> \param y_min ...
555 : !> \param y_max ...
556 : !> \param z_min ...
557 : !> \param z_max ...
558 : !> \param nkp_grid ...
559 : ! **************************************************************************************************
560 120 : SUBROUTINE get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
561 : x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
562 :
563 : TYPE(cell_type), POINTER :: cell
564 : TYPE(kpoint_type), POINTER :: kpoints
565 : INTEGER :: size_lattice_sum, factor
566 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
567 : INTEGER :: x_min, x_max, y_min, y_max, z_min, z_max
568 : INTEGER, DIMENSION(3) :: nkp_grid
569 :
570 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_factor_and_xyz_min_max'
571 :
572 : INTEGER :: handle, nkp
573 : INTEGER, DIMENSION(3) :: periodic
574 :
575 120 : CALL timeset(routineN, handle)
576 :
577 120 : CALL get_kpoint_info(kpoints, nkp_grid=nkp_grid, nkp=nkp)
578 120 : CALL get_cell(cell=cell, h=hmat, periodic=periodic)
579 :
580 120 : IF (periodic(1) == 0) THEN
581 94 : CPASSERT(nkp_grid(1) == 1)
582 : END IF
583 120 : IF (periodic(2) == 0) THEN
584 24 : CPASSERT(nkp_grid(2) == 1)
585 : END IF
586 120 : IF (periodic(3) == 0) THEN
587 28 : CPASSERT(nkp_grid(3) == 1)
588 : END IF
589 :
590 120 : IF (MODULO(nkp_grid(1), 2) == 1) THEN
591 94 : factor = 3**(size_lattice_sum - 1)
592 26 : ELSE IF (MODULO(nkp_grid(1), 2) == 0) THEN
593 26 : factor = 2**(size_lattice_sum - 1)
594 : END IF
595 :
596 120 : IF (MODULO(nkp_grid(1), 2) == 1) THEN
597 94 : x_min = -(factor*nkp_grid(1) - 1)/2
598 94 : x_max = (factor*nkp_grid(1) - 1)/2
599 26 : ELSE IF (MODULO(nkp_grid(1), 2) == 0) THEN
600 26 : x_min = -factor*nkp_grid(1)/2
601 26 : x_max = factor*nkp_grid(1)/2 - 1
602 : END IF
603 120 : IF (periodic(1) == 0) THEN
604 94 : x_min = 0
605 94 : x_max = 0
606 : END IF
607 :
608 120 : IF (MODULO(nkp_grid(2), 2) == 1) THEN
609 24 : y_min = -(factor*nkp_grid(2) - 1)/2
610 24 : y_max = (factor*nkp_grid(2) - 1)/2
611 96 : ELSE IF (MODULO(nkp_grid(2), 2) == 0) THEN
612 96 : y_min = -factor*nkp_grid(2)/2
613 96 : y_max = factor*nkp_grid(2)/2 - 1
614 : END IF
615 120 : IF (periodic(2) == 0) THEN
616 24 : y_min = 0
617 24 : y_max = 0
618 : END IF
619 :
620 120 : IF (MODULO(nkp_grid(3), 2) == 1) THEN
621 28 : z_min = -(factor*nkp_grid(3) - 1)/2
622 28 : z_max = (factor*nkp_grid(3) - 1)/2
623 92 : ELSE IF (MODULO(nkp_grid(3), 2) == 0) THEN
624 92 : z_min = -factor*nkp_grid(3)/2
625 92 : z_max = factor*nkp_grid(3)/2 - 1
626 : END IF
627 120 : IF (periodic(3) == 0) THEN
628 28 : z_min = 0
629 28 : z_max = 0
630 : END IF
631 :
632 120 : CALL timestop(handle)
633 :
634 120 : END SUBROUTINE get_factor_and_xyz_min_max
635 :
636 : ! **************************************************************************************************
637 : !> \brief ...
638 : !> \param matrix_v_L_tmp ...
639 : !> \param matrix_v_kp ...
640 : !> \param ikp_start ...
641 : ! **************************************************************************************************
642 104 : SUBROUTINE allocate_tmp(matrix_v_L_tmp, matrix_v_kp, ikp_start)
643 :
644 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
645 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp
646 : INTEGER :: ikp_start
647 :
648 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_tmp'
649 :
650 : INTEGER :: handle
651 :
652 104 : CALL timeset(routineN, handle)
653 :
654 104 : NULLIFY (matrix_v_L_tmp)
655 104 : CALL dbcsr_init_p(matrix_v_L_tmp)
656 : CALL dbcsr_create(matrix=matrix_v_L_tmp, &
657 : template=matrix_v_kp(ikp_start, 1)%matrix, &
658 104 : matrix_type=dbcsr_type_no_symmetry)
659 104 : CALL dbcsr_reserve_all_blocks(matrix_v_L_tmp)
660 104 : CALL dbcsr_set(matrix_v_L_tmp, 0.0_dp)
661 :
662 104 : CALL timestop(handle)
663 :
664 104 : END SUBROUTINE allocate_tmp
665 :
666 : ! **************************************************************************************************
667 : !> \brief ...
668 : !> \param matrix_v_L_tmp ...
669 : ! **************************************************************************************************
670 104 : SUBROUTINE deallocate_tmp(matrix_v_L_tmp)
671 :
672 : TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp
673 :
674 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_tmp'
675 :
676 : INTEGER :: handle
677 :
678 104 : CALL timeset(routineN, handle)
679 :
680 104 : CALL dbcsr_release_p(matrix_v_L_tmp)
681 :
682 104 : CALL timestop(handle)
683 :
684 104 : END SUBROUTINE deallocate_tmp
685 :
686 : ! **************************************************************************************************
687 : !> \brief ...
688 : !> \param V_k ...
689 : !> \param qs_env ...
690 : !> \param kpoints ...
691 : !> \param size_lattice_sum ...
692 : !> \param basis_type ...
693 : !> \param ikp_start ...
694 : !> \param ikp_end ...
695 : ! **************************************************************************************************
696 16 : SUBROUTINE build_2c_coulomb_matrix_kp_small_cell(V_k, qs_env, kpoints, size_lattice_sum, &
697 : basis_type, ikp_start, ikp_end)
698 : COMPLEX(KIND=dp), DIMENSION(:, :, :) :: V_k
699 : TYPE(qs_environment_type), POINTER :: qs_env
700 : TYPE(kpoint_type), POINTER :: kpoints
701 : INTEGER :: size_lattice_sum
702 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
703 : INTEGER :: ikp_start, ikp_end
704 :
705 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_2c_coulomb_matrix_kp_small_cell'
706 :
707 : INTEGER :: factor, handle, handle2, i_cell, i_x, i_x_inner, i_x_outer, ik, ikp_local, j_y, &
708 : j_y_inner, j_y_outer, k_z, k_z_inner, k_z_outer, n_atom, n_bf, x_max, x_min, y_max, &
709 : y_min, z_max, z_min
710 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bf_end_from_atom, bf_start_from_atom
711 : INTEGER, DIMENSION(3) :: nkp_grid
712 : REAL(KIND=dp) :: coskl, sinkl
713 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: V_L
714 : REAL(KIND=dp), DIMENSION(3) :: vec_L, vec_s
715 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
716 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
717 : TYPE(cell_type), POINTER :: cell
718 : TYPE(mp_para_env_type), POINTER :: para_env
719 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
720 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
721 :
722 16 : CALL timeset(routineN, handle)
723 :
724 : CALL get_qs_env(qs_env=qs_env, &
725 : para_env=para_env, &
726 : particle_set=particle_set, &
727 : cell=cell, &
728 : qs_kind_set=qs_kind_set, &
729 16 : atomic_kind_set=atomic_kind_set)
730 :
731 : CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
732 16 : x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
733 :
734 16 : CALL get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
735 :
736 64 : ALLOCATE (V_L(n_bf, n_bf))
737 :
738 232 : DO i_x_inner = 0, 2*nkp_grid(1) - 1
739 11624 : DO j_y_inner = 0, 2*nkp_grid(2) - 1
740 93528 : DO k_z_inner = 0, 2*nkp_grid(3) - 1
741 :
742 3276800 : V_L(:, :) = 0.0_dp
743 81920 : i_cell = 0
744 :
745 204800 : DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
746 696320 : DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
747 2334720 : DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
748 :
749 1720320 : i_x = i_x_inner + i_x_outer
750 1720320 : j_y = j_y_inner + j_y_outer
751 1720320 : k_z = k_z_inner + k_z_outer
752 :
753 : IF (i_x > x_max .OR. i_x < x_min .OR. &
754 : j_y > y_max .OR. j_y < y_min .OR. &
755 1720320 : k_z > z_max .OR. k_z < z_min) CYCLE
756 :
757 663040 : i_cell = i_cell + 1
758 :
759 2652160 : vec_s = [REAL(i_x, dp), REAL(j_y, dp), REAL(k_z, dp)]
760 :
761 663040 : IF (MODULO(i_cell, para_env%num_pe) /= para_env%mepos) CYCLE
762 :
763 4309760 : vec_L = MATMUL(hmat, vec_s)
764 :
765 : ! Compute (P 0 | Q vec_L) and add it to V_R
766 : CALL add_V_L(V_L, vec_L, n_atom, bf_start_from_atom, bf_end_from_atom, &
767 2211840 : particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
768 :
769 : END DO
770 : END DO
771 : END DO
772 :
773 81920 : CALL para_env%sync()
774 81920 : CALL para_env%sum(V_L)
775 :
776 81920 : CALL timeset(routineN//"_R_to_k", handle2)
777 :
778 81920 : ikp_local = 0
779 :
780 : ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
781 44122112 : DO ik = 1, ikp_end
782 :
783 44040192 : IF (MODULO(ik, para_env%num_pe) /= para_env%mepos) CYCLE
784 :
785 22020096 : ikp_local = ikp_local + 1
786 :
787 22020096 : IF (ik < ikp_start) CYCLE
788 :
789 : ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
790 71303168 : coskl = COS(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
791 71303168 : sinkl = SIN(twopi*DOT_PRODUCT(vec_s(1:3), kpoints%xkp(1:3, ik)))
792 :
793 : V_k(:, :, ikp_local) = V_k(:, :, ikp_local) + z_one*coskl*V_L(:, :) + &
794 717307904 : gaussi*sinkl*V_L(:, :)
795 :
796 : END DO
797 :
798 175232 : CALL timestop(handle2)
799 :
800 : END DO
801 : END DO
802 : END DO
803 :
804 16 : CALL timestop(handle)
805 :
806 32 : END SUBROUTINE build_2c_coulomb_matrix_kp_small_cell
807 :
808 : ! **************************************************************************************************
809 : !> \brief ...
810 : !> \param qs_env ...
811 : !> \param n_atom ...
812 : !> \param basis_type ...
813 : !> \param bf_start_from_atom ...
814 : !> \param bf_end_from_atom ...
815 : !> \param n_bf ...
816 : ! **************************************************************************************************
817 16 : SUBROUTINE get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
818 :
819 : TYPE(qs_environment_type), POINTER :: qs_env
820 : INTEGER :: n_atom
821 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
822 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bf_start_from_atom, bf_end_from_atom
823 : INTEGER :: n_bf
824 :
825 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_basis_sizes'
826 :
827 : INTEGER :: handle, iatom, ikind, n_kind, nsgf
828 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
829 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
830 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
831 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
832 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
833 :
834 16 : CALL timeset(routineN, handle)
835 :
836 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
837 16 : qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
838 16 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
839 :
840 16 : n_atom = SIZE(particle_set)
841 16 : n_kind = SIZE(qs_kind_set)
842 :
843 48 : DO ikind = 1, n_kind
844 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
845 32 : basis_type=basis_type)
846 48 : CPASSERT(ASSOCIATED(basis_set_a))
847 : END DO
848 :
849 80 : ALLOCATE (bf_start_from_atom(n_atom), bf_end_from_atom(n_atom))
850 :
851 16 : n_bf = 0
852 60 : DO iatom = 1, n_atom
853 44 : bf_start_from_atom(iatom) = n_bf + 1
854 44 : ikind = kind_of(iatom)
855 44 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type=basis_type)
856 44 : n_bf = n_bf + nsgf
857 60 : bf_end_from_atom(iatom) = n_bf
858 : END DO
859 :
860 16 : CALL timestop(handle)
861 :
862 32 : END SUBROUTINE get_basis_sizes
863 :
864 : ! **************************************************************************************************
865 : !> \brief ...
866 : !> \param V_L ...
867 : !> \param vec_L ...
868 : !> \param n_atom ...
869 : !> \param bf_start_from_atom ...
870 : !> \param bf_end_from_atom ...
871 : !> \param particle_set ...
872 : !> \param qs_kind_set ...
873 : !> \param atomic_kind_set ...
874 : !> \param basis_type ...
875 : !> \param cell ...
876 : ! **************************************************************************************************
877 331520 : SUBROUTINE add_V_L(V_L, vec_L, n_atom, bf_start_from_atom, bf_end_from_atom, &
878 : particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
879 :
880 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: V_L
881 : REAL(KIND=dp), DIMENSION(3) :: vec_L
882 : INTEGER :: n_atom
883 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bf_start_from_atom, bf_end_from_atom
884 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
885 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
886 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
887 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
888 : TYPE(cell_type), POINTER :: cell
889 :
890 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_V_L'
891 :
892 : INTEGER :: a_1, a_2, atom_a, atom_b, b_1, b_2, &
893 : handle, kind_a, kind_b
894 331520 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
895 : REAL(dp), DIMENSION(3) :: ra, rab_L, rb
896 331520 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: V_L_ab
897 331520 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: contr_a, contr_b
898 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
899 :
900 331520 : CALL timeset(routineN, handle)
901 :
902 331520 : NULLIFY (basis_set_a, basis_set_b)
903 :
904 331520 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
905 :
906 1305600 : DO atom_a = 1, n_atom
907 :
908 4186880 : DO atom_b = 1, n_atom
909 :
910 2881280 : kind_a = kind_of(atom_a)
911 2881280 : kind_b = kind_of(atom_b)
912 :
913 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, &
914 2881280 : basis_type=basis_type)
915 : CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, &
916 2881280 : basis_type=basis_type)
917 :
918 2881280 : ra(1:3) = pbc(particle_set(atom_a)%r(1:3), cell)
919 2881280 : rb(1:3) = pbc(particle_set(atom_b)%r(1:3), cell)
920 :
921 11525120 : rab_L(1:3) = rb(1:3) - ra(1:3) + vec_L(1:3)
922 :
923 2881280 : CALL contraction_matrix_shg(basis_set_a, contr_a)
924 2881280 : CALL contraction_matrix_shg(basis_set_b, contr_b)
925 :
926 2881280 : a_1 = bf_start_from_atom(atom_a)
927 2881280 : a_2 = bf_end_from_atom(atom_a)
928 2881280 : b_1 = bf_start_from_atom(atom_b)
929 2881280 : b_2 = bf_end_from_atom(atom_b)
930 :
931 11525120 : ALLOCATE (V_L_ab(a_2 - a_1 + 1, b_2 - b_1 + 1))
932 :
933 : CALL int_operators_r12_ab_shg(operator_coulomb, V_L_ab, rab=rab_L, &
934 : fba=basis_set_a, fbb=basis_set_b, &
935 : scona_shg=contr_a, sconb_shg=contr_b, &
936 2881280 : calculate_forces=.FALSE.)
937 :
938 20394240 : V_L(a_1:a_2, b_1:b_2) = V_L(a_1:a_2, b_1:b_2) + V_L_ab(:, :)
939 :
940 3855360 : DEALLOCATE (contr_a, contr_b, V_L_ab)
941 :
942 : END DO
943 :
944 : END DO
945 :
946 331520 : DEALLOCATE (kind_of)
947 :
948 331520 : CALL timestop(handle)
949 :
950 331520 : END SUBROUTINE add_V_L
951 :
952 0 : END MODULE kpoint_coulomb_2c
|
