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