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