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