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