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 RPA and 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_cfm_basic_linalg, ONLY: cp_cfm_column_scale
20 : USE cp_cfm_types, ONLY: cp_cfm_get_info,&
21 : cp_cfm_set_all,&
22 : cp_cfm_type
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, dbcsr_distribution_release, &
25 : dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, &
26 : dbcsr_p_type, dbcsr_release, dbcsr_release_p, dbcsr_type, dbcsr_type_no_symmetry
27 : USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
28 : dbcsr_get_diag,&
29 : dbcsr_reserve_all_blocks,&
30 : dbcsr_set_diag
31 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
32 : cp_dbcsr_m_by_n_from_row_template
33 : USE cp_dbcsr_redistribution, ONLY: cp_dbcsr_from_subgroups
34 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_info,&
37 : cp_fm_release,&
38 : cp_fm_set_all,&
39 : cp_fm_set_element,&
40 : cp_fm_to_fm,&
41 : cp_fm_type
42 : USE dbt_api, ONLY: &
43 : dbt_contract, dbt_copy, dbt_copy_matrix_to_tensor, dbt_create, dbt_default_distvec, &
44 : dbt_destroy, dbt_get_info, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
45 : USE hfx_types, ONLY: alloc_containers,&
46 : block_ind_type,&
47 : hfx_compression_type
48 : USE kinds, ONLY: dp,&
49 : int_8
50 : USE mathconstants, ONLY: twopi,&
51 : z_one,&
52 : z_zero
53 : USE message_passing, ONLY: mp_dims_create,&
54 : mp_para_env_type
55 : USE parallel_gemm_api, ONLY: parallel_gemm
56 : USE particle_methods, ONLY: get_particle_set
57 : USE particle_types, ONLY: particle_type
58 : USE qs_environment_types, ONLY: get_qs_env,&
59 : qs_environment_type
60 : USE qs_integral_utils, ONLY: basis_set_list_setup
61 : USE qs_kind_types, ONLY: qs_kind_type
62 : USE qs_tensors, ONLY: compress_tensor,&
63 : decompress_tensor,&
64 : get_tensor_occupancy
65 : USE qs_tensors_types, ONLY: create_2c_tensor,&
66 : create_3c_tensor,&
67 : pgf_block_sizes,&
68 : split_block_sizes
69 : #include "./base/base_uses.f90"
70 :
71 : IMPLICIT NONE
72 :
73 : PRIVATE
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_im_time_util'
76 :
77 : PUBLIC :: get_tensor_3c_overl_int_gw, compute_weight_re_im, get_atom_index_from_basis_function_index
78 : PUBLIC :: build_ao_propagator
79 :
80 : INTERFACE build_ao_propagator
81 : MODULE PROCEDURE build_ao_propagator_fm, build_ao_propagator_cfm
82 : END INTERFACE build_ao_propagator
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Builds G = X X^T, where X is the selected, column-scaled MO window.
88 : !> Only the selected workspace columns are overwritten; coefficients are unchanged.
89 : !> \param coefficients AO-by-MO coefficient matrix
90 : !> \param column_factors factors for consecutive MOs starting at first_mo
91 : !> \param first_mo first global MO column in the window
92 : !> \param work coefficient workspace with the same shape and distribution as coefficients
93 : !> \param propagator AO-by-AO output matrix
94 : ! **************************************************************************************************
95 14100 : SUBROUTINE build_ao_propagator_fm(coefficients, column_factors, first_mo, work, propagator)
96 : TYPE(cp_fm_type), INTENT(IN) :: coefficients
97 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: column_factors
98 : INTEGER, INTENT(IN) :: first_mo
99 : TYPE(cp_fm_type), INTENT(INOUT) :: work, propagator
100 :
101 : INTEGER :: nao, nmo
102 :
103 14100 : IF (SIZE(column_factors) == 0) THEN
104 0 : CALL cp_fm_set_all(propagator, 0.0_dp)
105 : RETURN
106 : END IF
107 :
108 14100 : CALL cp_fm_get_info(coefficients, nrow_global=nao, ncol_global=nmo)
109 14100 : CPASSERT(first_mo + SIZE(column_factors) - 1 <= nmo)
110 :
111 14100 : CALL cp_fm_column_scale(work, column_factors, first_col=first_mo, source=coefficients)
112 :
113 : CALL parallel_gemm("N", "T", nao, nao, SIZE(column_factors), 1.0_dp, work, work, 0.0_dp, propagator, &
114 14100 : a_first_col=first_mo, b_first_col=first_mo)
115 :
116 : END SUBROUTINE build_ao_propagator_fm
117 :
118 : ! **************************************************************************************************
119 : !> \brief Builds G = X X^H, where X is the selected, column-scaled MO window.
120 : !> Only the selected workspace columns are overwritten; coefficients are unchanged.
121 : !> \param coefficients complex AO-by-MO coefficient matrix
122 : !> \param column_factors real factors for consecutive MOs starting at first_mo
123 : !> \param first_mo first global MO column in the window
124 : !> \param work coefficient workspace with the same shape and distribution as coefficients
125 : !> \param propagator complex AO-by-AO output matrix
126 : ! **************************************************************************************************
127 6912 : SUBROUTINE build_ao_propagator_cfm(coefficients, column_factors, first_mo, work, propagator)
128 : TYPE(cp_cfm_type), INTENT(IN) :: coefficients
129 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: column_factors
130 : INTEGER, INTENT(IN) :: first_mo
131 : TYPE(cp_cfm_type), INTENT(INOUT) :: work, propagator
132 :
133 : INTEGER :: nao, nmo
134 :
135 6912 : IF (SIZE(column_factors) == 0) THEN
136 0 : CALL cp_cfm_set_all(propagator, z_zero)
137 : RETURN
138 : END IF
139 :
140 6912 : CALL cp_cfm_get_info(coefficients, nrow_global=nao, ncol_global=nmo)
141 6912 : CPASSERT(first_mo + SIZE(column_factors) - 1 <= nmo)
142 :
143 : CALL cp_cfm_column_scale(work, CMPLX(column_factors, 0.0_dp, KIND=dp), &
144 51328 : first_col=first_mo, source=coefficients)
145 :
146 : CALL parallel_gemm("N", "C", nao, nao, SIZE(column_factors), z_one, work, work, z_zero, propagator, &
147 6912 : a_first_col=first_mo, b_first_col=first_mo)
148 :
149 : END SUBROUTINE build_ao_propagator_cfm
150 :
151 : ! **************************************************************************************************
152 : !> \brief ...
153 : !> \param t_3c_overl_int ...
154 : !> \param t_3c_O_compressed ...
155 : !> \param t_3c_O_ind ...
156 : !> \param t_3c_overl_int_ao_mo ...
157 : !> \param t_3c_O_mo_compressed ...
158 : !> \param t_3c_O_mo_ind ...
159 : !> \param t_3c_overl_int_gw_RI ...
160 : !> \param t_3c_overl_int_gw_AO ...
161 : !> \param starts_array_mc ...
162 : !> \param ends_array_mc ...
163 : !> \param mo_coeff ...
164 : !> \param matrix_s ...
165 : !> \param gw_corr_lev_occ ...
166 : !> \param gw_corr_lev_virt ...
167 : !> \param homo ...
168 : !> \param nmo ...
169 : !> \param para_env ...
170 : !> \param do_ic_model ...
171 : !> \param t_3c_overl_nnP_ic ...
172 : !> \param t_3c_overl_nnP_ic_reflected ...
173 : !> \param qs_env ...
174 : !> \param unit_nr ...
175 : !> \param do_alpha ...
176 : ! **************************************************************************************************
177 54 : SUBROUTINE get_tensor_3c_overl_int_gw(t_3c_overl_int, t_3c_O_compressed, t_3c_O_ind, &
178 : t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
179 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
180 108 : starts_array_mc, ends_array_mc, &
181 : mo_coeff, matrix_s, &
182 : gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
183 : para_env, &
184 : do_ic_model, &
185 : t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
186 : qs_env, unit_nr, do_alpha)
187 :
188 : TYPE(dbt_type), DIMENSION(:, :) :: t_3c_overl_int
189 : TYPE(hfx_compression_type), DIMENSION(:, :, :) :: t_3c_O_compressed
190 : TYPE(block_ind_type), DIMENSION(:, :, :) :: t_3c_O_ind
191 : TYPE(dbt_type) :: t_3c_overl_int_ao_mo
192 : TYPE(hfx_compression_type) :: t_3c_O_mo_compressed
193 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: t_3c_O_mo_ind
194 : TYPE(dbt_type) :: t_3c_overl_int_gw_RI, &
195 : t_3c_overl_int_gw_AO
196 : INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc
197 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
198 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
199 : INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo, &
200 : nmo
201 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
202 : LOGICAL, INTENT(IN) :: do_ic_model
203 : TYPE(dbt_type) :: t_3c_overl_nnP_ic, &
204 : t_3c_overl_nnP_ic_reflected
205 : TYPE(qs_environment_type), POINTER :: qs_env
206 : INTEGER, INTENT(IN) :: unit_nr
207 : LOGICAL, INTENT(IN), OPTIONAL :: do_alpha
208 :
209 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_tensor_3c_overl_int_gw'
210 :
211 : INTEGER :: cut_memory, handle, i_mem, icol_global, imo, irow_global, min_bsize, &
212 : min_bsize_mo, nkind, nmo_blk_gw, npcols, nprows, size_MO, unit_nr_prv
213 : INTEGER(int_8) :: nze
214 108 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist1, dist2, dist3, sizes_AO, &
215 54 : sizes_AO_split, sizes_MO, sizes_MO_1, &
216 108 : sizes_RI, sizes_RI_split, tmp
217 : INTEGER, DIMENSION(2) :: pdims_2d
218 : INTEGER, DIMENSION(2, 1) :: bounds
219 : INTEGER, DIMENSION(2, 3) :: ibounds
220 : INTEGER, DIMENSION(3) :: bounds_3c, pdims
221 108 : INTEGER, DIMENSION(:), POINTER :: distp_1, distp_2, sizes_MO_blocked, &
222 54 : sizes_MO_p1, sizes_MO_p2
223 : LOGICAL :: memory_info, my_do_alpha
224 : REAL(dp) :: compression_factor, memory_3c, occ
225 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: norm
226 54 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
227 : TYPE(cp_fm_type) :: fm_mat_mo_coeff_gw
228 : TYPE(dbcsr_distribution_type) :: dist, dist_templ
229 : TYPE(dbcsr_type) :: mat_mo_coeff_gw_reflected_norm, &
230 : mat_norm, mat_norm_diag, mat_work
231 : TYPE(dbcsr_type), POINTER :: mat_mo_coeff_gw, &
232 : mat_mo_coeff_gw_reflected
233 486 : TYPE(dbt_pgrid_type) :: pgrid_2d, pgrid_AO, pgrid_ic, pgrid_MO
234 702 : TYPE(dbt_type) :: mo_coeff_gw_t, mo_coeff_gw_t_tmp, &
235 378 : t_3c_overl_int_ao_ao, &
236 378 : t_3c_overl_int_mo_ao, &
237 378 : t_3c_overl_int_mo_mo
238 54 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ao, basis_set_ri_aux
239 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
240 54 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
241 :
242 54 : memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info
243 54 : IF (memory_info) THEN
244 0 : unit_nr_prv = unit_nr
245 : ELSE
246 54 : unit_nr_prv = 0
247 : END IF
248 :
249 54 : my_do_alpha = .FALSE.
250 54 : IF (PRESENT(do_alpha)) my_do_alpha = do_alpha
251 :
252 54 : CALL timeset(routineN, handle)
253 :
254 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)
255 :
256 54 : CALL cp_fm_create(fm_mat_mo_coeff_gw, mo_coeff%matrix_struct)
257 54 : CALL cp_fm_to_fm(mo_coeff, fm_mat_mo_coeff_gw)
258 :
259 : ! set MO coeffs to zero where
260 1396 : DO irow_global = 1, nmo
261 3094 : DO icol_global = 1, homo - gw_corr_lev_occ
262 3094 : CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
263 : END DO
264 32506 : DO icol_global = homo + gw_corr_lev_virt + 1, nmo
265 32452 : CALL cp_fm_set_element(fm_mat_mo_coeff_gw, irow_global, icol_global, 0.0_dp)
266 : END DO
267 : END DO
268 :
269 54 : NULLIFY (mat_mo_coeff_gw)
270 54 : CALL dbcsr_init_p(mat_mo_coeff_gw)
271 :
272 : CALL cp_dbcsr_m_by_n_from_row_template(mat_mo_coeff_gw, template=matrix_s(1)%matrix, n=nmo, &
273 54 : sym=dbcsr_type_no_symmetry)
274 :
275 : CALL copy_fm_to_dbcsr(fm_mat_mo_coeff_gw, &
276 : mat_mo_coeff_gw, &
277 54 : keep_sparsity=.FALSE.)
278 :
279 : ! just remove the blocks which have been set to zero
280 54 : CALL dbcsr_filter(mat_mo_coeff_gw, 1.0E-20_dp)
281 :
282 54 : min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize
283 54 : min_bsize_mo = qs_env%mp2_env%ri_rpa_im_time%min_bsize_mo
284 :
285 108 : CALL split_block_sizes([gw_corr_lev_occ + gw_corr_lev_virt], sizes_MO, min_bsize_mo)
286 162 : ALLOCATE (sizes_MO_1(nmo))
287 1396 : sizes_MO_1(:) = 1
288 :
289 54 : nmo_blk_gw = SIZE(sizes_MO)
290 54 : CALL move_alloc(sizes_MO, tmp)
291 162 : ALLOCATE (sizes_MO(nmo_blk_gw + 2))
292 54 : sizes_MO(1) = homo - gw_corr_lev_occ
293 108 : sizes_MO(2:SIZE(tmp) + 1) = tmp(:)
294 54 : sizes_MO(SIZE(tmp) + 2) = nmo - (homo + gw_corr_lev_virt)
295 :
296 432 : ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
297 54 : CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
298 54 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux)
299 54 : CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
300 54 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao)
301 :
302 54 : CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_AO_split)
303 54 : CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_RI_split)
304 :
305 54 : DEALLOCATE (basis_set_ao, basis_set_ri_aux)
306 :
307 54 : pdims = 0
308 : CALL dbt_pgrid_create(para_env, pdims, pgrid_AO, &
309 216 : tensor_dims=[SIZE(sizes_RI_split), SIZE(sizes_AO_split), SIZE(sizes_AO_split)])
310 :
311 54 : pdims_2d = 0
312 54 : CALL mp_dims_create(para_env%num_pe, pdims_2d)
313 :
314 : ! we iterate over MO blocks for saving memory during contraction, thus we should not parallelize over MO dimension
315 216 : pdims = [pdims_2d(1), pdims_2d(2), 1]
316 : CALL dbt_pgrid_create(para_env, pdims, pgrid_MO, &
317 216 : tensor_dims=[SIZE(sizes_RI_split), SIZE(sizes_AO_split), 1])
318 :
319 54 : pdims_2d = 0
320 : CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d, &
321 162 : tensor_dims=[SIZE(sizes_AO_split), nmo])
322 :
323 : CALL create_3c_tensor(t_3c_overl_int_ao_ao, dist1, dist2, dist3, pgrid_AO, &
324 54 : sizes_RI_split, sizes_AO_split, sizes_AO_split, [1, 2], [3], name="(RI AO | AO)")
325 54 : DEALLOCATE (dist1, dist2, dist3)
326 :
327 54 : IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
328 : CALL create_3c_tensor(t_3c_overl_int_ao_mo, dist1, dist2, dist3, pgrid_AO, &
329 36 : sizes_RI_split, sizes_AO_split, sizes_MO_1, [1, 2], [3], name="(RI AO | MO)")
330 36 : DEALLOCATE (dist1, dist2, dist3)
331 : END IF
332 :
333 : CALL create_3c_tensor(t_3c_overl_int_gw_RI, dist1, dist2, dist3, pgrid_MO, &
334 54 : sizes_RI_split, sizes_AO_split, sizes_MO, [1], [2, 3], name="(RI | AO MO)")
335 54 : DEALLOCATE (dist1, dist2, dist3)
336 :
337 : CALL create_3c_tensor(t_3c_overl_int_gw_AO, dist1, dist2, dist3, pgrid_MO, &
338 54 : sizes_AO_split, sizes_RI_split, sizes_MO, [1], [2, 3], name="(AO | RI MO)")
339 54 : DEALLOCATE (dist1, dist2, dist3)
340 :
341 54 : CALL dbt_pgrid_destroy(pgrid_AO)
342 54 : CALL dbt_pgrid_destroy(pgrid_MO)
343 :
344 54 : CALL create_2c_tensor(mo_coeff_gw_t, dist1, dist2, pgrid_2d, sizes_AO_split, sizes_MO_1, name="(AO|MO)")
345 54 : DEALLOCATE (dist1, dist2)
346 54 : CALL dbt_pgrid_destroy(pgrid_2d)
347 :
348 54 : CALL dbt_create(mat_mo_coeff_gw, mo_coeff_gw_t_tmp, name="MO coeffs")
349 54 : CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw, mo_coeff_gw_t_tmp)
350 :
351 54 : CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
352 :
353 54 : bounds(1, 1) = homo - gw_corr_lev_occ + 1
354 54 : bounds(2, 1) = homo + gw_corr_lev_virt
355 :
356 54 : CALL dbt_get_info(t_3c_overl_int_ao_ao, nfull_total=bounds_3c)
357 :
358 162 : ibounds(:, 1) = [1, bounds_3c(1)]
359 162 : ibounds(:, 3) = [1, bounds_3c(3)]
360 :
361 54 : cut_memory = SIZE(starts_array_mc)
362 :
363 54 : IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
364 72 : DO i_mem = 1, cut_memory
365 : 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), &
366 36 : qs_env%mp2_env%ri_rpa_im_time%eps_compress)
367 :
368 108 : ibounds(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
369 :
370 36 : CALL dbt_copy(t_3c_overl_int(1, 1), t_3c_overl_int_ao_ao, move_data=.TRUE.)
371 :
372 : CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_ao_ao, 1.0_dp, &
373 : t_3c_overl_int_ao_mo, contract_1=[1], notcontract_1=[2], &
374 : contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
375 72 : bounds_2=ibounds, move_data=.FALSE., unit_nr=unit_nr_prv)
376 :
377 : END DO
378 : END IF
379 :
380 54 : CALL cp_fm_release(fm_mat_mo_coeff_gw)
381 :
382 54 : IF (do_ic_model) THEN
383 2 : pdims = 0
384 : CALL dbt_pgrid_create(para_env, pdims, pgrid_ic, &
385 8 : tensor_dims=[SIZE(sizes_RI_split), nmo, nmo])
386 :
387 : CALL create_3c_tensor(t_3c_overl_int_mo_ao, dist1, dist2, dist3, pgrid_ic, &
388 2 : sizes_RI_split, sizes_MO_1, sizes_AO_split, [1, 2], [3], name="(RI MO | AO)")
389 2 : DEALLOCATE (dist1, dist2, dist3)
390 : CALL create_3c_tensor(t_3c_overl_int_mo_mo, dist1, dist2, dist3, pgrid_ic, &
391 2 : sizes_RI_split, sizes_MO_1, sizes_MO_1, [1, 2], [3], name="(RI MO | MO)")
392 2 : DEALLOCATE (dist1, dist2, dist3)
393 2 : CALL dbt_create(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic)
394 : CALL create_3c_tensor(t_3c_overl_nnP_ic_reflected, dist1, dist2, dist3, pgrid_ic, &
395 2 : sizes_RI_split, sizes_MO_1, sizes_MO_1, [1], [2, 3], name="(RI | MO MO)")
396 2 : DEALLOCATE (dist1, dist2, dist3)
397 :
398 2 : CALL dbt_pgrid_destroy(pgrid_ic)
399 :
400 2 : CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_mo_ao, order=[1, 3, 2])
401 : CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_mo_ao, 0.0_dp, &
402 : t_3c_overl_int_mo_mo, contract_1=[1], notcontract_1=[2], &
403 : contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
404 2 : bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
405 2 : CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic)
406 :
407 2 : NULLIFY (mat_mo_coeff_gw_reflected)
408 2 : CALL dbcsr_init_p(mat_mo_coeff_gw_reflected)
409 :
410 : CALL cp_dbcsr_m_by_n_from_row_template(mat_mo_coeff_gw_reflected, template=matrix_s(1)%matrix, n=nmo, &
411 2 : sym=dbcsr_type_no_symmetry)
412 :
413 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)
414 :
415 : ! normalize reflected MOs (they are not properly normalized since high angular momentum basis functions
416 : ! of the image molecule are not exactly reflected at the image plane (sign problem in p_z function)
417 2 : CALL dbcsr_create(matrix=mat_work, template=mat_mo_coeff_gw_reflected, matrix_type=dbcsr_type_no_symmetry)
418 :
419 2 : CALL dbcsr_get_info(mat_work, distribution=dist_templ, nblkcols_total=size_MO, col_blk_size=sizes_MO_blocked)
420 :
421 2 : CALL dbcsr_distribution_get(dist_templ, nprows=nprows, npcols=npcols)
422 :
423 8 : ALLOCATE (distp_1(size_MO), distp_2(size_MO))
424 2 : CALL dbt_default_distvec(size_MO, nprows, sizes_MO_blocked, distp_1)
425 2 : CALL dbt_default_distvec(size_MO, npcols, sizes_MO_blocked, distp_2)
426 2 : CALL dbcsr_distribution_new(dist, template=dist_templ, row_dist=distp_1, col_dist=distp_2, reuse_arrays=.TRUE.)
427 :
428 4 : ALLOCATE (sizes_MO_p1(size_MO))
429 4 : ALLOCATE (sizes_MO_p2(size_MO))
430 14 : sizes_MO_p1(:) = sizes_MO_blocked
431 14 : sizes_MO_p2(:) = sizes_MO_blocked
432 : CALL dbcsr_create(mat_norm, "mo norm", dist, dbcsr_type_no_symmetry, sizes_MO_p1, sizes_MO_p2, &
433 2 : reuse_arrays=.TRUE.)
434 2 : CALL dbcsr_distribution_release(dist)
435 :
436 2 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s(1)%matrix, mat_mo_coeff_gw_reflected, 0.0_dp, mat_work)
437 2 : CALL dbcsr_multiply("T", "N", 1.0_dp, mat_mo_coeff_gw_reflected, mat_work, 0.0_dp, mat_norm)
438 :
439 2 : CALL dbcsr_release(mat_work)
440 :
441 6 : ALLOCATE (norm(nmo))
442 2 : norm = 0.0_dp
443 :
444 2 : CALL dbcsr_get_diag(mat_norm, norm)
445 2 : CALL para_env%sum(norm)
446 :
447 14 : DO imo = bounds(1, 1), bounds(2, 1)
448 14 : norm(imo) = 1.0_dp/SQRT(norm(imo))
449 : END DO
450 :
451 2 : CALL dbcsr_create(mat_norm_diag, template=mat_norm)
452 2 : CALL dbcsr_release(mat_norm)
453 :
454 2 : CALL dbcsr_add_on_diag(mat_norm_diag, 1.0_dp)
455 :
456 2 : CALL dbcsr_set_diag(mat_norm_diag, norm)
457 :
458 2 : CALL dbcsr_create(mat_mo_coeff_gw_reflected_norm, template=mat_mo_coeff_gw_reflected)
459 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)
460 2 : CALL dbcsr_release(mat_norm_diag)
461 :
462 2 : CALL dbcsr_filter(mat_mo_coeff_gw_reflected_norm, 1.0E-20_dp)
463 :
464 2 : CALL dbt_copy_matrix_to_tensor(mat_mo_coeff_gw_reflected_norm, mo_coeff_gw_t_tmp)
465 2 : CALL dbcsr_release(mat_mo_coeff_gw_reflected_norm)
466 2 : CALL dbt_copy(mo_coeff_gw_t_tmp, mo_coeff_gw_t)
467 :
468 : CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_ao_ao, 0.0_dp, &
469 : t_3c_overl_int_ao_mo, contract_1=[1], notcontract_1=[2], &
470 : contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
471 2 : bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
472 :
473 2 : CALL dbt_copy(t_3c_overl_int_ao_mo, t_3c_overl_int_mo_ao, order=[1, 3, 2])
474 : CALL dbt_contract(1.0_dp, mo_coeff_gw_t, t_3c_overl_int_mo_ao, 0.0_dp, &
475 : t_3c_overl_int_mo_mo, contract_1=[1], notcontract_1=[2], &
476 : contract_2=[3], notcontract_2=[1, 2], map_1=[3], map_2=[1, 2], &
477 2 : bounds_2=bounds, move_data=.FALSE., unit_nr=unit_nr_prv)
478 2 : CALL dbt_copy(t_3c_overl_int_mo_mo, t_3c_overl_nnP_ic_reflected)
479 2 : CALL dbt_destroy(t_3c_overl_int_mo_ao)
480 2 : CALL dbt_destroy(t_3c_overl_int_mo_mo)
481 :
482 4 : CALL dbcsr_release_p(mat_mo_coeff_gw_reflected)
483 :
484 : END IF
485 :
486 54 : IF (.NOT. qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
487 36 : CALL alloc_containers(t_3c_O_mo_compressed, 1)
488 36 : CALL get_tensor_occupancy(t_3c_overl_int_ao_mo, nze, occ)
489 36 : memory_3c = 0.0_dp
490 :
491 : CALL compress_tensor(t_3c_overl_int_ao_mo, t_3c_O_mo_ind, t_3c_O_mo_compressed, &
492 36 : qs_env%mp2_env%ri_rpa_im_time%eps_compress, memory_3c)
493 :
494 36 : CALL para_env%sum(memory_3c)
495 36 : compression_factor = REAL(nze, dp)*1.0E-06_dp*8.0_dp/memory_3c
496 :
497 36 : IF (unit_nr > 0) THEN
498 : WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") &
499 18 : "MEMORY_INFO| Memory of MO-contracted tensor (compressed):", memory_3c, ' MiB'
500 :
501 : WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") &
502 18 : "MEMORY_INFO| Compression factor: ", compression_factor
503 : END IF
504 : END IF
505 :
506 54 : CALL dbcsr_release_p(mat_mo_coeff_gw)
507 :
508 54 : CALL dbt_destroy(t_3c_overl_int_ao_ao)
509 54 : CALL dbt_destroy(mo_coeff_gw_t)
510 54 : CALL dbt_destroy(mo_coeff_gw_t_tmp)
511 :
512 54 : CALL timestop(handle)
513 :
514 324 : END SUBROUTINE get_tensor_3c_overl_int_gw
515 :
516 : ! **************************************************************************************************
517 : !> \brief reflect from V = (A,B|B,A) to V_reflected = (B,A|A,B) where A belongs to the block of the molecule
518 : !> and B to the off diagonal block between molecule and image of the molecule
519 : !> \param mat_reflected ...
520 : !> \param mat_orig ...
521 : !> \param para_env ...
522 : !> \param qs_env ...
523 : !> \param unit_nr ...
524 : !> \param do_alpha ...
525 : ! **************************************************************************************************
526 2 : SUBROUTINE reflect_mat_row(mat_reflected, mat_orig, para_env, qs_env, unit_nr, do_alpha)
527 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_reflected
528 : TYPE(dbcsr_type), INTENT(IN) :: mat_orig
529 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
530 : TYPE(qs_environment_type), POINTER :: qs_env
531 : INTEGER, INTENT(IN) :: unit_nr
532 : LOGICAL, INTENT(IN) :: do_alpha
533 :
534 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reflect_mat_row'
535 :
536 : INTEGER :: handle, i_atom, j_atom, natom, &
537 : nblkcols_total, nblkrows_total
538 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: image_atom
539 2 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
540 : LOGICAL :: found_image_atom
541 : REAL(KIND=dp) :: avg_z_dist, delta, eps_dist2, &
542 : min_z_dist, ra(3), rb(3), sum_z, &
543 : z_reflection
544 : TYPE(cell_type), POINTER :: cell
545 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
546 :
547 2 : CALL timeset(routineN, handle)
548 :
549 2 : CALL dbcsr_reserve_all_blocks(mat_reflected)
550 :
551 : CALL get_qs_env(qs_env, cell=cell, &
552 2 : particle_set=particle_set)
553 :
554 : ! first check, whether we have an image molecule
555 : CALL dbcsr_get_info(mat_orig, &
556 : nblkrows_total=nblkrows_total, &
557 : nblkcols_total=nblkcols_total, &
558 : row_blk_size=row_blk_sizes, &
559 2 : col_blk_size=col_blk_sizes)
560 :
561 2 : natom = SIZE(particle_set)
562 2 : CPASSERT(natom == nblkrows_total)
563 :
564 2 : eps_dist2 = qs_env%mp2_env%ri_g0w0%eps_dist
565 2 : eps_dist2 = eps_dist2*eps_dist2
566 :
567 2 : sum_z = 0.0_dp
568 :
569 18 : DO i_atom = 1, natom
570 :
571 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
572 :
573 18 : sum_z = sum_z + ra(3)
574 :
575 : END DO
576 :
577 2 : z_reflection = sum_z/REAL(natom, KIND=dp)
578 :
579 2 : sum_z = 0.0_dp
580 :
581 18 : DO i_atom = 1, natom
582 :
583 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
584 :
585 18 : sum_z = sum_z + ABS(ra(3) - z_reflection)
586 :
587 : END DO
588 :
589 2 : avg_z_dist = sum_z/REAL(natom, KIND=dp)
590 :
591 2 : min_z_dist = avg_z_dist
592 :
593 18 : DO i_atom = 1, natom
594 :
595 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
596 :
597 18 : IF (ABS(ra(3) - z_reflection) < min_z_dist) THEN
598 : min_z_dist = ABS(ra(3) - z_reflection)
599 : END IF
600 :
601 : END DO
602 :
603 2 : IF (unit_nr > 0 .AND. do_alpha) THEN
604 1 : WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Average distance of the molecule to the image plane:', &
605 2 : avg_z_dist*0.529_dp, ' A'
606 1 : WRITE (unit_nr, '(T3,A,T70,F9.2,A2)') 'IC_MODEL| Minimum distance of the molecule to the image plane:', &
607 2 : min_z_dist*0.529_dp, ' A'
608 : END IF
609 :
610 6 : ALLOCATE (image_atom(nblkrows_total))
611 2 : image_atom = 0
612 :
613 18 : DO i_atom = 1, natom
614 :
615 16 : found_image_atom = .FALSE.
616 :
617 16 : ra(:) = pbc(particle_set(i_atom)%r, cell)
618 :
619 144 : DO j_atom = 1, natom
620 :
621 128 : rb(:) = pbc(particle_set(j_atom)%r, cell)
622 :
623 128 : delta = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) + rb(3) - 2.0_dp*z_reflection)**2
624 :
625 : ! SQRT(delta) < eps_dist
626 144 : IF (delta < eps_dist2) THEN
627 : ! this CPASSERT ensures that there is at most one image atom for each atom
628 16 : CPASSERT(.NOT. found_image_atom)
629 16 : image_atom(i_atom) = j_atom
630 16 : found_image_atom = .TRUE.
631 : ! check whether we have the same basis at the image atom
632 : ! if this is wrong, check whether you have the same basis sets for the molecule and the image
633 16 : CPASSERT(row_blk_sizes(i_atom) == row_blk_sizes(j_atom))
634 : END IF
635 :
636 : END DO
637 :
638 : ! this CPASSERT ensures that there is at least one image atom for each atom
639 18 : CPASSERT(found_image_atom)
640 :
641 : END DO
642 :
643 2 : CALL cp_dbcsr_from_subgroups(mat_orig, mat_reflected, para_env, row_mapping=image_atom, copy_only=.TRUE.)
644 :
645 2 : CALL timestop(handle)
646 :
647 6 : END SUBROUTINE reflect_mat_row
648 :
649 : ! **************************************************************************************************
650 : !> \brief ...
651 : !> \param qs_env ...
652 : !> \param atom_from_basis_index ...
653 : !> \param basis_size ...
654 : !> \param basis_type ...
655 : !> \param first_bf_from_atom ...
656 : ! **************************************************************************************************
657 14898 : SUBROUTINE get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, &
658 : basis_type, first_bf_from_atom)
659 : TYPE(qs_environment_type), POINTER :: qs_env
660 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_basis_index
661 : INTEGER :: basis_size
662 : CHARACTER(LEN=*) :: basis_type
663 : INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL :: first_bf_from_atom
664 :
665 : INTEGER :: iatom, LLL, natom, nkind
666 : INTEGER, DIMENSION(:), POINTER :: row_blk_end, row_blk_start
667 14898 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set
668 14898 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
669 14898 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
670 :
671 14898 : NULLIFY (qs_kind_set, particle_set)
672 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, nkind=nkind, &
673 14898 : particle_set=particle_set)
674 :
675 44694 : ALLOCATE (row_blk_start(natom))
676 29796 : ALLOCATE (row_blk_end(natom))
677 68512 : ALLOCATE (basis_set(nkind))
678 14898 : CALL basis_set_list_setup(basis_set, basis_type, qs_kind_set)
679 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, &
680 14898 : basis=basis_set)
681 399746 : DO LLL = 1, basis_size
682 1401550 : DO iatom = 1, natom
683 1386652 : IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN
684 384848 : atom_from_basis_index(LLL) = iatom
685 : END IF
686 : END DO
687 : END DO
688 :
689 14898 : IF (PRESENT(first_bf_from_atom)) first_bf_from_atom(1:natom) = row_blk_start(:)
690 :
691 14898 : DEALLOCATE (basis_set)
692 14898 : DEALLOCATE (row_blk_start)
693 14898 : DEALLOCATE (row_blk_end)
694 :
695 14898 : END SUBROUTINE get_atom_index_from_basis_function_index
696 :
697 : ! **************************************************************************************************
698 : !> \brief ...
699 : !> \param weight_re ...
700 : !> \param weight_im ...
701 : !> \param num_cells ...
702 : !> \param iatom ...
703 : !> \param jatom ...
704 : !> \param xkp ...
705 : !> \param wkp_W ...
706 : !> \param cell ...
707 : !> \param index_to_cell ...
708 : !> \param hmat ...
709 : !> \param particle_set ...
710 : ! **************************************************************************************************
711 409772 : SUBROUTINE compute_weight_re_im(weight_re, weight_im, &
712 : num_cells, iatom, jatom, xkp, wkp_W, &
713 : cell, index_to_cell, hmat, particle_set)
714 :
715 : REAL(KIND=dp) :: weight_re, weight_im
716 : INTEGER :: num_cells, iatom, jatom
717 : REAL(KIND=dp), DIMENSION(3) :: xkp
718 : REAL(KIND=dp) :: wkp_W
719 : TYPE(cell_type), POINTER :: cell
720 : INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
721 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
722 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
723 :
724 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_weight_re_im'
725 :
726 : INTEGER :: handle, icell, n_equidistant_cells, &
727 : xcell, ycell, zcell
728 : REAL(KIND=dp) :: abs_rab_cell, abs_rab_cell_min, arg
729 : REAL(KIND=dp), DIMENSION(3) :: cell_vector, rab_cell_i
730 :
731 409772 : CALL timeset(routineN, handle)
732 :
733 409772 : weight_re = 0.0_dp
734 409772 : weight_im = 0.0_dp
735 :
736 409772 : abs_rab_cell_min = 1.0E10_dp
737 :
738 409772 : n_equidistant_cells = 0
739 :
740 3087532 : DO icell = 1, num_cells
741 :
742 2677760 : xcell = index_to_cell(1, icell)
743 2677760 : ycell = index_to_cell(2, icell)
744 2677760 : zcell = index_to_cell(3, icell)
745 :
746 42844160 : cell_vector(1:3) = MATMUL(hmat, REAL([xcell, ycell, zcell], dp))
747 :
748 : rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - &
749 10711040 : (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
750 :
751 2677760 : abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
752 :
753 409772 : IF (abs_rab_cell < abs_rab_cell_min) THEN
754 : abs_rab_cell_min = abs_rab_cell
755 : END IF
756 :
757 : END DO
758 :
759 3087532 : DO icell = 1, num_cells
760 :
761 2677760 : xcell = index_to_cell(1, icell)
762 2677760 : ycell = index_to_cell(2, icell)
763 2677760 : zcell = index_to_cell(3, icell)
764 :
765 42844160 : cell_vector(1:3) = MATMUL(hmat, REAL([xcell, ycell, zcell], dp))
766 :
767 : rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - &
768 10711040 : (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3))
769 :
770 2677760 : abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
771 :
772 3087532 : IF (abs_rab_cell < abs_rab_cell_min + 0.1_dp) THEN
773 :
774 409932 : arg = REAL(xcell, dp)*xkp(1) + REAL(ycell, dp)*xkp(2) + REAL(zcell, dp)*xkp(3)
775 :
776 409932 : weight_re = weight_re + wkp_W*COS(twopi*arg)
777 409932 : weight_im = weight_im + wkp_W*SIN(twopi*arg)
778 :
779 409932 : n_equidistant_cells = n_equidistant_cells + 1
780 :
781 : END IF
782 :
783 : END DO
784 :
785 409772 : weight_re = weight_re/REAL(n_equidistant_cells, KIND=dp)
786 409772 : weight_im = weight_im/REAL(n_equidistant_cells, KIND=dp)
787 :
788 409772 : CALL timestop(handle)
789 :
790 409772 : END SUBROUTINE compute_weight_re_im
791 :
792 : END MODULE rpa_gw_im_time_util
|
