Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utility routines for GW with imaginary time
10 : !> \par History
11 : !> 06.2019 split from rpa_im_time.F [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE rpa_gw_im_time_util
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE basis_set_types, ONLY: gto_basis_set_p_type
17 : USE cell_types, ONLY: cell_type,&
18 : pbc
19 : USE cp_dbcsr_api, ONLY: &
20 : dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, dbcsr_distribution_release, &
21 : dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, &
22 : dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_type, dbcsr_type_no_symmetry
23 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
24 : dbcsr_get_diag,&
25 : dbcsr_reserve_all_blocks,&
26 : dbcsr_set_diag
27 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
28 : cp_dbcsr_m_by_n_from_row_template
29 : USE cp_dbcsr_redistribution, ONLY: cp_dbcsr_from_subgroups
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_release,&
32 : cp_fm_set_element,&
33 : cp_fm_to_fm,&
34 : cp_fm_type
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,&
39 : block_ind_type,&
40 : hfx_compression_type
41 : USE kinds, ONLY: dp,&
42 : int_8
43 : USE mathconstants, ONLY: twopi
44 : USE message_passing, ONLY: mp_dims_create,&
45 : mp_para_env_type
46 : USE particle_methods, ONLY: get_particle_set
47 : USE particle_types, ONLY: particle_type
48 : USE qs_environment_types, ONLY: get_qs_env,&
49 : qs_environment_type
50 : USE qs_integral_utils, ONLY: basis_set_list_setup
51 : USE qs_kind_types, ONLY: qs_kind_type
52 : USE qs_tensors, ONLY: compress_tensor,&
53 : decompress_tensor,&
54 : get_tensor_occupancy
55 : USE qs_tensors_types, ONLY: create_2c_tensor,&
56 : create_3c_tensor,&
57 : pgf_block_sizes,&
58 : split_block_sizes
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 :
67 : PUBLIC :: get_tensor_3c_overl_int_gw, compute_weight_re_im, get_atom_index_from_basis_function_index
68 :
69 : CONTAINS
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 54 : 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 108 : 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 108 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist1, dist2, dist3, sizes_AO, &
135 54 : sizes_AO_split, sizes_MO, sizes_MO_1, &
136 108 : 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 108 : INTEGER, DIMENSION(:), POINTER :: distp_1, distp_2, sizes_MO_blocked, &
142 54 : sizes_MO_p1, sizes_MO_p2
143 : LOGICAL :: memory_info, my_do_alpha
144 : REAL(dp) :: compression_factor, memory_3c, occ
145 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: norm
146 54 : 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 486 : TYPE(dbt_pgrid_type) :: pgrid_2d, pgrid_AO, pgrid_ic, pgrid_MO
154 702 : TYPE(dbt_type) :: mo_coeff_gw_t, mo_coeff_gw_t_tmp, &
155 378 : t_3c_overl_int_ao_ao, &
156 378 : t_3c_overl_int_mo_ao, &
157 378 : t_3c_overl_int_mo_mo
158 54 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ao, basis_set_ri_aux
159 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
160 54 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
161 :
162 54 : memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info
163 54 : IF (memory_info) THEN
164 0 : unit_nr_prv = unit_nr
165 : ELSE
166 54 : unit_nr_prv = 0
167 : END IF
168 :
169 54 : my_do_alpha = .FALSE.
170 54 : IF (PRESENT(do_alpha)) my_do_alpha = do_alpha
171 :
172 54 : CALL timeset(routineN, handle)
173 :
174 54 : 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 54 : CALL cp_fm_create(fm_mat_mo_coeff_gw, mo_coeff%matrix_struct)
177 54 : CALL cp_fm_to_fm(mo_coeff, fm_mat_mo_coeff_gw)
178 :
179 : ! set MO coeffs to zero where
180 1396 : DO irow_global = 1, nmo
181 3094 : DO icol_global = 1, homo - gw_corr_lev_occ
182 3094 : CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
183 : END DO
184 32506 : DO icol_global = homo + gw_corr_lev_virt + 1, nmo
185 32452 : 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 54 : NULLIFY (mat_mo_coeff_gw)
190 54 : 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 54 : sym=dbcsr_type_no_symmetry)
194 :
195 : CALL copy_fm_to_dbcsr(fm_mat_mo_coeff_gw, &
196 : mat_mo_coeff_gw, &
197 54 : keep_sparsity=.FALSE.)
198 :
199 : ! just remove the blocks which have been set to zero
200 54 : CALL dbcsr_filter(mat_mo_coeff_gw, 1.0E-20_dp)
201 :
202 54 : min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize
203 54 : min_bsize_mo = qs_env%mp2_env%ri_rpa_im_time%min_bsize_mo
204 :
205 108 : CALL split_block_sizes([gw_corr_lev_occ + gw_corr_lev_virt], sizes_MO, min_bsize_mo)
206 162 : ALLOCATE (sizes_MO_1(nmo))
207 1396 : sizes_MO_1(:) = 1
208 :
209 54 : nmo_blk_gw = SIZE(sizes_MO)
210 54 : CALL move_alloc(sizes_MO, tmp)
211 162 : ALLOCATE (sizes_MO(nmo_blk_gw + 2))
212 54 : sizes_MO(1) = homo - gw_corr_lev_occ
213 108 : sizes_MO(2:SIZE(tmp) + 1) = tmp(:)
214 54 : sizes_MO(SIZE(tmp) + 2) = nmo - (homo + gw_corr_lev_virt)
215 :
216 432 : ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
217 54 : CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
218 54 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux)
219 54 : CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
220 54 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao)
221 :
222 54 : CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_AO_split)
223 54 : CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_RI_split)
224 :
225 54 : DEALLOCATE (basis_set_ao, basis_set_ri_aux)
226 :
227 54 : pdims = 0
228 : CALL dbt_pgrid_create(para_env, pdims, pgrid_AO, &
229 216 : tensor_dims=[SIZE(sizes_RI_split), SIZE(sizes_AO_split), SIZE(sizes_AO_split)])
230 :
231 54 : pdims_2d = 0
232 54 : 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 216 : pdims = [pdims_2d(1), pdims_2d(2), 1]
236 : CALL dbt_pgrid_create(para_env, pdims, pgrid_MO, &
237 216 : tensor_dims=[SIZE(sizes_RI_split), SIZE(sizes_AO_split), 1])
238 :
239 54 : pdims_2d = 0
240 : CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d, &
241 162 : 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 54 : sizes_RI_split, sizes_AO_split, sizes_AO_split, [1, 2], [3], name="(RI AO | AO)")
245 54 : DEALLOCATE (dist1, dist2, dist3)
246 :
247 54 : 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 36 : sizes_RI_split, sizes_AO_split, sizes_MO_1, [1, 2], [3], name="(RI AO | MO)")
250 36 : 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 54 : sizes_RI_split, sizes_AO_split, sizes_MO, [1], [2, 3], name="(RI | AO MO)")
255 54 : DEALLOCATE (dist1, dist2, dist3)
256 :
257 : CALL create_3c_tensor(t_3c_overl_int_gw_AO, dist1, dist2, dist3, pgrid_MO, &
258 54 : sizes_AO_split, sizes_RI_split, sizes_MO, [1], [2, 3], name="(AO | RI MO)")
259 54 : DEALLOCATE (dist1, dist2, dist3)
260 :
261 54 : CALL dbt_pgrid_destroy(pgrid_AO)
262 54 : CALL dbt_pgrid_destroy(pgrid_MO)
263 :
264 54 : CALL create_2c_tensor(mo_coeff_gw_t, dist1, dist2, pgrid_2d, sizes_AO_split, sizes_MO_1, name="(AO|MO)")
265 54 : DEALLOCATE (dist1, dist2)
266 54 : CALL dbt_pgrid_destroy(pgrid_2d)
267 :
268 54 : CALL dbt_create(mat_mo_coeff_gw, mo_coeff_gw_t_tmp, name="MO coeffs")
269 54 : CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw, mo_coeff_gw_t_tmp)
270 :
271 54 : CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
272 :
273 54 : bounds(1, 1) = homo - gw_corr_lev_occ + 1
274 54 : bounds(2, 1) = homo + gw_corr_lev_virt
275 :
276 54 : CALL dbt_get_info(t_3c_overl_int_ao_ao, nfull_total=bounds_3c)
277 :
278 162 : ibounds(:, 1) = [1, bounds_3c(1)]
279 162 : ibounds(:, 3) = [1, bounds_3c(3)]
280 :
281 54 : cut_memory = SIZE(starts_array_mc)
282 :
283 54 : IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
284 72 : 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 36 : qs_env%mp2_env%ri_rpa_im_time%eps_compress)
287 :
288 108 : ibounds(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
289 :
290 36 : 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 72 : bounds_2=ibounds, move_data=.FALSE., unit_nr=unit_nr_prv)
296 :
297 : END DO
298 : END IF
299 :
300 54 : CALL cp_fm_release(fm_mat_mo_coeff_gw)
301 :
302 54 : IF (do_ic_model) THEN
303 2 : pdims = 0
304 : CALL dbt_pgrid_create(para_env, pdims, pgrid_ic, &
305 8 : 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 2 : sizes_RI_split, sizes_MO_1, sizes_AO_split, [1, 2], [3], name="(RI MO | AO)")
309 2 : DEALLOCATE (dist1, dist2, dist3)
310 : CALL create_3c_tensor(t_3c_overl_int_mo_mo, dist1, dist2, dist3, pgrid_ic, &
311 2 : sizes_RI_split, sizes_MO_1, sizes_MO_1, [1, 2], [3], name="(RI MO | MO)")
312 2 : DEALLOCATE (dist1, dist2, dist3)
313 2 : 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 2 : sizes_RI_split, sizes_MO_1, sizes_MO_1, [1], [2, 3], name="(RI | MO MO)")
316 2 : DEALLOCATE (dist1, dist2, dist3)
317 :
318 2 : CALL dbt_pgrid_destroy(pgrid_ic)
319 :
320 2 : 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 2 : bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
325 2 : CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic)
326 :
327 2 : NULLIFY (mat_mo_coeff_gw_reflected)
328 2 : 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 2 : sym=dbcsr_type_no_symmetry)
332 :
333 2 : 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 2 : CALL dbcsr_create(matrix=mat_work, template=mat_mo_coeff_gw_reflected, matrix_type=dbcsr_type_no_symmetry)
338 :
339 2 : CALL dbcsr_get_info(mat_work, distribution=dist_templ, nblkcols_total=size_MO, col_blk_size=sizes_MO_blocked)
340 :
341 2 : CALL dbcsr_distribution_get(dist_templ, nprows=nprows, npcols=npcols)
342 :
343 8 : ALLOCATE (distp_1(size_MO), distp_2(size_MO))
344 2 : CALL dbt_default_distvec(size_MO, nprows, sizes_MO_blocked, distp_1)
345 2 : CALL dbt_default_distvec(size_MO, npcols, sizes_MO_blocked, distp_2)
346 2 : CALL dbcsr_distribution_new(dist, template=dist_templ, row_dist=distp_1, col_dist=distp_2, reuse_arrays=.TRUE.)
347 :
348 4 : ALLOCATE (sizes_MO_p1(size_MO))
349 4 : ALLOCATE (sizes_MO_p2(size_MO))
350 14 : sizes_MO_p1(:) = sizes_MO_blocked
351 14 : 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 2 : reuse_arrays=.TRUE.)
354 2 : CALL dbcsr_distribution_release(dist)
355 :
356 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s(1)%matrix, mat_mo_coeff_gw_reflected, 0.0_dp, mat_work)
357 2 : CALL dbcsr_multiply("T", "N", 1.0_dp, mat_mo_coeff_gw_reflected, mat_work, 0.0_dp, mat_norm)
358 :
359 2 : CALL dbcsr_release(mat_work)
360 :
361 6 : ALLOCATE (norm(nmo))
362 2 : norm = 0.0_dp
363 :
364 2 : CALL dbcsr_get_diag(mat_norm, norm)
365 2 : CALL para_env%sum(norm)
366 :
367 14 : DO imo = bounds(1, 1), bounds(2, 1)
368 14 : norm(imo) = 1.0_dp/SQRT(norm(imo))
369 : END DO
370 :
371 2 : CALL dbcsr_create(mat_norm_diag, template=mat_norm)
372 2 : CALL dbcsr_release(mat_norm)
373 :
374 2 : CALL dbcsr_add_on_diag(mat_norm_diag, 1.0_dp)
375 :
376 2 : CALL dbcsr_set_diag(mat_norm_diag, norm)
377 :
378 2 : CALL dbcsr_create(mat_mo_coeff_gw_reflected_norm, template=mat_mo_coeff_gw_reflected)
379 2 : 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 2 : CALL dbcsr_release(mat_norm_diag)
381 :
382 2 : CALL dbcsr_filter(mat_mo_coeff_gw_reflected_norm, 1.0E-20_dp)
383 :
384 2 : CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw_reflected_norm, mo_coeff_gw_t_tmp)
385 2 : CALL dbcsr_release(mat_mo_coeff_gw_reflected_norm)
386 2 : 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 2 : bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
392 :
393 2 : 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 2 : bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
398 2 : CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic_reflected)
399 2 : CALL dbt_destroy(t_3c_overl_int_mo_ao)
400 2 : CALL dbt_destroy(t_3c_overl_int_mo_mo)
401 :
402 4 : CALL dbcsr_release_p(mat_mo_coeff_gw_reflected)
403 :
404 : END IF
405 :
406 54 : IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
407 36 : CALL alloc_containers(t_3c_O_mo_compressed, 1)
408 36 : CALL get_tensor_occupancy(t_3c_overl_int_ao_mo, nze, occ)
409 36 : 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 36 : qs_env%mp2_env%ri_rpa_im_time%eps_compress, memory_3c)
413 :
414 36 : CALL para_env%sum(memory_3c)
415 36 : compression_factor = REAL(nze, dp)*1.0E-06_dp*8.0_dp/memory_3c
416 :
417 36 : IF (unit_nr > 0) THEN
418 : WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") &
419 18 : "MEMORY_INFO| Memory of MO-contracted tensor (compressed):", memory_3c, ' MiB'
420 :
421 : WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") &
422 18 : "MEMORY_INFO| Compression factor: ", compression_factor
423 : END IF
424 : END IF
425 :
426 54 : CALL dbcsr_release_p(mat_mo_coeff_gw)
427 :
428 54 : CALL dbt_destroy(t_3c_overl_int_ao_ao)
429 54 : CALL dbt_destroy(mo_coeff_gw_t)
430 54 : CALL dbt_destroy(mo_coeff_gw_t_tmp)
431 :
432 54 : CALL timestop(handle)
433 :
434 324 : 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 2 : 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 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: image_atom
459 2 : 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 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
466 :
467 2 : CALL timeset(routineN, handle)
468 :
469 2 : CALL dbcsr_reserve_all_blocks(mat_reflected)
470 :
471 : CALL get_qs_env(qs_env, cell=cell, &
472 2 : 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 2 : col_blk_size=col_blk_sizes)
480 :
481 2 : natom = SIZE(particle_set)
482 2 : CPASSERT(natom == nblkrows_total)
483 :
484 2 : eps_dist2 = qs_env%mp2_env%ri_g0w0%eps_dist
485 2 : eps_dist2 = eps_dist2*eps_dist2
486 :
487 2 : sum_z = 0.0_dp
488 :
489 18 : DO i_atom = 1, natom
490 :
491 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
492 :
493 18 : sum_z = sum_z + ra(3)
494 :
495 : END DO
496 :
497 2 : z_reflection = sum_z/REAL(natom, KIND=dp)
498 :
499 2 : sum_z = 0.0_dp
500 :
501 18 : DO i_atom = 1, natom
502 :
503 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
504 :
505 18 : sum_z = sum_z + ABS(ra(3) - z_reflection)
506 :
507 : END DO
508 :
509 2 : avg_z_dist = sum_z/REAL(natom, KIND=dp)
510 :
511 2 : min_z_dist = avg_z_dist
512 :
513 18 : DO i_atom = 1, natom
514 :
515 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
516 :
517 18 : 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 2 : IF (unit_nr > 0 .AND. do_alpha) THEN
524 1 : WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Average distance of the molecule to the image plane:', &
525 2 : avg_z_dist*0.529_dp, ' A'
526 1 : WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Minimum distance of the molecule to the image plane:', &
527 2 : min_z_dist*0.529_dp, ' A'
528 : END IF
529 :
530 6 : ALLOCATE (image_atom(nblkrows_total))
531 2 : image_atom = 0
532 :
533 18 : DO i_atom = 1, natom
534 :
535 16 : found_image_atom = .FALSE.
536 :
537 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
538 :
539 144 : DO j_atom = 1, natom
540 :
541 128 : rb(:) = pbc(particle_set(j_atom)%r, cell)
542 :
543 128 : 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 144 : IF (delta < eps_dist2) THEN
547 : ! this CPASSERT ensures that there is at most one image atom for each atom
548 16 : CPASSERT(.NOT. found_image_atom)
549 16 : image_atom(i_atom) = j_atom
550 16 : 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 16 : 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 18 : CPASSERT(found_image_atom)
560 :
561 : END DO
562 :
563 2 : CALL cp_dbcsr_from_subgroups(mat_orig, mat_reflected, para_env, row_mapping=image_atom, copy_only=.TRUE.)
564 :
565 2 : CALL timestop(handle)
566 :
567 6 : 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 14898 : 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 14898 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
588 14898 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
589 14898 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
590 :
591 14898 : NULLIFY (qs_kind_set, particle_set)
592 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, nkind=nkind, &
593 14898 : particle_set=particle_set)
594 :
595 44694 : ALLOCATE (row_blk_start(natom))
596 29796 : ALLOCATE (row_blk_end(natom))
597 68512 : ALLOCATE (basis_set(nkind))
598 14898 : 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 14898 : basis=basis_set)
601 399746 : DO LLL = 1, basis_size
602 1401550 : DO iatom = 1, natom
603 1386652 : IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN
604 384848 : atom_from_basis_index(LLL) = iatom
605 : END IF
606 : END DO
607 : END DO
608 :
609 14898 : IF (PRESENT(first_bf_from_atom)) first_bf_from_atom(1:natom) = row_blk_start(:)
610 :
611 14898 : DEALLOCATE (basis_set)
612 14898 : DEALLOCATE (row_blk_start)
613 14898 : DEALLOCATE (row_blk_end)
614 :
615 14898 : END SUBROUTINE get_atom_index_from_basis_function_index
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 409772 : 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 409772 : CALL timeset(routineN, handle)
652 :
653 409772 : weight_re = 0.0_dp
654 409772 : weight_im = 0.0_dp
655 :
656 409772 : abs_rab_cell_min = 1.0E10_dp
657 :
658 409772 : n_equidistant_cells = 0
659 :
660 3087532 : DO icell = 1, num_cells
661 :
662 2677760 : xcell = index_to_cell(1, icell)
663 2677760 : ycell = index_to_cell(2, icell)
664 2677760 : zcell = index_to_cell(3, icell)
665 :
666 42844160 : 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 10711040 : (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
670 :
671 2677760 : abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
672 :
673 409772 : 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 3087532 : DO icell = 1, num_cells
680 :
681 2677760 : xcell = index_to_cell(1, icell)
682 2677760 : ycell = index_to_cell(2, icell)
683 2677760 : zcell = index_to_cell(3, icell)
684 :
685 42844160 : 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 10711040 : (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
689 :
690 2677760 : abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
691 :
692 3087532 : IF (abs_rab_cell < abs_rab_cell_min + 0.1_dp) THEN
693 :
694 409932 : arg = REAL(xcell, dp)*xkp(1) + REAL(ycell, dp)*xkp(2) + REAL(zcell, dp)*xkp(3)
695 :
696 409932 : weight_re = weight_re + wkp_W*COS(twopi*arg)
697 409932 : weight_im = weight_im + wkp_W*SIN(twopi*arg)
698 :
699 409932 : n_equidistant_cells = n_equidistant_cells + 1
700 :
701 : END IF
702 :
703 : END DO
704 :
705 409772 : weight_re = weight_re/REAL(n_equidistant_cells, KIND=dp)
706 409772 : weight_im = weight_im/REAL(n_equidistant_cells, KIND=dp)
707 :
708 409772 : CALL timestop(handle)
709 :
710 409772 : END SUBROUTINE compute_weight_re_im
711 :
712 : END MODULE rpa_gw_im_time_util
|
