Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Data storage and other types for propagation via RT-BSE method.
10 : !> \author Stepan Marek (01.24)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_bse_types
14 :
15 : USE kinds, ONLY: dp
16 : USE cp_fm_types, ONLY: cp_fm_type, &
17 : cp_fm_p_type, &
18 : cp_fm_release, &
19 : cp_fm_create, &
20 : cp_fm_set_all, &
21 : cp_fm_write_formatted, &
22 : cp_fm_to_fm
23 : USE cp_cfm_types, ONLY: cp_cfm_p_type, &
24 : cp_cfm_type, &
25 : cp_cfm_set_all, &
26 : cp_cfm_create, &
27 : cp_fm_to_cfm, &
28 : cp_cfm_to_fm, &
29 : cp_cfm_to_cfm, &
30 : cp_cfm_release
31 : USE cp_fm_basic_linalg, ONLY: cp_fm_invert, &
32 : cp_fm_transpose, &
33 : cp_fm_column_scale, &
34 : cp_fm_scale_and_add
35 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale
36 : USE cp_dbcsr_api, ONLY: dbcsr_type, &
37 : dbcsr_p_type, &
38 : dbcsr_create, &
39 : dbcsr_release, &
40 : dbcsr_print, &
41 : dbcsr_get_info, &
42 : dbcsr_copy, &
43 : dbcsr_set, &
44 : dbcsr_scale
45 : USE parallel_gemm_api, ONLY: parallel_gemm
46 : USE dbt_api, ONLY: dbt_type, &
47 : dbt_pgrid_type, &
48 : dbt_pgrid_create, &
49 : dbt_pgrid_destroy, &
50 : dbt_mp_environ_pgrid, &
51 : dbt_default_distvec, &
52 : dbt_distribution_type, &
53 : dbt_distribution_new, &
54 : dbt_distribution_destroy, &
55 : dbt_create, &
56 : dbt_copy_matrix_to_tensor, &
57 : dbt_get_num_blocks, &
58 : dbt_destroy
59 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm, &
60 : copy_fm_to_dbcsr
61 : USE qs_mo_types, ONLY: mo_set_type
62 : USE basis_set_types, ONLY: gto_basis_set_p_type
63 : USE basis_set_types, ONLY: gto_basis_set_p_type
64 : USE cp_control_types, ONLY: dft_control_type
65 : USE qs_environment_types, ONLY: qs_environment_type, &
66 : get_qs_env
67 : USE force_env_types, ONLY: force_env_type
68 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
69 : USE post_scf_bandstructure_utils, ONLY: create_and_init_bs_env
70 : USE rt_propagation_types, ONLY: rt_prop_type, &
71 : get_rtp
72 : USE rt_propagation_utils, ONLY: warn_section_unused
73 : USE rt_propagator_init, ONLY: rt_initialize_rho_from_mos
74 : USE rt_propagation_methods, ONLY: s_matrices_create
75 : USE qs_moments, ONLY: build_local_moment_matrix
76 : USE moments_utils, ONLY: get_reference_point
77 : USE matrix_exp, ONLY: get_nsquare_norder
78 : USE gw_integrals, ONLY: build_3c_integral_block
79 : USE gw_large_cell_Gamma, ONLY: compute_3c_integrals
80 : USE qs_tensors, ONLY: neighbor_list_3c_destroy
81 : ! USE gw_utils, ONLY: create_and_init_bs_env_for_gw,&
82 : ! setup_AO_and_RI_basis_set,&
83 : ! get_RI_basis_and_basis_function_indices,&
84 : ! set_heuristic_parameters,&
85 : ! set_parallelization_parameters,&
86 : ! allocate_and_fill_matrices_and_arrays,&
87 : ! create_tensors
88 : USE libint_wrapper, ONLY: cp_libint_static_init
89 : USE libint_2c_3c, ONLY: libint_potential_type
90 : USE input_constants, ONLY: use_mom_ref_coac, &
91 : use_mom_ref_zero, &
92 : use_mom_ref_user, &
93 : use_mom_ref_com, &
94 : rtp_bse_ham_g0w0, &
95 : rtp_bse_ham_ks, &
96 : do_taylor, &
97 : do_bch, &
98 : do_exact
99 : USE physcon, ONLY: angstrom
100 : USE mathconstants, ONLY: z_zero
101 : USE input_section_types, ONLY: section_vals_type, &
102 : section_vals_val_get, &
103 : section_vals_get_subs_vals
104 :
105 : #include "../base/base_uses.f90"
106 :
107 : IMPLICIT NONE
108 :
109 : PRIVATE
110 :
111 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse"
112 :
113 : #:include "rt_bse_macros.fypp"
114 :
115 : PUBLIC :: rtbse_env_type, &
116 : create_rtbse_env, &
117 : release_rtbse_env, &
118 : multiply_cfm_fm, &
119 : multiply_fm_cfm, &
120 : create_hartree_ri_3c, &
121 : create_sigma_workspace_qs_only
122 :
123 : ! Created so that we can have an array of pointers to arrays
124 : TYPE series_real_type
125 : REAL(kind=dp), DIMENSION(:), POINTER :: series => NULL()
126 : END TYPE series_real_type
127 : TYPE series_complex_type
128 : COMPLEX(kind=dp), DIMENSION(:), POINTER :: series => NULL()
129 : END TYPE series_complex_type
130 :
131 : ! **************************************************************************************************
132 : !> \param n_spin Number of spin channels that are present
133 : !> \param n_ao Number of atomic orbitals
134 : !> \param n_RI Number of RI orbitals
135 : !> \param n_occ Number of occupied orbitals, spin dependent
136 : !> \param spin_degeneracy Number of electrons per orbital
137 : !> \param field Electric field calculated at the given timestep
138 : !> \param moments Moment operators along cartesian directions - centered at zero charge - used for plotting
139 : !> \param moments_field Moment operators along cartesian directions - used to coupling to the field -
140 : !> origin bound to unit cell
141 : !> \param sim_step Current step of the simulation
142 : !> \param sim_start Starting step of the simulation
143 : !> \param sim_nsteps Number of steps of the simulation
144 : !> \param sim_time Current time of the simulation
145 : !> \param sim_dt Timestep of the simulation
146 : !> \param etrs_threshold Self-consistency threshold for enforced time reversal symmetry propagation
147 : !> \param exp_accuracy Threshold for matrix exponential calculation
148 : !> \param dft_control DFT control parameters
149 : !> \param ham_effective Real and imaginary part of the effective Hamiltonian used to propagate
150 : !> the density matrix
151 : !> \param ham_reference Reference Hamiltonian, which does not change in the propagation = DFT+G0W0 - initial Hartree - initial COHSEX
152 : !> \param ham_workspace Workspace matrices for use with the Hamiltonian propagation - storage of
153 : !> exponential propagators etc.
154 : !> \param rho Density matrix at the current time step
155 : !> \param rho_new Density matrix - workspace in ETRS
156 : !> \param rho_last Density matrix - workspace in ETRS
157 : !> \param rho_new_last Density matrix - workspace in ETRS
158 : !> \param rho_M Density matrix - workspace in ETRS
159 : !> \param S_inv_fm Inverse overlap matrix, full matrix
160 : !> \param S_fm Overlap matrix, full matrix
161 : !> \param S_inv Inverse overlap matrix, sparse matrix
162 : !> \param rho_dbcsr Density matrix, sparse matrix
163 : !> \param rho_workspace Matrices for storage of density matrix at different timesteps for
164 : !> interpolation and self-consistency checks etc.
165 : !> \param complex_workspace Workspace for complex density (exact diagonalisation)
166 : !> \param complex_s Complex overlap matrix (exact diagonalisation)
167 : !> \param real_eigvals Eigenvalues of hermitian matrix (exact diagonalisation)
168 : !> \param exp_eigvals Exponentiated eigenvalues (exact diagonalisation)
169 : !> \param v_dbcsr Sparse matrix with bare Coulomb in RI basis
170 : !> \param w_dbcsr Sparse matrix with correlation part of dressed Coulomb in RI basis (without bare Coulomb)
171 : !> \param screened_dbt Tensor for screened Coulomb interaction
172 : !> \param greens_dbt Tensor for greens function/density matrix
173 : !> \param t_3c_w Tensor containing 3c integrals
174 : !> \param t_3c_work_RI_AO__AO Tensor sigma contraction
175 : !> \param t_3c_work2_RI_AO__AO Tensor sigma contraction
176 : !> \param sigma_SEX Screened exchange self-energy
177 : !> \param sigma_COH Coulomb hole self-energy
178 : !> \param hartree_curr Current Hartree matrix
179 : !> \param etrs_max_iter Maximum number of ETRS iterations
180 : !> \param ham_reference_type Which Hamiltonian to use as single particle basis
181 : !> \param mat_exp_method Which method to use for matrix exponentiation
182 : !> \param unit_nr Number of output unit
183 : !> \param int_3c_array Array containing the local 3c integrals
184 : !> \author Stepan Marek (01.24)
185 : ! **************************************************************************************************
186 : TYPE rtbse_env_type
187 : INTEGER :: n_spin = 1, &
188 : n_ao = -1, &
189 : n_RI = -1
190 : INTEGER, DIMENSION(2) :: n_occ = -1
191 : REAL(kind=dp) :: spin_degeneracy = 2
192 : REAL(kind=dp), DIMENSION(3) :: field = 0.0_dp
193 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: moments => NULL(), &
194 : moments_field => NULL()
195 : INTEGER :: sim_step = 0, &
196 : sim_start = 0, &
197 : ! Needed for continuation runs for loading of previous moments trace
198 : sim_start_orig = 0, &
199 : sim_nsteps = -1
200 : REAL(kind=dp) :: sim_time = 0.0_dp, &
201 : sim_dt = 0.1_dp, &
202 : etrs_threshold = 1.0e-7_dp, &
203 : exp_accuracy = 1.0e-10_dp, &
204 : ft_damping = 0.0_dp, &
205 : ft_start = 0.0_dp
206 : ! Which element of polarizability to print out
207 : INTEGER, DIMENSION(2) :: pol_element = [1, 1]
208 : TYPE(dft_control_type), POINTER :: dft_control => NULL()
209 : ! DEBUG : Trying keeping the reference to previous environments inside this one
210 : TYPE(qs_environment_type), POINTER :: qs_env => NULL()
211 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env => NULL()
212 : ! Stores data needed for reading/writing to the restart files
213 : TYPE(section_vals_type), POINTER :: restart_section => NULL(), &
214 : field_section => NULL(), &
215 : rho_section => NULL(), &
216 : ft_section => NULL(), &
217 : pol_section => NULL(), &
218 : moments_section => NULL()
219 : LOGICAL :: restart_extracted = .FALSE.
220 :
221 : ! Different indices signify different spins
222 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_effective => NULL()
223 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_reference => NULL()
224 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_workspace => NULL()
225 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: sigma_SEX => NULL()
226 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: sigma_COH => NULL(), &
227 : hartree_curr => NULL()
228 :
229 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho => NULL(), &
230 : rho_new => NULL(), &
231 : rho_new_last => NULL(), &
232 : rho_M => NULL(), &
233 : rho_orig => NULL()
234 : TYPE(cp_fm_type) :: S_inv_fm = cp_fm_type(), &
235 : S_fm = cp_fm_type()
236 : ! Many routines require overlap in the complex format
237 : TYPE(cp_cfm_type) :: S_cfm = cp_cfm_type()
238 : TYPE(dbcsr_type) :: rho_dbcsr = dbcsr_type(), &
239 : v_ao_dbcsr = dbcsr_type()
240 : ! Indices only correspond to different workspaces
241 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_workspace => NULL()
242 : ! Many methods use real and imaginary parts separately - prevent unnecessary reallocation
243 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: real_workspace => NULL()
244 : ! Workspace required for exact matrix exponentiation
245 : REAL(kind=dp), DIMENSION(:), POINTER :: real_eigvals => NULL()
246 : COMPLEX(kind=dp), DIMENSION(:), POINTER :: exp_eigvals => NULL()
247 : ! Workspace for saving the values for FT
248 : TYPE(series_complex_type), DIMENSION(3) :: moments_trace = series_complex_type()
249 : TYPE(series_real_type) :: time_trace = series_real_type()
250 : TYPE(series_real_type), DIMENSION(3) :: field_trace = series_real_type()
251 : ! Workspace required for hartree_pw
252 : TYPE(dbcsr_type) :: v_dbcsr = dbcsr_type(), &
253 : w_dbcsr = dbcsr_type()
254 : #if defined(FTN_NO_DEFAULT_INIT)
255 : TYPE(dbt_type) :: screened_dbt, &
256 : greens_dbt, &
257 : t_3c_w, &
258 : t_3c_work_RI_AO__AO, &
259 : t_3c_work2_RI_AO__AO
260 : #else
261 : TYPE(dbt_type) :: screened_dbt = dbt_type(), &
262 : greens_dbt = dbt_type(), &
263 : t_3c_w = dbt_type(), &
264 : t_3c_work_RI_AO__AO = dbt_type(), &
265 : t_3c_work2_RI_AO__AO = dbt_type()
266 : #endif
267 : ! These matrices are always real
268 : INTEGER :: etrs_max_iter = 10
269 : INTEGER :: ham_reference_type = 2
270 : INTEGER :: mat_exp_method = 4
271 : INTEGER :: unit_nr = -1
272 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c_array => NULL()
273 :
274 : END TYPE rtbse_env_type
275 :
276 : CONTAINS
277 :
278 : ! **************************************************************************************************
279 : !> \brief Allocates structures and prepares rtbse_env for run
280 : !> \param rtbse_env rtbse_env_type that is initialised
281 : !> \param qs_env Entry point of the calculation
282 : !> \author Stepan Marek
283 : !> \date 02.2024
284 : ! **************************************************************************************************
285 12 : SUBROUTINE create_rtbse_env(rtbse_env, qs_env, force_env)
286 : TYPE(rtbse_env_type), POINTER :: rtbse_env
287 : TYPE(qs_environment_type), POINTER :: qs_env
288 : TYPE(force_env_type), POINTER :: force_env
289 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
290 : TYPE(rt_prop_type), POINTER :: rtp
291 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
292 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
293 : INTEGER :: i, k
294 : TYPE(section_vals_type), POINTER :: input, bs_sec, md_sec
295 12 : INTEGER, DIMENSION(:), POINTER :: pol_tmp
296 :
297 : ! Allocate the storage for the gwbse environment
298 : NULLIFY (rtbse_env)
299 684 : ALLOCATE (rtbse_env)
300 : ! Extract the other types first
301 : CALL get_qs_env(qs_env, &
302 : bs_env=bs_env, &
303 : rtp=rtp, &
304 : matrix_s=matrix_s, &
305 : mos=mos, &
306 : dft_control=rtbse_env%dft_control, &
307 12 : input=input)
308 12 : bs_sec => section_vals_get_subs_vals(input, "PROPERTIES%BANDSTRUCTURE")
309 12 : IF (.NOT. ASSOCIATED(bs_env)) THEN
310 0 : CPABORT("Cannot run RT-BSE without running GW calculation (PROPERTIES) before")
311 : END IF
312 : ! Number of spins
313 12 : rtbse_env%n_spin = bs_env%n_spin
314 : ! Number of atomic orbitals
315 12 : rtbse_env%n_ao = bs_env%n_ao
316 : ! Number of auxiliary basis orbitals
317 12 : rtbse_env%n_RI = bs_env%n_RI
318 : ! Number of occupied orbitals - for closed shell equals to half the number of electrons
319 72 : rtbse_env%n_occ(:) = bs_env%n_occ(:)
320 : ! Spin degeneracy - number of spins per orbital
321 12 : rtbse_env%spin_degeneracy = bs_env%spin_degeneracy
322 : ! Default field is zero
323 48 : rtbse_env%field(:) = 0.0_dp
324 : ! Default time is zero
325 12 : rtbse_env%sim_step = 0
326 12 : rtbse_env%sim_time = 0
327 : ! Time step is taken from rtp
328 12 : md_sec => section_vals_get_subs_vals(force_env%root_section, "MOTION%MD")
329 12 : CALL section_vals_val_get(md_sec, "TIMESTEP", r_val=rtbse_env%sim_dt)
330 : ! rtbse_env%sim_dt = rtp%dt
331 : ! Threshold for etrs is taken from the eps_energy from RT propagation
332 12 : rtbse_env%etrs_threshold = rtbse_env%dft_control%rtp_control%eps_ener
333 12 : rtbse_env%exp_accuracy = rtbse_env%dft_control%rtp_control%eps_exp
334 : ! Recover custom options
335 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%RTBSE_HAMILTONIAN", &
336 12 : i_val=rtbse_env%ham_reference_type)
337 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAX_ITER", &
338 12 : i_val=rtbse_env%etrs_max_iter)
339 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAT_EXP", &
340 12 : i_val=rtbse_env%mat_exp_method)
341 : ! Output unit number, recovered from the post_scf_bandstructure_type
342 12 : rtbse_env%unit_nr = bs_env%unit_nr
343 : ! Sim start index and total number of steps as well
344 12 : CALL section_vals_val_get(md_sec, "STEP_START_VAL", i_val=rtbse_env%sim_start)
345 : ! Copy this value to sim_start_orig for continuation runs
346 12 : rtbse_env%sim_start_orig = rtbse_env%sim_start
347 12 : CALL section_vals_val_get(md_sec, "STEPS", i_val=rtbse_env%sim_nsteps)
348 : ! Get the values for the FT
349 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS_FT%DAMPING", &
350 12 : r_val=rtbse_env%ft_damping)
351 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS_FT%START_TIME", &
352 12 : r_val=rtbse_env%ft_start)
353 : ! TODO : Units for starting time and damping - so far give atomic units
354 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%PRINT%POLARIZABILITY%ELEMENT", &
355 12 : i_vals=pol_tmp)
356 12 : IF (SIZE(pol_tmp) < 2) CPABORT("Less than two elements provided for polarizability")
357 72 : rtbse_env%pol_element(:) = pol_tmp(:)
358 : ! Do basic sanity checks for pol_element
359 36 : DO k = 1, 2
360 24 : IF (rtbse_env%pol_element(k) > 3 .OR. rtbse_env%pol_element(k) < 1) &
361 12 : CPABORT("Polarisation tensor element not 1,2 or 3 in at least one index")
362 : END DO
363 :
364 : ! Get the restart section
365 12 : rtbse_env%restart_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%RESTART")
366 12 : rtbse_env%restart_extracted = .FALSE.
367 12 : rtbse_env%field_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%FIELD")
368 12 : rtbse_env%moments_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS")
369 12 : rtbse_env%rho_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%DENSITY_MATRIX")
370 12 : rtbse_env%ft_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS_FT")
371 12 : rtbse_env%pol_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%POLARIZABILITY")
372 : ! Warn the user about print sections which are not yet implemented in the RTBSE run
373 : CALL warn_section_unused(input, "DFT%REAL_TIME_PROPAGATION%PRINT%CURRENT", &
374 12 : "CURRENT print section not yet implemented for RTBSE.")
375 : CALL warn_section_unused(input, "DFT%REAL_TIME_PROPAGATION%PRINT%E_CONSTITUENTS", &
376 12 : "E_CONSTITUENTS print section not yet implemented for RTBSE.")
377 : CALL warn_section_unused(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROGRAM_RUN_INFO", &
378 12 : "PROGRAM_RUN_INFO print section not yet implemented for RTBSE.")
379 : CALL warn_section_unused(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO", &
380 12 : "PROJECTION_MO print section not yet implemented for RTBSE.")
381 : CALL warn_section_unused(input, "DFT%REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY", &
382 12 : "RESTART_HISTORY print section not yet implemented for RTBSE.")
383 : ! DEBUG : References to previous environments
384 12 : rtbse_env%qs_env => qs_env
385 12 : rtbse_env%bs_env => bs_env
386 :
387 : ! Allocate moments matrices
388 12 : NULLIFY (rtbse_env%moments)
389 48 : ALLOCATE (rtbse_env%moments(3))
390 12 : NULLIFY (rtbse_env%moments_field)
391 48 : ALLOCATE (rtbse_env%moments_field(3))
392 48 : DO k = 1, 3
393 : ! Matrices are created from overlap template
394 : ! Values are initialized in initialize_rtbse_env
395 36 : CALL cp_fm_create(rtbse_env%moments(k), bs_env%fm_s_Gamma%matrix_struct)
396 48 : CALL cp_fm_create(rtbse_env%moments_field(k), bs_env%fm_s_Gamma%matrix_struct)
397 : END DO
398 :
399 : ! Allocate space for density propagation and other operations
400 12 : NULLIFY (rtbse_env%rho_workspace)
401 60 : ALLOCATE (rtbse_env%rho_workspace(4))
402 60 : DO i = 1, SIZE(rtbse_env%rho_workspace)
403 48 : CALL cp_cfm_create(rtbse_env%rho_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
404 60 : CALL cp_cfm_set_all(rtbse_env%rho_workspace(i), CMPLX(0.0, 0.0, kind=dp))
405 : END DO
406 : ! Allocate real workspace
407 12 : NULLIFY (rtbse_env%real_workspace)
408 12 : SELECT CASE (rtbse_env%mat_exp_method)
409 : CASE (do_exact)
410 0 : ALLOCATE (rtbse_env%real_workspace(4))
411 : CASE (do_bch)
412 36 : ALLOCATE (rtbse_env%real_workspace(2))
413 : CASE DEFAULT
414 12 : CPABORT("Only exact and BCH matrix propagation implemented in RT-BSE")
415 : END SELECT
416 36 : DO i = 1, SIZE(rtbse_env%real_workspace)
417 24 : CALL cp_fm_create(rtbse_env%real_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
418 36 : CALL cp_fm_set_all(rtbse_env%real_workspace(i), 0.0_dp)
419 : END DO
420 : ! Allocate density matrix
421 12 : NULLIFY (rtbse_env%rho)
422 48 : ALLOCATE (rtbse_env%rho(rtbse_env%n_spin))
423 24 : DO i = 1, rtbse_env%n_spin
424 24 : CALL cp_cfm_create(rtbse_env%rho(i), matrix_struct=bs_env%fm_s_Gamma%matrix_struct)
425 : END DO
426 : ! Create the inverse overlap matrix, for use in density propagation
427 : ! Start by creating the actual overlap matrix
428 12 : CALL cp_fm_create(rtbse_env%S_fm, bs_env%fm_s_Gamma%matrix_struct)
429 12 : CALL cp_fm_create(rtbse_env%S_inv_fm, bs_env%fm_s_Gamma%matrix_struct)
430 12 : CALL cp_cfm_create(rtbse_env%S_cfm, bs_env%fm_s_Gamma%matrix_struct)
431 :
432 : ! Create the single particle hamiltonian
433 : ! Allocate workspace
434 12 : NULLIFY (rtbse_env%ham_workspace)
435 48 : ALLOCATE (rtbse_env%ham_workspace(rtbse_env%n_spin))
436 24 : DO i = 1, rtbse_env%n_spin
437 12 : CALL cp_cfm_create(rtbse_env%ham_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
438 24 : CALL cp_cfm_set_all(rtbse_env%ham_workspace(i), CMPLX(0.0, 0.0, kind=dp))
439 : END DO
440 : ! Now onto the Hamiltonian itself
441 12 : NULLIFY (rtbse_env%ham_reference)
442 48 : ALLOCATE (rtbse_env%ham_reference(rtbse_env%n_spin))
443 24 : DO i = 1, rtbse_env%n_spin
444 24 : CALL cp_cfm_create(rtbse_env%ham_reference(i), bs_env%fm_ks_Gamma(i)%matrix_struct)
445 : END DO
446 :
447 : ! Create the matrices and workspaces for ETRS propagation
448 12 : NULLIFY (rtbse_env%ham_effective)
449 12 : NULLIFY (rtbse_env%rho_new)
450 12 : NULLIFY (rtbse_env%rho_new_last)
451 12 : NULLIFY (rtbse_env%rho_M)
452 12 : NULLIFY (rtbse_env%rho_orig)
453 48 : ALLOCATE (rtbse_env%ham_effective(rtbse_env%n_spin))
454 48 : ALLOCATE (rtbse_env%rho_new(rtbse_env%n_spin))
455 48 : ALLOCATE (rtbse_env%rho_new_last(rtbse_env%n_spin))
456 48 : ALLOCATE (rtbse_env%rho_M(rtbse_env%n_spin))
457 48 : ALLOCATE (rtbse_env%rho_orig(rtbse_env%n_spin))
458 24 : DO i = 1, rtbse_env%n_spin
459 12 : CALL cp_cfm_create(rtbse_env%ham_effective(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
460 12 : CALL cp_cfm_set_all(rtbse_env%ham_effective(i), CMPLX(0.0, 0.0, kind=dp))
461 12 : CALL cp_cfm_create(rtbse_env%rho_new(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
462 12 : CALL cp_cfm_set_all(rtbse_env%rho_new(i), CMPLX(0.0, 0.0, kind=dp))
463 12 : CALL cp_cfm_create(rtbse_env%rho_new_last(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
464 12 : CALL cp_cfm_set_all(rtbse_env%rho_new_last(i), CMPLX(0.0, 0.0, kind=dp))
465 12 : CALL cp_cfm_create(rtbse_env%rho_M(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
466 12 : CALL cp_cfm_set_all(rtbse_env%rho_M(i), CMPLX(0.0, 0.0, kind=dp))
467 24 : CALL cp_cfm_create(rtbse_env%rho_orig(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
468 : END DO
469 :
470 : ! Fields for exact diagonalisation
471 12 : NULLIFY (rtbse_env%real_eigvals)
472 36 : ALLOCATE (rtbse_env%real_eigvals(rtbse_env%n_ao))
473 36 : rtbse_env%real_eigvals(:) = 0.0_dp
474 12 : NULLIFY (rtbse_env%exp_eigvals)
475 36 : ALLOCATE (rtbse_env%exp_eigvals(rtbse_env%n_ao))
476 36 : rtbse_env%exp_eigvals(:) = CMPLX(0.0, 0.0, kind=dp)
477 :
478 : ! Workspace for FT - includes (in principle) the zeroth step and the extra last step
479 48 : DO i = 1, 3
480 36 : NULLIFY (rtbse_env%moments_trace(i)%series)
481 2580 : ALLOCATE (rtbse_env%moments_trace(i)%series(rtbse_env%sim_nsteps + 2), source=z_zero)
482 36 : NULLIFY (rtbse_env%field_trace(i)%series)
483 2592 : ALLOCATE (rtbse_env%field_trace(i)%series(rtbse_env%sim_nsteps + 2), source=0.0_dp)
484 : END DO
485 12 : NULLIFY (rtbse_env%time_trace%series)
486 860 : ALLOCATE (rtbse_env%time_trace%series(rtbse_env%sim_nsteps + 2), source=0.0_dp)
487 :
488 : ! Allocate self-energy parts and dynamic Hartree potential
489 12 : NULLIFY (rtbse_env%hartree_curr)
490 12 : NULLIFY (rtbse_env%sigma_SEX)
491 12 : NULLIFY (rtbse_env%sigma_COH)
492 48 : ALLOCATE (rtbse_env%hartree_curr(rtbse_env%n_spin))
493 48 : ALLOCATE (rtbse_env%sigma_SEX(rtbse_env%n_spin))
494 48 : ALLOCATE (rtbse_env%sigma_COH(rtbse_env%n_spin))
495 24 : DO i = 1, rtbse_env%n_spin
496 12 : CALL cp_fm_create(rtbse_env%sigma_COH(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
497 12 : CALL cp_cfm_create(rtbse_env%sigma_SEX(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
498 12 : CALL cp_fm_create(rtbse_env%hartree_curr(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
499 12 : CALL cp_fm_set_all(rtbse_env%sigma_COH(i), 0.0_dp)
500 12 : CALL cp_cfm_set_all(rtbse_env%sigma_SEX(i), CMPLX(0.0, 0.0, kind=dp))
501 24 : CALL cp_fm_set_all(rtbse_env%hartree_curr(i), 0.0_dp)
502 : END DO
503 :
504 : ! Allocate workspaces for get_sigma
505 12 : CALL create_sigma_workspace(rtbse_env, qs_env)
506 :
507 : ! Depending on the chosen methods, allocate extra workspace
508 12 : CALL create_hartree_ri_workspace(rtbse_env, qs_env)
509 :
510 12 : END SUBROUTINE
511 :
512 : ! **************************************************************************************************
513 : !> \brief Simple reimplementation of cp_fm_release_pp1 for complex matrices
514 : !> \param matrices cp_cfm_p_type(:)
515 : !> \author Stepan Marek
516 : !> \date 02.2024
517 : ! **************************************************************************************************
518 120 : SUBROUTINE cp_cfm_release_pa1(matrices)
519 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: matrices
520 : INTEGER :: i
521 :
522 276 : DO i = 1, SIZE(matrices)
523 276 : CALL cp_cfm_release(matrices(i))
524 : END DO
525 120 : DEALLOCATE (matrices)
526 : NULLIFY (matrices)
527 120 : END SUBROUTINE cp_cfm_release_pa1
528 :
529 : ! **************************************************************************************************
530 : !> \brief Releases the environment allocated structures
531 : !> \param rtbse_env
532 : !> \author Stepan Marek
533 : !> \date 02.2024
534 : ! **************************************************************************************************
535 12 : SUBROUTINE release_rtbse_env(rtbse_env)
536 : TYPE(rtbse_env_type), POINTER :: rtbse_env
537 : INTEGER :: i
538 :
539 12 : CALL cp_cfm_release_pa1(rtbse_env%ham_effective)
540 12 : CALL cp_cfm_release_pa1(rtbse_env%ham_workspace)
541 12 : CALL cp_fm_release(rtbse_env%sigma_COH)
542 12 : CALL cp_cfm_release_pa1(rtbse_env%sigma_SEX)
543 12 : CALL cp_fm_release(rtbse_env%hartree_curr)
544 12 : CALL cp_cfm_release_pa1(rtbse_env%ham_reference)
545 12 : CALL cp_cfm_release_pa1(rtbse_env%rho)
546 12 : CALL cp_cfm_release_pa1(rtbse_env%rho_workspace)
547 12 : CALL cp_cfm_release_pa1(rtbse_env%rho_new)
548 12 : CALL cp_cfm_release_pa1(rtbse_env%rho_new_last)
549 12 : CALL cp_cfm_release_pa1(rtbse_env%rho_M)
550 12 : CALL cp_cfm_release_pa1(rtbse_env%rho_orig)
551 12 : CALL cp_fm_release(rtbse_env%real_workspace)
552 12 : CALL cp_fm_release(rtbse_env%S_inv_fm)
553 12 : CALL cp_fm_release(rtbse_env%S_fm)
554 12 : CALL cp_cfm_release(rtbse_env%S_cfm)
555 :
556 : ! DO i = 1, 3
557 : ! CALL cp_fm_release(rtbse_env%moments(i)%matrix)
558 : ! CALL cp_fm_release(rtbse_env%moments_field(i)%matrix)
559 : ! END DO
560 12 : CALL cp_fm_release(rtbse_env%moments)
561 12 : CALL cp_fm_release(rtbse_env%moments_field)
562 :
563 12 : CALL release_sigma_workspace(rtbse_env)
564 :
565 12 : CALL release_hartree_ri_workspace(rtbse_env)
566 :
567 12 : DEALLOCATE (rtbse_env%real_eigvals)
568 12 : DEALLOCATE (rtbse_env%exp_eigvals)
569 48 : DO i = 1, 3
570 36 : DEALLOCATE (rtbse_env%moments_trace(i)%series)
571 48 : DEALLOCATE (rtbse_env%field_trace(i)%series)
572 : END DO
573 12 : DEALLOCATE (rtbse_env%time_trace%series)
574 :
575 : ! Deallocate the neighbour list that is not deallocated in gw anymore
576 12 : IF (ASSOCIATED(rtbse_env%bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(rtbse_env%bs_env%nl_3c)
577 : ! Deallocate the storage for the environment itself
578 12 : DEALLOCATE (rtbse_env)
579 : ! Nullify to make sure it is not used again
580 : NULLIFY (rtbse_env)
581 :
582 12 : END SUBROUTINE
583 : ! **************************************************************************************************
584 : !> \brief Allocates the workspaces for Hartree RI method
585 : !> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors
586 : !> \param rtbse_env
587 : !> \param qs_env Quickstep environment - entry point of calculation
588 : !> \author Stepan Marek
589 : !> \date 05.2024
590 : ! **************************************************************************************************
591 12 : SUBROUTINE create_hartree_ri_workspace(rtbse_env, qs_env)
592 : TYPE(rtbse_env_type) :: rtbse_env
593 : TYPE(qs_environment_type), POINTER :: qs_env
594 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
595 :
596 12 : CALL get_qs_env(qs_env, bs_env=bs_env)
597 :
598 12 : CALL dbcsr_create(rtbse_env%rho_dbcsr, name="Sparse density", template=bs_env%mat_ao_ao%matrix)
599 12 : CALL dbcsr_create(rtbse_env%v_ao_dbcsr, name="Sparse Hartree", template=bs_env%mat_ao_ao%matrix)
600 :
601 : CALL create_hartree_ri_3c(rtbse_env%rho_dbcsr, rtbse_env%int_3c_array, rtbse_env%n_ao, rtbse_env%n_RI, &
602 : bs_env%basis_set_AO, bs_env%basis_set_RI, bs_env%i_RI_start_from_atom, &
603 12 : bs_env%ri_metric, qs_env, rtbse_env%unit_nr)
604 12 : END SUBROUTINE create_hartree_ri_workspace
605 : ! **************************************************************************************************
606 : !> \brief Separated method for allocating the 3c integrals for RI Hartree
607 : !> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors
608 : !> \param rho_dbcsr matrix used for the description of shape of 3c array
609 : !> \param int_3c 3-center integral array to be allocated and filled
610 : !> \param n_ao Number of atomic orbitals
611 : !> \param n_RI Number of auxiliary RI orbitals
612 : !> \param basis_set_AO AO basis set
613 : !> \param basis_set_RI RI auxiliary basis set
614 : !> \param i_RI_start_from_atom Array of indices where functions of a given atom in RI basis start
615 : !> \param unit_nr Unit number used for printing information about the size of int_3c
616 : !> \author Stepan Marek
617 : !> \date 01.2025
618 : ! **************************************************************************************************
619 12 : SUBROUTINE create_hartree_ri_3c(rho_dbcsr, int_3c, n_ao, n_RI, basis_set_AO, basis_set_RI, &
620 12 : i_RI_start_from_atom, ri_metric, qs_env, unit_nr)
621 : TYPE(dbcsr_type) :: rho_dbcsr
622 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c
623 : INTEGER :: n_ao, n_RI
624 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_set_AO, &
625 : basis_set_RI
626 : INTEGER, DIMENSION(:) :: i_RI_start_from_atom
627 : TYPE(libint_potential_type) :: ri_metric
628 : TYPE(qs_environment_type), POINTER :: qs_env
629 : INTEGER :: unit_nr
630 : REAL(kind=dp) :: size_mb
631 : INTEGER :: nblkrows_local, &
632 : nblkcols_local, &
633 : i_blk_local, &
634 : j_blk_local, &
635 : nrows_local, &
636 : ncols_local, &
637 : col_local_offset, &
638 : row_local_offset, &
639 : start_col_index, &
640 : end_col_index, &
641 : start_row_index, &
642 : end_row_index
643 12 : INTEGER, DIMENSION(:), POINTER :: local_blk_rows, &
644 12 : local_blk_cols, &
645 12 : row_blk_size, &
646 12 : col_blk_size
647 : ! TODO : Implement option/decision to not precompute all the 3c integrals
648 : size_mb = REAL(n_ao, kind=dp)*REAL(n_ao, kind=dp)*REAL(n_RI, kind=dp)* &
649 12 : REAL(STORAGE_SIZE(size_mb), kind=dp)/8.0_dp/1024.0_dp/1024.0_dp
650 12 : IF (unit_nr > 0) WRITE (unit_nr, '(A44,E32.2E3,A4)') &
651 6 : " RTBSE| Approximate size of the 3c integrals", size_mb, " MiB"
652 :
653 : ! Get the number of block rows and columns
654 12 : CALL dbcsr_get_info(rho_dbcsr, nblkrows_local=nblkrows_local, nblkcols_local=nblkcols_local)
655 : ! Get the global indices of local rows and columns
656 12 : CALL dbcsr_get_info(rho_dbcsr, local_rows=local_blk_rows, local_cols=local_blk_cols)
657 : ! Get the sizes of all blocks
658 12 : CALL dbcsr_get_info(rho_dbcsr, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
659 :
660 : ! Get the total required local rows and cols
661 12 : nrows_local = 0
662 24 : DO i_blk_local = 1, nblkrows_local
663 24 : nrows_local = nrows_local + row_blk_size(local_blk_rows(i_blk_local))
664 : END DO
665 12 : ncols_local = 0
666 36 : DO j_blk_local = 1, nblkcols_local
667 36 : ncols_local = ncols_local + col_blk_size(local_blk_cols(j_blk_local))
668 : END DO
669 :
670 : ! Allocate the appropriate storage
671 60 : ALLOCATE (int_3c(nrows_local, ncols_local, n_RI))
672 :
673 : ! Fill the storage with appropriate values, block by block
674 12 : row_local_offset = 1
675 24 : DO i_blk_local = 1, nblkrows_local
676 : col_local_offset = 1
677 36 : DO j_blk_local = 1, nblkcols_local
678 24 : start_row_index = row_local_offset
679 24 : end_row_index = start_row_index + row_blk_size(local_blk_rows(i_blk_local)) - 1
680 24 : start_col_index = col_local_offset
681 24 : end_col_index = start_col_index + col_blk_size(local_blk_cols(j_blk_local)) - 1
682 : CALL build_3c_integral_block(int_3c(start_row_index:end_row_index, &
683 : start_col_index:end_col_index, &
684 : 1:n_RI), &
685 : qs_env, potential_parameter=ri_metric, &
686 : basis_j=basis_set_AO, basis_k=basis_set_AO, &
687 : basis_i=basis_set_RI, &
688 : atom_j=local_blk_rows(i_blk_local), &
689 : atom_k=local_blk_cols(j_blk_local), &
690 24 : i_bf_start_from_atom=i_RI_start_from_atom)
691 36 : col_local_offset = col_local_offset + col_blk_size(local_blk_cols(j_blk_local))
692 : END DO
693 24 : row_local_offset = row_local_offset + row_blk_size(local_blk_rows(i_blk_local))
694 : END DO
695 18 : END SUBROUTINE create_hartree_ri_3c
696 : ! **************************************************************************************************
697 : !> \brief Releases the workspace for the Hartree RI method
698 : !> \param rtbse_env RT-BSE Environment, containing specific RI Hartree storage
699 : !> \author Stepan Marek
700 : !> \date 09.2024
701 : ! **************************************************************************************************
702 12 : SUBROUTINE release_hartree_ri_workspace(rtbse_env)
703 : TYPE(rtbse_env_type) :: rtbse_env
704 :
705 12 : DEALLOCATE (rtbse_env%int_3c_array)
706 :
707 12 : CALL dbcsr_release(rtbse_env%rho_dbcsr)
708 :
709 12 : CALL dbcsr_release(rtbse_env%v_dbcsr)
710 :
711 12 : CALL dbcsr_release(rtbse_env%v_ao_dbcsr)
712 :
713 12 : END SUBROUTINE release_hartree_ri_workspace
714 : ! **************************************************************************************************
715 : !> \brief Allocates the workspaces for self-energy determination routine
716 : !> \param rtbse_env Structure for holding information and workspace structures
717 : !> \param qs_env Quickstep environment - entry point of calculation
718 : !> \author Stepan Marek
719 : !> \date 02.2024
720 : ! **************************************************************************************************
721 12 : SUBROUTINE create_sigma_workspace(rtbse_env, qs_env)
722 : TYPE(rtbse_env_type) :: rtbse_env
723 : TYPE(qs_environment_type), POINTER :: qs_env
724 :
725 : CALL create_sigma_workspace_qs_only(qs_env, rtbse_env%screened_dbt, rtbse_env%w_dbcsr, &
726 : rtbse_env%t_3c_w, rtbse_env%t_3c_work_RI_AO__AO, &
727 12 : rtbse_env%t_3c_work2_RI_AO__AO, rtbse_env%greens_dbt)
728 12 : END SUBROUTINE create_sigma_workspace
729 : ! **************************************************************************************************
730 : !> \brief Allocates the workspaces for self-energy determination routine
731 : !> \note Does so without referencing the rtbse_env
732 : !> \note References bs_env
733 : !> \param rtbse_env Structure for holding information and workspace structures
734 : !> \param qs_env Quickstep environment - entry point of calculation
735 : !> \author Stepan Marek
736 : !> \date 02.2024
737 : ! **************************************************************************************************
738 12 : SUBROUTINE create_sigma_workspace_qs_only(qs_env, screened_dbt, screened_dbcsr, int_3c_dbt, &
739 : work_dbt_3c_1, work_dbt_3c_2, work_dbt_2c)
740 : TYPE(qs_environment_type), POINTER :: qs_env
741 : TYPE(dbcsr_type) :: screened_dbcsr
742 : TYPE(dbt_type) :: screened_dbt, &
743 : int_3c_dbt, &
744 : work_dbt_3c_1, &
745 : work_dbt_3c_2, &
746 : work_dbt_2c
747 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
748 :
749 12 : CALL get_qs_env(qs_env, bs_env=bs_env)
750 :
751 : ! t_3c_w
752 12 : CALL dbt_create(bs_env%t_RI__AO_AO, int_3c_dbt)
753 : ! TODO : Provide option/decision whether to store the 3c integrals precomputed
754 12 : CALL compute_3c_integrals(qs_env, bs_env, int_3c_dbt)
755 : ! t_3c_work_RI_AO__AO
756 12 : CALL dbt_create(bs_env%t_RI_AO__AO, work_dbt_3c_1)
757 : ! t_3c_work2_RI_AO__AO
758 12 : CALL dbt_create(bs_env%t_RI_AO__AO, work_dbt_3c_2)
759 : ! t_W
760 : ! Populate screened_dbt from gw run
761 12 : CALL dbcsr_create(screened_dbcsr, name="W", template=bs_env%mat_RI_RI%matrix)
762 12 : CALL dbt_create(screened_dbcsr, screened_dbt)
763 : ! greens_dbt
764 12 : CALL dbt_create(bs_env%mat_ao_ao%matrix, work_dbt_2c)
765 12 : END SUBROUTINE
766 : ! **************************************************************************************************
767 : !> \brief Releases the workspaces for self-energy determination
768 : !> \param rtbse_env
769 : !> \author Stepan Marek
770 : !> \date 02.2024
771 : ! **************************************************************************************************
772 12 : SUBROUTINE release_sigma_workspace(rtbse_env)
773 : TYPE(rtbse_env_type) :: rtbse_env
774 :
775 12 : CALL dbt_destroy(rtbse_env%t_3c_w)
776 12 : CALL dbt_destroy(rtbse_env%t_3c_work_RI_AO__AO)
777 12 : CALL dbt_destroy(rtbse_env%t_3c_work2_RI_AO__AO)
778 12 : CALL dbt_destroy(rtbse_env%screened_dbt)
779 12 : CALL dbt_destroy(rtbse_env%greens_dbt)
780 12 : CALL dbcsr_release(rtbse_env%w_dbcsr)
781 12 : END SUBROUTINE
782 : ! **************************************************************************************************
783 : !> \brief Multiplies real matrix by a complex matrix from the right
784 : !> \note So far only converts the real matrix to complex one, potentially doubling the work
785 : !> \param rtbse_env
786 : !> \author Stepan Marek
787 : !> \date 09.2024
788 : ! **************************************************************************************************
789 12960 : SUBROUTINE multiply_fm_cfm(trans_r, trans_c, na, nb, nc, &
790 : alpha, matrix_r, matrix_c, beta, res)
791 : ! Transposition
792 : CHARACTER(len=1) :: trans_r, trans_c
793 : INTEGER :: na, nb, nc
794 : ! accept real numbers
795 : ! TODO : Just use complex numbers and import z_one, z_zero etc.
796 : REAL(kind=dp) :: alpha, beta
797 : TYPE(cp_fm_type) :: matrix_r
798 : TYPE(cp_cfm_type) :: matrix_c, res
799 : TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im
800 : REAL(kind=dp) :: i_unit
801 : CHARACTER(len=1) :: trans_cr
802 :
803 3240 : CALL cp_fm_create(work_re, matrix_c%matrix_struct)
804 3240 : CALL cp_fm_create(work_im, matrix_c%matrix_struct)
805 3240 : CALL cp_fm_create(res_re, res%matrix_struct)
806 3240 : CALL cp_fm_create(res_im, res%matrix_struct)
807 3240 : CALL cp_cfm_to_fm(matrix_c, work_re, work_im)
808 0 : SELECT CASE (trans_c)
809 : CASE ("C")
810 0 : i_unit = -1.0_dp
811 0 : trans_cr = "T"
812 : CASE ("T")
813 0 : i_unit = 1.0_dp
814 0 : trans_cr = "T"
815 : CASE default
816 3240 : i_unit = 1.0_dp
817 3240 : trans_cr = "N"
818 : END SELECT
819 : ! Actual multiplication
820 : CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, &
821 3240 : alpha, matrix_r, work_re, beta, res_re)
822 : CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, &
823 3240 : i_unit*alpha, matrix_r, work_im, beta, res_im)
824 3240 : CALL cp_fm_to_cfm(res_re, res_im, res)
825 3240 : CALL cp_fm_release(work_re)
826 3240 : CALL cp_fm_release(work_im)
827 3240 : CALL cp_fm_release(res_re)
828 3240 : CALL cp_fm_release(res_im)
829 :
830 3240 : END SUBROUTINE multiply_fm_cfm
831 : ! **************************************************************************************************
832 : !> \brief Multiplies complex matrix by a real matrix from the right
833 : !> \note So far only converts the real matrix to complex one, potentially doubling the work
834 : !> \param rtbse_env
835 : !> \author Stepan Marek
836 : !> \date 09.2024
837 : ! **************************************************************************************************
838 4864 : SUBROUTINE multiply_cfm_fm(trans_c, trans_r, na, nb, nc, &
839 : alpha, matrix_c, matrix_r, beta, res)
840 : ! Transposition
841 : CHARACTER(len=1) :: trans_c, trans_r
842 : INTEGER :: na, nb, nc
843 : ! accept real numbers
844 : ! TODO : complex number support via interface?
845 : REAL(kind=dp) :: alpha, beta
846 : TYPE(cp_cfm_type) :: matrix_c, res
847 : TYPE(cp_fm_type) :: matrix_r
848 : TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im
849 : REAL(kind=dp) :: i_unit
850 : CHARACTER(len=1) :: trans_cr
851 :
852 1216 : CALL cp_fm_create(work_re, matrix_c%matrix_struct)
853 1216 : CALL cp_fm_create(work_im, matrix_c%matrix_struct)
854 1216 : CALL cp_fm_create(res_re, res%matrix_struct)
855 1216 : CALL cp_fm_create(res_im, res%matrix_struct)
856 1216 : CALL cp_cfm_to_fm(matrix_c, work_re, work_im)
857 0 : SELECT CASE (trans_c)
858 : CASE ("C")
859 0 : i_unit = -1.0_dp
860 0 : trans_cr = "T"
861 : CASE ("T")
862 0 : i_unit = 1.0_dp
863 0 : trans_cr = "T"
864 : CASE default
865 1216 : i_unit = 1.0_dp
866 1216 : trans_cr = "N"
867 : END SELECT
868 : ! Actual multiplication
869 : CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, &
870 1216 : alpha, work_re, matrix_r, beta, res_re)
871 : CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, &
872 1216 : i_unit*alpha, work_im, matrix_r, beta, res_im)
873 1216 : CALL cp_fm_to_cfm(res_re, res_im, res)
874 1216 : CALL cp_fm_release(work_re)
875 1216 : CALL cp_fm_release(work_im)
876 1216 : CALL cp_fm_release(res_re)
877 1216 : CALL cp_fm_release(res_im)
878 :
879 1216 : END SUBROUTINE multiply_cfm_fm
880 0 : END MODULE
|
