(git:936074a)
Loading...
Searching...
No Matches
kpoint_coulomb_2c.F
Go to the documentation of this file.
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! **************************************************************************************************
19 USE cell_types, ONLY: cell_type,&
20 get_cell,&
21 pbc
23 USE cp_dbcsr_api, ONLY: &
26 dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
30 USE kinds, ONLY: dp
31 USE kpoint_types, ONLY: get_kpoint_info,&
33 USE mathconstants, ONLY: gaussi,&
34 twopi,&
35 z_one
40 USE qs_kind_types, ONLY: get_qs_kind,&
42#include "./base/base_uses.f90"
43
44 IMPLICIT NONE
45
46 PRIVATE
47
48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_coulomb_2c'
49
51
52! **************************************************************************************************
53
54 TYPE two_d_util_type
55 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: block
56 END TYPE two_d_util_type
57
58CONTAINS
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 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 CALL timeset(routinen, handle)
93
94 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 operator_type, ikp_start, ikp_end)
99
100 CALL deallocate_tmp(matrix_v_l_tmp)
101
102 CALL timestop(handle)
103
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 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 TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_v_l, blocks_v_l_store
145 TYPE(two_d_util_type), ALLOCATABLE, &
146 DIMENSION(:, :, :) :: blocks_v_kp
147
148 CALL timeset(routinen, handle)
149
150 CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
151 x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
152
153 CALL allocate_blocks_v_kp(blocks_v_kp, matrix_v_kp, ikp_start, ikp_end)
154 CALL allocate_blocks_v_l(blocks_v_l, matrix_v_l_tmp)
155 CALL allocate_blocks_v_l(blocks_v_l_store, matrix_v_l_tmp)
156
157 DO i_x_inner = 0, 2*nkp_grid(1) - 1
158 DO j_y_inner = 0, 2*nkp_grid(2) - 1
159 DO k_z_inner = 0, 2*nkp_grid(3) - 1
160
161 DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
162 DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
163 DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
164
165 i_x = i_x_inner + i_x_outer
166 j_y = j_y_inner + j_y_outer
167 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 k_z > z_max .OR. k_z < z_min) cycle
172
173 vec_s = [real(i_x, dp), real(j_y, dp), real(k_z, dp)]
174
175 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 operator_type)
181
182!$OMP PARALLEL DO DEFAULT(NONE) &
183!$OMP SHARED(blocks_v_L, blocks_v_L_store) &
184!$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 CALL timeset(routinen//"_R_to_k", handle2)
196
197 ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
198 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 coskl = cos(twopi*dot_product(vec_s(1:3), kpoints%xkp(1:3, ik)))
202 sinkl = sin(twopi*dot_product(vec_s(1:3), kpoints%xkp(1:3, ik)))
203
204 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 + 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 + sinkl*blocks_v_l_store(i_block)%block(:, :)
210
211 END DO
212
213 END DO
214
215 DO i_block = 1, SIZE(blocks_v_l)
216
217 blocks_v_l_store(i_block)%block(:, :) = 0.0_dp
218
219 END DO
220
221 CALL timestop(handle2)
222
223 END DO
224 END DO
225 END DO
226
227 CALL set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_v_kp, ikp_start, ikp_end)
228
229 CALL deallocate_blocks_v_kp(blocks_v_kp)
230 CALL deallocate_blocks_v_l(blocks_v_l)
231 CALL deallocate_blocks_v_l(blocks_v_l_store)
232
233 CALL timestop(handle)
234
235 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
255 TYPE(dbcsr_iterator_type) :: iter
256
257 CALL timeset(routinen, handle)
258
259 DO ik = ikp_start, ikp_end
260
261 DO i_real_im = 1, 2
262
263 i_block = 1
264
265 CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_im)%matrix)
266
267 DO WHILE (dbcsr_iterator_blocks_left(iter))
268
269 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
270
271 data_block(:, :) = blocks_v_kp(ik, i_real_im, i_block)%block(:, :)
272
273 i_block = i_block + 1
274
275 END DO
276
277 CALL dbcsr_iterator_stop(iter)
278
279 END DO
280
281 END DO
282
283 CALL timestop(handle)
284
285 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
316 REAL(dp), DIMENSION(3) :: ra, rab_l, rb
317 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
318 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 CALL timeset(routinen, handle)
323
324 NULLIFY (basis_set_a, basis_set_b, data_block)
325
326 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
327
328 CALL dbcsr_set(matrix_v_l_tmp, 0.0_dp)
329
330 i_block = 1
331
332 CALL dbcsr_iterator_start(iter, matrix_v_l_tmp)
333
334 DO WHILE (dbcsr_iterator_blocks_left(iter))
335
336 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
337
338 kind_a = kind_of(row)
339 kind_b = kind_of(col)
340
341 CALL get_qs_kind(qs_kind=qs_kind_set(kind_a), basis_set=basis_set_a, basis_type=basis_type)
342 CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, basis_type=basis_type)
343
344 ra(1:3) = pbc(particle_set(row)%r(1:3), cell)
345 rb(1:3) = pbc(particle_set(col)%r(1:3), cell)
346
347 rab_l(1:3) = rb(1:3) - ra(1:3) + vec_l(1:3)
348
349 CALL contraction_matrix_shg(basis_set_a, contr_a)
350 CALL contraction_matrix_shg(basis_set_b, contr_b)
351
352 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 calculate_forces=.false.)
357
358 i_block = i_block + 1
359
360 DEALLOCATE (contr_a, contr_b)
361
362 END DO
363
364 CALL dbcsr_iterator_stop(iter)
365
366 DEALLOCATE (kind_of)
367
368 CALL timestop(handle)
369
370 END SUBROUTINE compute_v_transl
371
372! **************************************************************************************************
373!> \brief ...
374!> \param blocks_v_kp ...
375! **************************************************************************************************
376 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 CALL timeset(routinen, handle)
385
386 DO ik = lbound(blocks_v_kp, 1), ubound(blocks_v_kp, 1)
387 DO i_real_img = 1, SIZE(blocks_v_kp, 2)
388 DO i_block = 1, SIZE(blocks_v_kp, 3)
389 DEALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block)
390 END DO
391 END DO
392 END DO
393
394 DEALLOCATE (blocks_v_kp)
395
396 CALL timestop(handle)
397
398 END SUBROUTINE deallocate_blocks_v_kp
399
400! **************************************************************************************************
401!> \brief ...
402!> \param blocks_v_L ...
403! **************************************************************************************************
404 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 CALL timeset(routinen, handle)
412
413 DO i_block = 1, SIZE(blocks_v_l, 1)
414 DEALLOCATE (blocks_v_l(i_block)%block)
415 END DO
416
417 DEALLOCATE (blocks_v_l)
418
419 CALL timestop(handle)
420
421 END SUBROUTINE deallocate_blocks_v_l
422
423! **************************************************************************************************
424!> \brief ...
425!> \param blocks_v_L ...
426!> \param matrix_v_L_tmp ...
427! **************************************************************************************************
428 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
436 TYPE(dbcsr_iterator_type) :: iter
437
438 CALL timeset(routinen, handle)
439
440 nblocks = 0
441
442 CALL dbcsr_iterator_start(iter, matrix_v_l_tmp)
443
444 DO WHILE (dbcsr_iterator_blocks_left(iter))
445
446 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
447
448 nblocks = nblocks + 1
449
450 END DO
451
452 CALL dbcsr_iterator_stop(iter)
453
454 ALLOCATE (blocks_v_l(nblocks))
455
456 i_block = 1
457
458 CALL dbcsr_iterator_start(iter, matrix_v_l_tmp)
459
460 DO WHILE (dbcsr_iterator_blocks_left(iter))
461
462 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
463
464 ALLOCATE (blocks_v_l(i_block)%block(SIZE(data_block, 1), SIZE(data_block, 2)))
465 blocks_v_l(i_block)%block = 0.0_dp
466
467 i_block = i_block + 1
468
469 END DO
470
471 CALL dbcsr_iterator_stop(iter)
472
473 CALL timestop(handle)
474
475 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
495 TYPE(dbcsr_iterator_type) :: iter
496
497 CALL timeset(routinen, handle)
498
499 nblocks = 0
500
501 CALL dbcsr_iterator_start(iter, matrix_v_kp(ikp_start, 1)%matrix)
502
503 DO WHILE (dbcsr_iterator_blocks_left(iter))
504
505 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
506
507 nblocks = nblocks + 1
508
509 END DO
510
511 CALL dbcsr_iterator_stop(iter)
512
513 ALLOCATE (blocks_v_kp(ikp_start:ikp_end, SIZE(matrix_v_kp, 2), nblocks))
514
515 DO ik = ikp_start, ikp_end
516
517 DO i_real_img = 1, SIZE(matrix_v_kp, 2)
518
519 i_block = 1
520
521 CALL dbcsr_iterator_start(iter, matrix_v_kp(ik, i_real_img)%matrix)
522
523 DO WHILE (dbcsr_iterator_blocks_left(iter))
524
525 CALL dbcsr_iterator_next_block(iter, row, col, data_block)
526
527 ALLOCATE (blocks_v_kp(ik, i_real_img, i_block)%block(SIZE(data_block, 1), &
528 SIZE(data_block, 2)))
529 blocks_v_kp(ik, i_real_img, i_block)%block = 0.0_dp
530
531 i_block = i_block + 1
532
533 END DO
534
535 CALL dbcsr_iterator_stop(iter)
536
537 END DO
538
539 END DO
540
541 CALL timestop(handle)
542
543 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 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 CALL timeset(routinen, handle)
576
577 CALL get_kpoint_info(kpoints, nkp_grid=nkp_grid, nkp=nkp)
578 CALL get_cell(cell=cell, h=hmat, periodic=periodic)
579
580 IF (periodic(1) == 0) THEN
581 cpassert(nkp_grid(1) == 1)
582 END IF
583 IF (periodic(2) == 0) THEN
584 cpassert(nkp_grid(2) == 1)
585 END IF
586 IF (periodic(3) == 0) THEN
587 cpassert(nkp_grid(3) == 1)
588 END IF
589
590 IF (modulo(nkp_grid(1), 2) == 1) THEN
591 factor = 3**(size_lattice_sum - 1)
592 ELSE IF (modulo(nkp_grid(1), 2) == 0) THEN
593 factor = 2**(size_lattice_sum - 1)
594 END IF
595
596 IF (modulo(nkp_grid(1), 2) == 1) THEN
597 x_min = -(factor*nkp_grid(1) - 1)/2
598 x_max = (factor*nkp_grid(1) - 1)/2
599 ELSE IF (modulo(nkp_grid(1), 2) == 0) THEN
600 x_min = -factor*nkp_grid(1)/2
601 x_max = factor*nkp_grid(1)/2 - 1
602 END IF
603 IF (periodic(1) == 0) THEN
604 x_min = 0
605 x_max = 0
606 END IF
607
608 IF (modulo(nkp_grid(2), 2) == 1) THEN
609 y_min = -(factor*nkp_grid(2) - 1)/2
610 y_max = (factor*nkp_grid(2) - 1)/2
611 ELSE IF (modulo(nkp_grid(2), 2) == 0) THEN
612 y_min = -factor*nkp_grid(2)/2
613 y_max = factor*nkp_grid(2)/2 - 1
614 END IF
615 IF (periodic(2) == 0) THEN
616 y_min = 0
617 y_max = 0
618 END IF
619
620 IF (modulo(nkp_grid(3), 2) == 1) THEN
621 z_min = -(factor*nkp_grid(3) - 1)/2
622 z_max = (factor*nkp_grid(3) - 1)/2
623 ELSE IF (modulo(nkp_grid(3), 2) == 0) THEN
624 z_min = -factor*nkp_grid(3)/2
625 z_max = factor*nkp_grid(3)/2 - 1
626 END IF
627 IF (periodic(3) == 0) THEN
628 z_min = 0
629 z_max = 0
630 END IF
631
632 CALL timestop(handle)
633
634 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 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 CALL timeset(routinen, handle)
653
654 NULLIFY (matrix_v_l_tmp)
655 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 matrix_type=dbcsr_type_no_symmetry)
659 CALL dbcsr_reserve_all_blocks(matrix_v_l_tmp)
660 CALL dbcsr_set(matrix_v_l_tmp, 0.0_dp)
661
662 CALL timestop(handle)
663
664 END SUBROUTINE allocate_tmp
665
666! **************************************************************************************************
667!> \brief ...
668!> \param matrix_v_L_tmp ...
669! **************************************************************************************************
670 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 CALL timeset(routinen, handle)
679
680 CALL dbcsr_release_p(matrix_v_l_tmp)
681
682 CALL timestop(handle)
683
684 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 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 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 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
720 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
721
722 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 atomic_kind_set=atomic_kind_set)
730
731 CALL get_factor_and_xyz_min_max(cell, kpoints, size_lattice_sum, factor, hmat, &
732 x_min, x_max, y_min, y_max, z_min, z_max, nkp_grid)
733
734 CALL get_basis_sizes(qs_env, n_atom, basis_type, bf_start_from_atom, bf_end_from_atom, n_bf)
735
736 ALLOCATE (v_l(n_bf, n_bf))
737
738 DO i_x_inner = 0, 2*nkp_grid(1) - 1
739 DO j_y_inner = 0, 2*nkp_grid(2) - 1
740 DO k_z_inner = 0, 2*nkp_grid(3) - 1
741
742 v_l(:, :) = 0.0_dp
743 i_cell = 0
744
745 DO i_x_outer = x_min, x_max + nkp_grid(1), 2*nkp_grid(1)
746 DO j_y_outer = y_min, y_max + nkp_grid(2), 2*nkp_grid(2)
747 DO k_z_outer = z_min, z_max + nkp_grid(3), 2*nkp_grid(3)
748
749 i_x = i_x_inner + i_x_outer
750 j_y = j_y_inner + j_y_outer
751 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 k_z > z_max .OR. k_z < z_min) cycle
756
757 i_cell = i_cell + 1
758
759 vec_s = [real(i_x, dp), real(j_y, dp), real(k_z, dp)]
760
761 IF (modulo(i_cell, para_env%num_pe) /= para_env%mepos) cycle
762
763 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 particle_set, qs_kind_set, atomic_kind_set, basis_type, cell)
768
769 END DO
770 END DO
771 END DO
772
773 CALL para_env%sync()
774 CALL para_env%sum(v_l)
775
776 CALL timeset(routinen//"_R_to_k", handle2)
777
778 ikp_local = 0
779
780 ! add exp(iq*vec_L) * (P 0 | Q vec_L) to V_PQ(q)
781 DO ik = 1, ikp_end
782
783 IF (modulo(ik, para_env%num_pe) /= para_env%mepos) cycle
784
785 ikp_local = ikp_local + 1
786
787 IF (ik < ikp_start) cycle
788
789 ! coskl and sinkl are identical for all i_x_outer, j_y_outer, k_z_outer
790 coskl = cos(twopi*dot_product(vec_s(1:3), kpoints%xkp(1:3, ik)))
791 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 gaussi*sinkl*v_l(:, :)
795
796 END DO
797
798 CALL timestop(handle2)
799
800 END DO
801 END DO
802 END DO
803
804 CALL timestop(handle)
805
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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
829 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
830 TYPE(gto_basis_set_type), POINTER :: basis_set_a
831 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
832 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
833
834 CALL timeset(routinen, handle)
835
836 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
837 qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
838 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
839
840 n_atom = SIZE(particle_set)
841 n_kind = SIZE(qs_kind_set)
842
843 DO ikind = 1, n_kind
844 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
845 basis_type=basis_type)
846 cpassert(ASSOCIATED(basis_set_a))
847 END DO
848
849 ALLOCATE (bf_start_from_atom(n_atom), bf_end_from_atom(n_atom))
850
851 n_bf = 0
852 DO iatom = 1, n_atom
853 bf_start_from_atom(iatom) = n_bf + 1
854 ikind = kind_of(iatom)
855 CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type=basis_type)
856 n_bf = n_bf + nsgf
857 bf_end_from_atom(iatom) = n_bf
858 END DO
859
860 CALL timestop(handle)
861
862 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
895 REAL(dp), DIMENSION(3) :: ra, rab_l, rb
896 REAL(kind=dp), DIMENSION(:, :), POINTER :: v_l_ab
897 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: contr_a, contr_b
898 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
899
900 CALL timeset(routinen, handle)
901
902 NULLIFY (basis_set_a, basis_set_b)
903
904 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
905
906 DO atom_a = 1, n_atom
907
908 DO atom_b = 1, n_atom
909
910 kind_a = kind_of(atom_a)
911 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 basis_type=basis_type)
915 CALL get_qs_kind(qs_kind=qs_kind_set(kind_b), basis_set=basis_set_b, &
916 basis_type=basis_type)
917
918 ra(1:3) = pbc(particle_set(atom_a)%r(1:3), cell)
919 rb(1:3) = pbc(particle_set(atom_b)%r(1:3), cell)
920
921 rab_l(1:3) = rb(1:3) - ra(1:3) + vec_l(1:3)
922
923 CALL contraction_matrix_shg(basis_set_a, contr_a)
924 CALL contraction_matrix_shg(basis_set_b, contr_b)
925
926 a_1 = bf_start_from_atom(atom_a)
927 a_2 = bf_end_from_atom(atom_a)
928 b_1 = bf_start_from_atom(atom_b)
929 b_2 = bf_end_from_atom(atom_b)
930
931 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 calculate_forces=.false.)
937
938 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 DEALLOCATE (contr_a, contr_b, v_l_ab)
941
942 END DO
943
944 END DO
945
946 DEALLOCATE (kind_of)
947
948 CALL timestop(handle)
949
950 END SUBROUTINE add_v_l
951
952END MODULE kpoint_coulomb_2c
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:200
constants for the different operators of the 2c-integrals
integer, parameter, public operator_coulomb
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation of contracte...
subroutine, public contraction_matrix_shg(basis, scon_shg)
contraction matrix for SHG integrals
Calculation of contracted, spherical Gaussian integrals using the solid harmonic Gaussian (SHG) integ...
subroutine, public int_operators_r12_ab_shg(r12_operator, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, omega, calculate_forces)
Calcululates the two-center integrals of the type (a|O(r12)|b) using the SHG scheme.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice summation in...
subroutine, public build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
...
subroutine, public build_2c_coulomb_matrix_kp_small_cell(v_k, qs_env, kpoints, size_lattice_sum, basis_type, ikp_start, ikp_end)
...
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition of mathematical constants and functions.
complex(kind=dp), parameter, public z_one
complex(kind=dp), parameter, public gaussi
real(kind=dp), parameter, public twopi
Interface to the message passing library MPI.
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.