(git:0e87843)
Loading...
Searching...
No Matches
rpa_gw_im_time_util.F
Go to the documentation of this file.
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 Utility routines for GW with imaginary time
10!> \par History
11!> 06.2019 split from rpa_im_time.F [Frederick Stein]
12! **************************************************************************************************
14
17 USE cell_types, ONLY: cell_type,&
18 pbc
19 USE cp_dbcsr_api, ONLY: &
22 dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_type, dbcsr_type_no_symmetry
30 USE cp_fm_types, ONLY: cp_fm_create,&
35 USE dbt_api, ONLY: &
36 dbt_contract, dbt_copy, dbt_copy_matrix_to_tensor, dbt_create, dbt_default_distvec, &
37 dbt_destroy, dbt_get_info, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
38 USE hfx_types, ONLY: alloc_containers,&
41 USE kinds, ONLY: dp,&
42 int_8
43 USE mathconstants, ONLY: twopi
52 USE qs_tensors, ONLY: compress_tensor,&
59#include "./base/base_uses.f90"
60
61 IMPLICIT NONE
62
63 PRIVATE
64
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_im_time_util'
66
68
69CONTAINS
70
71! **************************************************************************************************
72!> \brief ...
73!> \param t_3c_overl_int ...
74!> \param t_3c_O_compressed ...
75!> \param t_3c_O_ind ...
76!> \param t_3c_overl_int_ao_mo ...
77!> \param t_3c_O_mo_compressed ...
78!> \param t_3c_O_mo_ind ...
79!> \param t_3c_overl_int_gw_RI ...
80!> \param t_3c_overl_int_gw_AO ...
81!> \param starts_array_mc ...
82!> \param ends_array_mc ...
83!> \param mo_coeff ...
84!> \param matrix_s ...
85!> \param gw_corr_lev_occ ...
86!> \param gw_corr_lev_virt ...
87!> \param homo ...
88!> \param nmo ...
89!> \param para_env ...
90!> \param do_ic_model ...
91!> \param t_3c_overl_nnP_ic ...
92!> \param t_3c_overl_nnP_ic_reflected ...
93!> \param qs_env ...
94!> \param unit_nr ...
95!> \param do_alpha ...
96! **************************************************************************************************
97 SUBROUTINE get_tensor_3c_overl_int_gw(t_3c_overl_int, t_3c_O_compressed, t_3c_O_ind, &
98 t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
99 t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
100 starts_array_mc, ends_array_mc, &
101 mo_coeff, matrix_s, &
102 gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
103 para_env, &
104 do_ic_model, &
105 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
106 qs_env, unit_nr, do_alpha)
107
108 TYPE(dbt_type), DIMENSION(:, :) :: t_3c_overl_int
109 TYPE(hfx_compression_type), DIMENSION(:, :, :) :: t_3c_o_compressed
110 TYPE(block_ind_type), DIMENSION(:, :, :) :: t_3c_o_ind
111 TYPE(dbt_type) :: t_3c_overl_int_ao_mo
112 TYPE(hfx_compression_type) :: t_3c_o_mo_compressed
113 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: t_3c_o_mo_ind
114 TYPE(dbt_type) :: t_3c_overl_int_gw_ri, &
115 t_3c_overl_int_gw_ao
116 INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc
117 TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
118 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
119 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo, &
120 nmo
121 TYPE(mp_para_env_type), INTENT(IN) :: para_env
122 LOGICAL, INTENT(IN) :: do_ic_model
123 TYPE(dbt_type) :: t_3c_overl_nnp_ic, &
124 t_3c_overl_nnp_ic_reflected
125 TYPE(qs_environment_type), POINTER :: qs_env
126 INTEGER, INTENT(IN) :: unit_nr
127 LOGICAL, INTENT(IN), OPTIONAL :: do_alpha
128
129 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_tensor_3c_overl_int_gw'
130
131 INTEGER :: cut_memory, handle, i_mem, icol_global, imo, irow_global, min_bsize, &
132 min_bsize_mo, nkind, nmo_blk_gw, npcols, nprows, size_mo, unit_nr_prv
133 INTEGER(int_8) :: nze
134 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist1, dist2, dist3, sizes_ao, &
135 sizes_ao_split, sizes_mo, sizes_mo_1, &
136 sizes_ri, sizes_ri_split, tmp
137 INTEGER, DIMENSION(2) :: pdims_2d
138 INTEGER, DIMENSION(2, 1) :: bounds
139 INTEGER, DIMENSION(2, 3) :: ibounds
140 INTEGER, DIMENSION(3) :: bounds_3c, pdims
141 INTEGER, DIMENSION(:), POINTER :: distp_1, distp_2, sizes_mo_blocked, &
142 sizes_mo_p1, sizes_mo_p2
143 LOGICAL :: memory_info, my_do_alpha
144 REAL(dp) :: compression_factor, memory_3c, occ
145 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: norm
146 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
147 TYPE(cp_fm_type) :: fm_mat_mo_coeff_gw
148 TYPE(dbcsr_distribution_type) :: dist, dist_templ
149 TYPE(dbcsr_type) :: mat_mo_coeff_gw_reflected_norm, &
150 mat_norm, mat_norm_diag, mat_work
151 TYPE(dbcsr_type), POINTER :: mat_mo_coeff_gw, &
152 mat_mo_coeff_gw_reflected
153 TYPE(dbt_pgrid_type) :: pgrid_2d, pgrid_ao, pgrid_ic, pgrid_mo
154 TYPE(dbt_type) :: mo_coeff_gw_t, mo_coeff_gw_t_tmp, &
155 t_3c_overl_int_ao_ao, &
156 t_3c_overl_int_mo_ao, &
157 t_3c_overl_int_mo_mo
158 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ao, basis_set_ri_aux
159 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
160 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
161
162 memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info
163 IF (memory_info) THEN
164 unit_nr_prv = unit_nr
165 ELSE
166 unit_nr_prv = 0
167 END IF
168
169 my_do_alpha = .false.
170 IF (PRESENT(do_alpha)) my_do_alpha = do_alpha
171
172 CALL timeset(routinen, handle)
173
174 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, atomic_kind_set=atomic_kind_set)
175
176 CALL cp_fm_create(fm_mat_mo_coeff_gw, mo_coeff%matrix_struct)
177 CALL cp_fm_to_fm(mo_coeff, fm_mat_mo_coeff_gw)
178
179 ! set MO coeffs to zero where
180 DO irow_global = 1, nmo
181 DO icol_global = 1, homo - gw_corr_lev_occ
182 CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
183 END DO
184 DO icol_global = homo + gw_corr_lev_virt + 1, nmo
185 CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
186 END DO
187 END DO
188
189 NULLIFY (mat_mo_coeff_gw)
190 CALL dbcsr_init_p(mat_mo_coeff_gw)
191
192 CALL cp_dbcsr_m_by_n_from_row_template(mat_mo_coeff_gw, template=matrix_s(1)%matrix, n=nmo, &
193 sym=dbcsr_type_no_symmetry)
194
195 CALL copy_fm_to_dbcsr(fm_mat_mo_coeff_gw, &
196 mat_mo_coeff_gw, &
197 keep_sparsity=.false.)
198
199 ! just remove the blocks which have been set to zero
200 CALL dbcsr_filter(mat_mo_coeff_gw, 1.0e-20_dp)
201
202 min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize
203 min_bsize_mo = qs_env%mp2_env%ri_rpa_im_time%min_bsize_mo
204
205 CALL split_block_sizes([gw_corr_lev_occ + gw_corr_lev_virt], sizes_mo, min_bsize_mo)
206 ALLOCATE (sizes_mo_1(nmo))
207 sizes_mo_1(:) = 1
208
209 nmo_blk_gw = SIZE(sizes_mo)
210 CALL move_alloc(sizes_mo, tmp)
211 ALLOCATE (sizes_mo(nmo_blk_gw + 2))
212 sizes_mo(1) = homo - gw_corr_lev_occ
213 sizes_mo(2:SIZE(tmp) + 1) = tmp(:)
214 sizes_mo(SIZE(tmp) + 2) = nmo - (homo + gw_corr_lev_virt)
215
216 ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
217 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
218 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
219 CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
220 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ao, basis=basis_set_ao)
221
222 CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_ao_split)
223 CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_ri_split)
224
225 DEALLOCATE (basis_set_ao, basis_set_ri_aux)
226
227 pdims = 0
228 CALL dbt_pgrid_create(para_env, pdims, pgrid_ao, &
229 tensor_dims=[SIZE(sizes_ri_split), SIZE(sizes_ao_split), SIZE(sizes_ao_split)])
230
231 pdims_2d = 0
232 CALL mp_dims_create(para_env%num_pe, pdims_2d)
233
234 ! we iterate over MO blocks for saving memory during contraction, thus we should not parallelize over MO dimension
235 pdims = [pdims_2d(1), pdims_2d(2), 1]
236 CALL dbt_pgrid_create(para_env, pdims, pgrid_mo, &
237 tensor_dims=[SIZE(sizes_ri_split), SIZE(sizes_ao_split), 1])
238
239 pdims_2d = 0
240 CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d, &
241 tensor_dims=[SIZE(sizes_ao_split), nmo])
242
243 CALL create_3c_tensor(t_3c_overl_int_ao_ao, dist1, dist2, dist3, pgrid_ao, &
244 sizes_ri_split, sizes_ao_split, sizes_ao_split, [1, 2], [3], name="(RI AO | AO)")
245 DEALLOCATE (dist1, dist2, dist3)
246
247 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
248 CALL create_3c_tensor(t_3c_overl_int_ao_mo, dist1, dist2, dist3, pgrid_ao, &
249 sizes_ri_split, sizes_ao_split, sizes_mo_1, [1, 2], [3], name="(RI AO | MO)")
250 DEALLOCATE (dist1, dist2, dist3)
251 END IF
252
253 CALL create_3c_tensor(t_3c_overl_int_gw_ri, dist1, dist2, dist3, pgrid_mo, &
254 sizes_ri_split, sizes_ao_split, sizes_mo, [1], [2, 3], name="(RI | AO MO)")
255 DEALLOCATE (dist1, dist2, dist3)
256
257 CALL create_3c_tensor(t_3c_overl_int_gw_ao, dist1, dist2, dist3, pgrid_mo, &
258 sizes_ao_split, sizes_ri_split, sizes_mo, [1], [2, 3], name="(AO | RI MO)")
259 DEALLOCATE (dist1, dist2, dist3)
260
261 CALL dbt_pgrid_destroy(pgrid_ao)
262 CALL dbt_pgrid_destroy(pgrid_mo)
263
264 CALL create_2c_tensor(mo_coeff_gw_t, dist1, dist2, pgrid_2d, sizes_ao_split, sizes_mo_1, name="(AO|MO)")
265 DEALLOCATE (dist1, dist2)
266 CALL dbt_pgrid_destroy(pgrid_2d)
267
268 CALL dbt_create(mat_mo_coeff_gw, mo_coeff_gw_t_tmp, name="MO coeffs")
269 CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw, mo_coeff_gw_t_tmp)
270
271 CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
272
273 bounds(1, 1) = homo - gw_corr_lev_occ + 1
274 bounds(2, 1) = homo + gw_corr_lev_virt
275
276 CALL dbt_get_info(t_3c_overl_int_ao_ao, nfull_total=bounds_3c)
277
278 ibounds(:, 1) = [1, bounds_3c(1)]
279 ibounds(:, 3) = [1, bounds_3c(3)]
280
281 cut_memory = SIZE(starts_array_mc)
282
283 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
284 DO i_mem = 1, cut_memory
285 CALL decompress_tensor(t_3c_overl_int(1, 1), t_3c_o_ind(1, 1, i_mem)%ind, t_3c_o_compressed(1, 1, i_mem), &
286 qs_env%mp2_env%ri_rpa_im_time%eps_compress)
287
288 ibounds(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
289
290 CALL dbt_copy(t_3c_overl_int(1, 1), t_3c_overl_int_ao_ao, move_data=.true.)
291
292 CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_ao_ao, 1.0_dp, &
293 t_3c_overl_int_ao_mo, contract_1=[1], notcontract_1=[2], &
294 contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
295 bounds_2=ibounds, move_data=.false., unit_nr=unit_nr_prv)
296
297 END DO
298 END IF
299
300 CALL cp_fm_release(fm_mat_mo_coeff_gw)
301
302 IF (do_ic_model) THEN
303 pdims = 0
304 CALL dbt_pgrid_create(para_env, pdims, pgrid_ic, &
305 tensor_dims=[SIZE(sizes_ri_split), nmo, nmo])
306
307 CALL create_3c_tensor(t_3c_overl_int_mo_ao, dist1, dist2, dist3, pgrid_ic, &
308 sizes_ri_split, sizes_mo_1, sizes_ao_split, [1, 2], [3], name="(RI MO | AO)")
309 DEALLOCATE (dist1, dist2, dist3)
310 CALL create_3c_tensor(t_3c_overl_int_mo_mo, dist1, dist2, dist3, pgrid_ic, &
311 sizes_ri_split, sizes_mo_1, sizes_mo_1, [1, 2], [3], name="(RI MO | MO)")
312 DEALLOCATE (dist1, dist2, dist3)
313 CALL dbt_create(t_3c_overl_int_mo_mo, t_3c_overl_nnp_ic)
314 CALL create_3c_tensor(t_3c_overl_nnp_ic_reflected, dist1, dist2, dist3, pgrid_ic, &
315 sizes_ri_split, sizes_mo_1, sizes_mo_1, [1], [2, 3], name="(RI | MO MO)")
316 DEALLOCATE (dist1, dist2, dist3)
317
318 CALL dbt_pgrid_destroy(pgrid_ic)
319
320 CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_mo_ao, order=[1, 3, 2])
321 CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_mo_ao, 0.0_dp, &
322 t_3c_overl_int_mo_mo, contract_1=[1], notcontract_1=[2], &
323 contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
324 bounds_2=bounds, move_data=.false., unit_nr=unit_nr_prv)
325 CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnp_ic)
326
327 NULLIFY (mat_mo_coeff_gw_reflected)
328 CALL dbcsr_init_p(mat_mo_coeff_gw_reflected)
329
330 CALL cp_dbcsr_m_by_n_from_row_template(mat_mo_coeff_gw_reflected, template=matrix_s(1)%matrix, n=nmo, &
331 sym=dbcsr_type_no_symmetry)
332
333 CALL reflect_mat_row(mat_mo_coeff_gw_reflected, mat_mo_coeff_gw, para_env, qs_env, unit_nr, do_alpha=my_do_alpha)
334
335 ! normalize reflected MOs (they are not properly normalized since high angular momentum basis functions
336 ! of the image molecule are not exactly reflected at the image plane (sign problem in p_z function)
337 CALL dbcsr_create(matrix=mat_work, template=mat_mo_coeff_gw_reflected, matrix_type=dbcsr_type_no_symmetry)
338
339 CALL dbcsr_get_info(mat_work, distribution=dist_templ, nblkcols_total=size_mo, col_blk_size=sizes_mo_blocked)
340
341 CALL dbcsr_distribution_get(dist_templ, nprows=nprows, npcols=npcols)
342
343 ALLOCATE (distp_1(size_mo), distp_2(size_mo))
344 CALL dbt_default_distvec(size_mo, nprows, sizes_mo_blocked, distp_1)
345 CALL dbt_default_distvec(size_mo, npcols, sizes_mo_blocked, distp_2)
346 CALL dbcsr_distribution_new(dist, template=dist_templ, row_dist=distp_1, col_dist=distp_2, reuse_arrays=.true.)
347
348 ALLOCATE (sizes_mo_p1(size_mo))
349 ALLOCATE (sizes_mo_p2(size_mo))
350 sizes_mo_p1(:) = sizes_mo_blocked
351 sizes_mo_p2(:) = sizes_mo_blocked
352 CALL dbcsr_create(mat_norm, "mo norm", dist, dbcsr_type_no_symmetry, sizes_mo_p1, sizes_mo_p2, &
353 reuse_arrays=.true.)
355
356 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s(1)%matrix, mat_mo_coeff_gw_reflected, 0.0_dp, mat_work)
357 CALL dbcsr_multiply("T", "N", 1.0_dp, mat_mo_coeff_gw_reflected, mat_work, 0.0_dp, mat_norm)
358
359 CALL dbcsr_release(mat_work)
360
361 ALLOCATE (norm(nmo))
362 norm = 0.0_dp
363
364 CALL dbcsr_get_diag(mat_norm, norm)
365 CALL para_env%sum(norm)
366
367 DO imo = bounds(1, 1), bounds(2, 1)
368 norm(imo) = 1.0_dp/sqrt(norm(imo))
369 END DO
370
371 CALL dbcsr_create(mat_norm_diag, template=mat_norm)
372 CALL dbcsr_release(mat_norm)
373
374 CALL dbcsr_add_on_diag(mat_norm_diag, 1.0_dp)
375
376 CALL dbcsr_set_diag(mat_norm_diag, norm)
377
378 CALL dbcsr_create(mat_mo_coeff_gw_reflected_norm, template=mat_mo_coeff_gw_reflected)
379 CALL dbcsr_multiply("N", "N", 1.0_dp, mat_mo_coeff_gw_reflected, mat_norm_diag, 0.0_dp, mat_mo_coeff_gw_reflected_norm)
380 CALL dbcsr_release(mat_norm_diag)
381
382 CALL dbcsr_filter(mat_mo_coeff_gw_reflected_norm, 1.0e-20_dp)
383
384 CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw_reflected_norm, mo_coeff_gw_t_tmp)
385 CALL dbcsr_release(mat_mo_coeff_gw_reflected_norm)
386 CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
387
388 CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_ao_ao, 0.0_dp, &
389 t_3c_overl_int_ao_mo, contract_1=[1], notcontract_1=[2], &
390 contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
391 bounds_2=bounds, move_data=.false., unit_nr=unit_nr_prv)
392
393 CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_mo_ao, order=[1, 3, 2])
394 CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_mo_ao, 0.0_dp, &
395 t_3c_overl_int_mo_mo, contract_1=[1], notcontract_1=[2], &
396 contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
397 bounds_2=bounds, move_data=.false., unit_nr=unit_nr_prv)
398 CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnp_ic_reflected)
399 CALL dbt_destroy(t_3c_overl_int_mo_ao)
400 CALL dbt_destroy(t_3c_overl_int_mo_mo)
401
402 CALL dbcsr_release_p(mat_mo_coeff_gw_reflected)
403
404 END IF
405
406 IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
407 CALL alloc_containers(t_3c_o_mo_compressed, 1)
408 CALL get_tensor_occupancy(t_3c_overl_int_ao_mo, nze, occ)
409 memory_3c = 0.0_dp
410
411 CALL compress_tensor(t_3c_overl_int_ao_mo, t_3c_o_mo_ind, t_3c_o_mo_compressed, &
412 qs_env%mp2_env%ri_rpa_im_time%eps_compress, memory_3c)
413
414 CALL para_env%sum(memory_3c)
415 compression_factor = real(nze, dp)*1.0e-06_dp*8.0_dp/memory_3c
416
417 IF (unit_nr > 0) THEN
418 WRITE (unit=unit_nr, fmt="((T3,A,T66,F11.2,A4))") &
419 "MEMORY_INFO| Memory of MO-contracted tensor (compressed):", memory_3c, ' MiB'
420
421 WRITE (unit=unit_nr, fmt="((T3,A,T60,F21.2))") &
422 "MEMORY_INFO| Compression factor: ", compression_factor
423 END IF
424 END IF
425
426 CALL dbcsr_release_p(mat_mo_coeff_gw)
427
428 CALL dbt_destroy(t_3c_overl_int_ao_ao)
429 CALL dbt_destroy(mo_coeff_gw_t)
430 CALL dbt_destroy(mo_coeff_gw_t_tmp)
431
432 CALL timestop(handle)
433
434 END SUBROUTINE get_tensor_3c_overl_int_gw
435
436! **************************************************************************************************
437!> \brief reflect from V = (A,B|B,A) to V_reflected = (B,A|A,B) where A belongs to the block of the molecule
438!> and B to the off diagonal block between molecule and image of the molecule
439!> \param mat_reflected ...
440!> \param mat_orig ...
441!> \param para_env ...
442!> \param qs_env ...
443!> \param unit_nr ...
444!> \param do_alpha ...
445! **************************************************************************************************
446 SUBROUTINE reflect_mat_row(mat_reflected, mat_orig, para_env, qs_env, unit_nr, do_alpha)
447 TYPE(dbcsr_type), INTENT(INOUT) :: mat_reflected
448 TYPE(dbcsr_type), INTENT(IN) :: mat_orig
449 TYPE(mp_para_env_type), INTENT(IN) :: para_env
450 TYPE(qs_environment_type), POINTER :: qs_env
451 INTEGER, INTENT(IN) :: unit_nr
452 LOGICAL, INTENT(IN) :: do_alpha
453
454 CHARACTER(LEN=*), PARAMETER :: routinen = 'reflect_mat_row'
455
456 INTEGER :: handle, i_atom, j_atom, natom, &
457 nblkcols_total, nblkrows_total
458 INTEGER, ALLOCATABLE, DIMENSION(:) :: image_atom
459 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
460 LOGICAL :: found_image_atom
461 REAL(kind=dp) :: avg_z_dist, delta, eps_dist2, &
462 min_z_dist, ra(3), rb(3), sum_z, &
463 z_reflection
464 TYPE(cell_type), POINTER :: cell
465 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
466
467 CALL timeset(routinen, handle)
468
469 CALL dbcsr_reserve_all_blocks(mat_reflected)
470
471 CALL get_qs_env(qs_env, cell=cell, &
472 particle_set=particle_set)
473
474 ! first check, whether we have an image molecule
475 CALL dbcsr_get_info(mat_orig, &
476 nblkrows_total=nblkrows_total, &
477 nblkcols_total=nblkcols_total, &
478 row_blk_size=row_blk_sizes, &
479 col_blk_size=col_blk_sizes)
480
481 natom = SIZE(particle_set)
482 cpassert(natom == nblkrows_total)
483
484 eps_dist2 = qs_env%mp2_env%ri_g0w0%eps_dist
485 eps_dist2 = eps_dist2*eps_dist2
486
487 sum_z = 0.0_dp
488
489 DO i_atom = 1, natom
490
491 ra(:) = pbc(particle_set(i_atom)%r, cell)
492
493 sum_z = sum_z + ra(3)
494
495 END DO
496
497 z_reflection = sum_z/real(natom, kind=dp)
498
499 sum_z = 0.0_dp
500
501 DO i_atom = 1, natom
502
503 ra(:) = pbc(particle_set(i_atom)%r, cell)
504
505 sum_z = sum_z + abs(ra(3) - z_reflection)
506
507 END DO
508
509 avg_z_dist = sum_z/real(natom, kind=dp)
510
511 min_z_dist = avg_z_dist
512
513 DO i_atom = 1, natom
514
515 ra(:) = pbc(particle_set(i_atom)%r, cell)
516
517 IF (abs(ra(3) - z_reflection) < min_z_dist) THEN
518 min_z_dist = abs(ra(3) - z_reflection)
519 END IF
520
521 END DO
522
523 IF (unit_nr > 0 .AND. do_alpha) THEN
524 WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Average distance of the molecule to the image plane:', &
525 avg_z_dist*0.529_dp, ' A'
526 WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Minimum distance of the molecule to the image plane:', &
527 min_z_dist*0.529_dp, ' A'
528 END IF
529
530 ALLOCATE (image_atom(nblkrows_total))
531 image_atom = 0
532
533 DO i_atom = 1, natom
534
535 found_image_atom = .false.
536
537 ra(:) = pbc(particle_set(i_atom)%r, cell)
538
539 DO j_atom = 1, natom
540
541 rb(:) = pbc(particle_set(j_atom)%r, cell)
542
543 delta = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) + rb(3) - 2.0_dp*z_reflection)**2
544
545 ! SQRT(delta) < eps_dist
546 IF (delta < eps_dist2) THEN
547 ! this CPASSERT ensures that there is at most one image atom for each atom
548 cpassert(.NOT. found_image_atom)
549 image_atom(i_atom) = j_atom
550 found_image_atom = .true.
551 ! check whether we have the same basis at the image atom
552 ! if this is wrong, check whether you have the same basis sets for the molecule and the image
553 cpassert(row_blk_sizes(i_atom) == row_blk_sizes(j_atom))
554 END IF
555
556 END DO
557
558 ! this CPASSERT ensures that there is at least one image atom for each atom
559 cpassert(found_image_atom)
560
561 END DO
562
563 CALL cp_dbcsr_from_subgroups(mat_orig, mat_reflected, para_env, row_mapping=image_atom, copy_only=.true.)
564
565 CALL timestop(handle)
566
567 END SUBROUTINE reflect_mat_row
568
569! **************************************************************************************************
570!> \brief ...
571!> \param qs_env ...
572!> \param atom_from_basis_index ...
573!> \param basis_size ...
574!> \param basis_type ...
575!> \param first_bf_from_atom ...
576! **************************************************************************************************
577 SUBROUTINE get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, &
578 basis_type, first_bf_from_atom)
579 TYPE(qs_environment_type), POINTER :: qs_env
580 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_basis_index
581 INTEGER :: basis_size
582 CHARACTER(LEN=*) :: basis_type
583 INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL :: first_bf_from_atom
584
585 INTEGER :: iatom, lll, natom, nkind
586 INTEGER, DIMENSION(:), POINTER :: row_blk_end, row_blk_start
587 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
588 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
589 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
590
591 NULLIFY (qs_kind_set, particle_set)
592 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, nkind=nkind, &
593 particle_set=particle_set)
594
595 ALLOCATE (row_blk_start(natom))
596 ALLOCATE (row_blk_end(natom))
597 ALLOCATE (basis_set(nkind))
598 CALL basis_set_list_setup(basis_set, basis_type, qs_kind_set)
599 CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, &
600 basis=basis_set)
601 DO lll = 1, basis_size
602 DO iatom = 1, natom
603 IF (lll >= row_blk_start(iatom) .AND. lll <= row_blk_end(iatom)) THEN
604 atom_from_basis_index(lll) = iatom
605 END IF
606 END DO
607 END DO
608
609 IF (PRESENT(first_bf_from_atom)) first_bf_from_atom(1:natom) = row_blk_start(:)
610
611 DEALLOCATE (basis_set)
612 DEALLOCATE (row_blk_start)
613 DEALLOCATE (row_blk_end)
614
616
617! **************************************************************************************************
618!> \brief ...
619!> \param weight_re ...
620!> \param weight_im ...
621!> \param num_cells ...
622!> \param iatom ...
623!> \param jatom ...
624!> \param xkp ...
625!> \param wkp_W ...
626!> \param cell ...
627!> \param index_to_cell ...
628!> \param hmat ...
629!> \param particle_set ...
630! **************************************************************************************************
631 SUBROUTINE compute_weight_re_im(weight_re, weight_im, &
632 num_cells, iatom, jatom, xkp, wkp_W, &
633 cell, index_to_cell, hmat, particle_set)
634
635 REAL(kind=dp) :: weight_re, weight_im
636 INTEGER :: num_cells, iatom, jatom
637 REAL(kind=dp), DIMENSION(3) :: xkp
638 REAL(kind=dp) :: wkp_w
639 TYPE(cell_type), POINTER :: cell
640 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
641 REAL(kind=dp), DIMENSION(3, 3) :: hmat
642 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
643
644 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_weight_re_im'
645
646 INTEGER :: handle, icell, n_equidistant_cells, &
647 xcell, ycell, zcell
648 REAL(kind=dp) :: abs_rab_cell, abs_rab_cell_min, arg
649 REAL(kind=dp), DIMENSION(3) :: cell_vector, rab_cell_i
650
651 CALL timeset(routinen, handle)
652
653 weight_re = 0.0_dp
654 weight_im = 0.0_dp
655
656 abs_rab_cell_min = 1.0e10_dp
657
658 n_equidistant_cells = 0
659
660 DO icell = 1, num_cells
661
662 xcell = index_to_cell(1, icell)
663 ycell = index_to_cell(2, icell)
664 zcell = index_to_cell(3, icell)
665
666 cell_vector(1:3) = matmul(hmat, real([xcell, ycell, zcell], dp))
667
668 rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - &
669 (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
670
671 abs_rab_cell = sqrt(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
672
673 IF (abs_rab_cell < abs_rab_cell_min) THEN
674 abs_rab_cell_min = abs_rab_cell
675 END IF
676
677 END DO
678
679 DO icell = 1, num_cells
680
681 xcell = index_to_cell(1, icell)
682 ycell = index_to_cell(2, icell)
683 zcell = index_to_cell(3, icell)
684
685 cell_vector(1:3) = matmul(hmat, real([xcell, ycell, zcell], dp))
686
687 rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - &
688 (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
689
690 abs_rab_cell = sqrt(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
691
692 IF (abs_rab_cell < abs_rab_cell_min + 0.1_dp) THEN
693
694 arg = real(xcell, dp)*xkp(1) + real(ycell, dp)*xkp(2) + real(zcell, dp)*xkp(3)
695
696 weight_re = weight_re + wkp_w*cos(twopi*arg)
697 weight_im = weight_im + wkp_w*sin(twopi*arg)
698
699 n_equidistant_cells = n_equidistant_cells + 1
700
701 END IF
702
703 END DO
704
705 weight_re = weight_re/real(n_equidistant_cells, kind=dp)
706 weight_im = weight_im/real(n_equidistant_cells, kind=dp)
707
708 CALL timestop(handle)
709
710 END SUBROUTINE compute_weight_re_im
711
712END MODULE rpa_gw_im_time_util
Define the atomic kind types and their sub types.
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public dbcsr_distribution_release(dist)
...
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
subroutine, public dbcsr_set_diag(matrix, diag)
Copies the diagonal elements from the given array into the given matrix.
subroutine, public dbcsr_get_diag(matrix, diag)
Copies the diagonal elements from the given matrix into the given array.
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Redistribution of DBCSR blocks and distributed matrix slices.
subroutine, public cp_dbcsr_from_subgroups(mat_local, mat_global, para_env, row_mapping, copy_only)
Sum subgroup blocks into the global layout, or copy mapped blocks.
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
Types and set/get functions for HFX.
Definition hfx_types.F:16
subroutine, public alloc_containers(data, bin_size)
...
Definition hfx_types.F:2989
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public twopi
Interface to the message passing library MPI.
subroutine, public mp_dims_create(nodes, dims)
wrapper to MPI_Dims_create
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
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, mimic, 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, matrix_vhxc, 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, xcint_weights, 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.
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Define the quickstep kind type and their sub types.
Utility methods to build 3-center integral tensors of various types.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
subroutine, public split_block_sizes(blk_sizes, blk_sizes_split, max_size)
...
subroutine, public pgf_block_sizes(atomic_kind_set, basis, min_blk_size, pgf_blk_sizes)
...
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
subroutine, public compress_tensor(tensor, blk_indices, compressed, eps, memory)
...
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Utility routines for GW with imaginary time.
subroutine, public compute_weight_re_im(weight_re, weight_im, num_cells, iatom, jatom, xkp, wkp_w, cell, index_to_cell, hmat, particle_set)
...
subroutine, public get_tensor_3c_overl_int_gw(t_3c_overl_int, t_3c_o_compressed, t_3c_o_ind, t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, starts_array_mc, ends_array_mc, mo_coeff, matrix_s, gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, para_env, do_ic_model, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, qs_env, unit_nr, do_alpha)
...
subroutine, public get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type, first_bf_from_atom)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a full matrix
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.