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