Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the propagation via RT-BSE method.
10 : !> \note The control is handed directly from cp2k_runs
11 : !> \author Stepan Marek (12.23)
12 : ! **************************************************************************************************
13 :
14 : MODULE rt_bse
15 : USE bibliography, ONLY: Marek2025, &
16 : cite_reference
17 : USE qs_environment_types, ONLY: get_qs_env
18 : USE force_env_types, ONLY: force_env_type
19 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
20 : USE cp_fm_types, ONLY: cp_fm_type, &
21 : cp_fm_to_fm, &
22 : cp_fm_create, &
23 : cp_fm_set_all, &
24 : cp_fm_release
25 : USE cp_cfm_types, ONLY: cp_cfm_type, &
26 : cp_fm_to_cfm, &
27 : cp_cfm_to_cfm, &
28 : cp_cfm_to_fm, &
29 : cp_cfm_create, &
30 : cp_cfm_get_info, &
31 : cp_cfm_release
32 : USE kinds, ONLY: dp
33 : USE cp_dbcsr_api, ONLY: dbcsr_p_type, &
34 : dbcsr_type, &
35 : dbcsr_create, &
36 : dbcsr_release, &
37 : dbcsr_copy, &
38 : dbcsr_add, &
39 : dbcsr_set, &
40 : dbcsr_clear, &
41 : dbcsr_iterator_type, &
42 : dbcsr_iterator_start, &
43 : dbcsr_iterator_stop, &
44 : dbcsr_iterator_next_block, &
45 : dbcsr_reserve_blocks, &
46 : dbcsr_get_num_blocks, &
47 : dbcsr_get_block_p, &
48 : dbcsr_get_info
49 : USE dbt_api, ONLY: dbt_clear, &
50 : dbt_contract, &
51 : dbt_copy_matrix_to_tensor, &
52 : dbt_copy_tensor_to_matrix, &
53 : dbt_type
54 : USE libint_2c_3c, ONLY: libint_potential_type
55 : USE qs_tensors, ONLY: build_2c_integrals, &
56 : build_2c_neighbor_lists
57 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type, &
58 : release_neighbor_list_sets
59 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm, &
60 : copy_fm_to_dbcsr, &
61 : cp_dbcsr_sm_fm_multiply
62 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale, &
63 : cp_fm_invert, &
64 : cp_fm_transpose, &
65 : cp_fm_column_scale, &
66 : cp_fm_scale_and_add
67 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add, &
68 : cp_cfm_scale, &
69 : cp_cfm_transpose, &
70 : cp_cfm_norm, &
71 : cp_cfm_trace, &
72 : cp_cfm_column_scale
73 : USE cp_cfm_diag, ONLY: cp_cfm_geeig
74 : USE parallel_gemm_api, ONLY: parallel_gemm
75 : USE qs_moments, ONLY: build_local_moment_matrix
76 : USE moments_utils, ONLY: get_reference_point
77 : USE force_env_methods, ONLY: force_env_calc_energy_force
78 : USE efield_utils, ONLY: make_field
79 : USE message_passing, ONLY: mp_para_env_type
80 : USE input_constants, ONLY: rtp_bse_ham_g0w0, &
81 : use_mom_ref_zero, &
82 : do_bch, &
83 : do_exact, &
84 : use_rt_restart
85 : USE rt_bse_types, ONLY: rtbse_env_type, &
86 : create_rtbse_env, &
87 : release_rtbse_env, &
88 : multiply_fm_cfm
89 : USE rt_bse_ri_rs, ONLY: compute_sigma_ri_rs_complex
90 : USE rt_bse_io, ONLY: output_moments, &
91 : output_field, &
92 : output_mos_contravariant, &
93 : read_field, &
94 : read_restart, &
95 : output_restart, &
96 : print_timestep_info, &
97 : print_etrs_info_header, &
98 : print_etrs_info, &
99 : print_rtbse_header_info
100 : USE cp_log_handling, ONLY: cp_logger_type, &
101 : cp_get_default_logger
102 : USE cp_output_handling, ONLY: cp_add_iter_level, &
103 : cp_rm_iter_level, &
104 : cp_iterate
105 : USE rt_propagation_output, ONLY: print_ft
106 : USE rt_propagation_utils, ONLY: read_moments
107 :
108 : #include "../base/base_uses.f90"
109 :
110 : IMPLICIT NONE
111 :
112 : PRIVATE
113 :
114 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse"
115 :
116 : #:include "rt_bse_macros.fypp"
117 :
118 : PUBLIC :: run_propagation_bse, &
119 : get_hartree, &
120 : get_sigma, &
121 : initialize_rtbse_env, &
122 : initialize_singleparticle_hamiltonian, &
123 : initialize_hartree_potential, &
124 : initialize_cohsex_selfenergy, &
125 : get_idempotence_deviation, &
126 : antiherm_metric, &
127 : init_hartree, &
128 : rho_metric, &
129 : propagate_density, &
130 : cp_cfm_gexp, &
131 : get_electron_number
132 :
133 : INTERFACE get_sigma
134 : MODULE PROCEDURE get_sigma_complex, &
135 : get_sigma_real, &
136 : get_sigma_dbcsr, &
137 : get_sigma_noenv
138 : END INTERFACE
139 : INTERFACE get_hartree
140 : MODULE PROCEDURE get_hartree_env, &
141 : get_hartree_noenv
142 : END INTERFACE
143 :
144 : CONTAINS
145 :
146 : ! **************************************************************************************************
147 : !> \brief Runs the electron-only real time BSE propagation
148 : !> \param force_env Force environment data, entry point of the calculation
149 : ! **************************************************************************************************
150 14 : SUBROUTINE run_propagation_bse(force_env)
151 : TYPE(force_env_type), POINTER :: force_env
152 : CHARACTER(len=*), PARAMETER :: routineN = 'run_propagation_bse'
153 : TYPE(rtbse_env_type), POINTER :: rtbse_env
154 : INTEGER :: i, j, k, handle
155 : LOGICAL :: converged
156 : REAL(kind=dp) :: metric, enum_re, enum_im, &
157 : idempotence_dev, a_metric_1, a_metric_2
158 : TYPE(cp_logger_type), POINTER :: logger
159 :
160 14 : CALL timeset(routineN, handle)
161 :
162 : ! Bibliography information
163 14 : CALL cite_reference(Marek2025)
164 :
165 14 : logger => cp_get_default_logger()
166 :
167 : ! Run the initial SCF calculation / read SCF restart information
168 14 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., consistent_energies=.FALSE.)
169 :
170 : ! Allocate all persistant storage and read input that does not need further processing
171 14 : CALL create_rtbse_env(rtbse_env, force_env)
172 :
173 14 : CALL print_rtbse_header_info(rtbse_env)
174 :
175 : ! Initiate iteration level "MD" in order to copy the structure of other RTP codes
176 14 : CALL cp_add_iter_level(logger%iter_info, "MD")
177 : ! Initialize non-trivial values
178 : ! - calculates the moment operators
179 14 : CALL initialize_moments(rtbse_env)
180 : ! - populates overlap and inverse overlap matrices
181 14 : CALL initialize_rtbse_env(rtbse_env)
182 :
183 : ! - populates the initial density matrix
184 : ! - reads the restart density if requested
185 14 : CALL initialize_density_matrix(rtbse_env)
186 : ! - reads the moment and field traces from previous runs (no-op if the files are absent)
187 : CALL read_moments(rtbse_env%moments_section, rtbse_env%sim_start_orig, &
188 14 : rtbse_env%sim_start, rtbse_env%moments_trace, rtbse_env%time_trace)
189 14 : CALL read_field(rtbse_env)
190 : ! - calculates/populates the G0W0/KS Hamiltonian, respectively
191 14 : CALL initialize_singleparticle_hamiltonian(rtbse_env)
192 : ! - calculates the Hartree reference potential
193 14 : CALL initialize_hartree_potential(rtbse_env)
194 : ! - calculates the COHSEX reference self-energy
195 14 : CALL initialize_cohsex_selfenergy(rtbse_env)
196 :
197 : ! Setup the time based on the starting step
198 : ! Assumes identical dt between two runs
199 14 : rtbse_env%sim_time = REAL(rtbse_env%sim_start, dp)*rtbse_env%sim_dt
200 : ! Output 0 time moments and field
201 14 : IF (.NOT. rtbse_env%restart_extracted) THEN
202 10 : CALL output_field(rtbse_env, append_opt=.FALSE.)
203 10 : CALL output_moments(rtbse_env, rtbse_env%rho)
204 : END IF
205 :
206 : ! Do not apply the delta kick if we are doing a restart calculation
207 14 : IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. (.NOT. rtbse_env%restart_extracted)) THEN
208 6 : CALL apply_delta_pulse(rtbse_env)
209 : END IF
210 :
211 : ! ********************** Start the time loop **********************
212 : ! NOTE : Time-loop starts at index sim_start = 0, unless restarted or configured otherwise
213 692 : DO i = rtbse_env%sim_start, rtbse_env%sim_nsteps - 1
214 :
215 : ! Update the simulation time
216 678 : rtbse_env%sim_time = REAL(i, dp)*rtbse_env%sim_dt
217 678 : rtbse_env%sim_step = i
218 : ! Carry out the ETRS self-consistent propagation - propagates rho to rho_new (through rho_M)
219 678 : CALL etrs_scf_loop(rtbse_env, rtbse_env%rho, rtbse_env%rho_M, rtbse_env%rho_new, converged, k, metric)
220 678 : CALL get_electron_number(rtbse_env, rtbse_env%rho_new, enum_re, enum_im)
221 : IF (.FALSE.) THEN
222 : ! Not all of these are used, but they are all good metrics to check the convergence in problematic cases
223 : ! TODO : Allow for conditional warning
224 : CALL get_idempotence_deviation(rtbse_env, rtbse_env%rho_new, idempotence_dev)
225 : DO j = 1, rtbse_env%n_spin
226 : CALL cp_cfm_to_fm(rtbse_env%sigma_SEX(j), rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
227 : CALL antiherm_metric(real_fm=rtbse_env%real_workspace(1), imag_fm=rtbse_env%real_workspace(2), &
228 : workspace=rtbse_env%rho_workspace, metric=a_metric_1)
229 : CALL antiherm_metric(real_fm=rtbse_env%hartree_curr(j), &
230 : workspace=rtbse_env%rho_workspace, metric=a_metric_2)
231 : END DO
232 : END IF
233 1356 : CALL print_timestep_info(rtbse_env, i, [enum_re], metric, k)
234 678 : IF (.NOT. converged) CPABORT("ETRS did not converge")
235 678 : CALL cp_iterate(logger%iter_info, iter_nr=i, last=(i == rtbse_env%sim_nsteps))
236 1356 : DO j = 1, rtbse_env%n_spin
237 1356 : CALL cp_cfm_to_cfm(rtbse_env%rho_new(j), rtbse_env%rho(j))
238 : END DO
239 : ! Print the updated field
240 678 : CALL output_field(rtbse_env)
241 : ! If needed, print out the density matrix in MO basis
242 678 : CALL output_mos_contravariant(rtbse_env, rtbse_env%rho, rtbse_env%rho_section)
243 : ! Also handles outputting to memory
244 678 : CALL output_moments(rtbse_env, rtbse_env%rho)
245 : ! Output restart files, so that the restart starts at the following time index
246 1370 : CALL output_restart(rtbse_env, rtbse_env%rho, i + 1)
247 : END DO
248 : ! ********************** End the time loop **********************
249 :
250 14 : CALL cp_rm_iter_level(logger%iter_info, "MD")
251 :
252 : ! Carry out the FT
253 : CALL print_ft(rtbse_env%rtp_section, &
254 : rtbse_env%moments_trace, &
255 : rtbse_env%time_trace, &
256 : rtbse_env%field_trace, &
257 : rtbse_env%dft_control%rtp_control, &
258 14 : info_opt=rtbse_env%unit_nr)
259 :
260 : ! Deallocate everything
261 14 : CALL release_rtbse_env(rtbse_env)
262 :
263 14 : CALL timestop(handle)
264 14 : END SUBROUTINE run_propagation_bse
265 :
266 : ! **************************************************************************************************
267 : !> \brief Calculates the initial values, based on restart/scf density, and other non-trivial values
268 : !> \param rtbse_env RT-BSE environment
269 : !> \author Stepan Marek (09.24)
270 : ! **************************************************************************************************
271 66 : SUBROUTINE initialize_rtbse_env(rtbse_env)
272 : TYPE(rtbse_env_type), POINTER :: rtbse_env
273 : CHARACTER(len=*), PARAMETER :: routineN = "initialize_rtbse_env"
274 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
275 66 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
276 : INTEGER :: handle
277 :
278 66 : CALL timeset(routineN, handle)
279 :
280 : ! Get pointers to parameters from qs_env
281 66 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, matrix_s=matrix_s)
282 :
283 : ! ****** START OVERLAP + INVERSE OVERLAP CALCULATION
284 66 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, rtbse_env%S_fm)
285 66 : CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%S_cfm)
286 66 : CALL cp_fm_invert(rtbse_env%S_fm, rtbse_env%S_inv_fm)
287 : ! ****** END OVERLAP + INVERSE OVERLAP CALCULATION
288 :
289 66 : CALL timestop(handle)
290 66 : END SUBROUTINE initialize_rtbse_env
291 :
292 : ! **************************************************************************************************
293 : !> \brief Calculates the moment operators
294 : !> \param rtbse_env RT-BSE environment
295 : !> \author Stepan Marek (09.24)
296 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
297 : ! **************************************************************************************************
298 14 : SUBROUTINE initialize_moments(rtbse_env)
299 : TYPE(rtbse_env_type), POINTER :: rtbse_env
300 : CHARACTER(len=*), PARAMETER :: routineN = "initialize_moments"
301 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
302 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments_dbcsr_p
303 : INTEGER :: i, k, handle
304 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
305 : REAL(kind=dp), DIMENSION(3) :: rpoint
306 :
307 14 : CALL timeset(routineN, handle)
308 : ! Get pointers to parameters from qs_env
309 14 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, matrix_s=matrix_s)
310 :
311 : ! ****** START MOMENTS OPERATOR CALCULATION
312 : ! Construct moments from dbcsr
313 : NULLIFY (moments_dbcsr_p)
314 56 : ALLOCATE (moments_dbcsr_p(3))
315 56 : DO k = 1, 3
316 : ! Make sure the pointer is empty
317 42 : NULLIFY (moments_dbcsr_p(k)%matrix)
318 : ! Allocate a new matrix that the pointer points to
319 42 : ALLOCATE (moments_dbcsr_p(k)%matrix)
320 : ! Create the matrix storage - matrix copies the structure of overlap matrix
321 56 : CALL dbcsr_copy(moments_dbcsr_p(k)%matrix, matrix_s(1)%matrix)
322 : END DO
323 : ! Run the moment calculation
324 : ! check for presence to prevent memory errors
325 14 : rpoint(:) = 0.0_dp
326 : CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, &
327 14 : reference=rtbse_env%moment_ref_type, ref_point=rtbse_env%user_moment_ref_point)
328 14 : CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint)
329 : ! Copy to full matrix
330 28 : DO i = 1, rtbse_env%n_spin
331 70 : DO k = 1, 3
332 : ! AO dipole is spin-independent; replicate into each spin slot
333 56 : CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, rtbse_env%moments(k, i))
334 : END DO
335 : END DO
336 : ! Now, repeat without reference point to get the moments for field
337 : CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, &
338 14 : reference=use_mom_ref_zero)
339 14 : CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint)
340 28 : DO i = 1, rtbse_env%n_spin
341 70 : DO k = 1, 3
342 56 : CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, rtbse_env%moments_field(k, i))
343 : END DO
344 : END DO
345 :
346 : ! Now can deallocate dbcsr matrices
347 56 : DO k = 1, 3
348 42 : CALL dbcsr_release(moments_dbcsr_p(k)%matrix)
349 56 : DEALLOCATE (moments_dbcsr_p(k)%matrix)
350 : END DO
351 14 : DEALLOCATE (moments_dbcsr_p)
352 : ! ****** END MOMENTS OPERATOR CALCULATION
353 :
354 14 : CALL timestop(handle)
355 14 : END SUBROUTINE initialize_moments
356 :
357 : ! **************************************************************************************************
358 : !> \brief Calculates the initial density matrix, based on the SCF density or restart density if requested
359 : !> \param rtbse_env RT-BSE environment
360 : !> \author Stepan Marek (09.24)
361 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
362 : ! **************************************************************************************************
363 14 : SUBROUTINE initialize_density_matrix(rtbse_env)
364 : TYPE(rtbse_env_type), POINTER :: rtbse_env
365 : CHARACTER(len=*), PARAMETER :: routineN = "initialize_density_matrix"
366 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
367 : REAL(kind=dp), DIMENSION(:), POINTER :: occupations
368 : INTEGER :: i, handle
369 :
370 14 : CALL timeset(routineN, handle)
371 : ! Get pointers to parameters from qs_env
372 14 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
373 :
374 : ! ****** START INITIAL DENSITY MATRIX CALCULATION
375 : ! Get the rho from fm_MOS
376 : ! Uses real orbitals only - no kpoints
377 42 : ALLOCATE (occupations(rtbse_env%n_ao))
378 : ! Iterate over both spins
379 28 : DO i = 1, rtbse_env%n_spin
380 42 : occupations(:) = 0.0_dp
381 28 : occupations(1:rtbse_env%n_occ(i)) = 1.0_dp
382 : ! Create real part
383 14 : CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(i), rtbse_env%real_workspace(1))
384 14 : CALL cp_fm_column_scale(rtbse_env%real_workspace(1), occupations)
385 : CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
386 : 1.0_dp, rtbse_env%real_workspace(1), bs_env%fm_mo_coeff_Gamma(i), &
387 14 : 0.0_dp, rtbse_env%real_workspace(2))
388 : ! Sets imaginary part to zero
389 14 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%rho(i))
390 : ! Save the reference value for the case of delta kick
391 28 : CALL cp_cfm_to_cfm(rtbse_env%rho(i), rtbse_env%rho_orig(i))
392 : END DO
393 14 : DEALLOCATE (occupations)
394 : ! If the restart field is provided, overwrite rho from restart
395 14 : IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
396 6 : CALL read_restart(rtbse_env)
397 : END IF
398 : ! ****** END INITIAL DENSITY MATRIX CALCULATION
399 :
400 14 : CALL timestop(handle)
401 14 : END SUBROUTINE initialize_density_matrix
402 :
403 : ! **************************************************************************************************
404 : !> \brief Calculates the single particle Hamiltonian
405 : !> \param rtbse_env RT-BSE environment
406 : !> \author Stepan Marek (09.24)
407 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
408 : ! **************************************************************************************************
409 14 : SUBROUTINE initialize_singleparticle_hamiltonian(rtbse_env)
410 : TYPE(rtbse_env_type), POINTER :: rtbse_env
411 : CHARACTER(len=*), PARAMETER :: routineN = "initialize_singleparticle_hamiltonian"
412 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
413 : INTEGER :: i, handle
414 :
415 14 : CALL timeset(routineN, handle)
416 : ! Get pointers to parameters from qs_env
417 14 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
418 :
419 : ! ****** START SINGLE PARTICLE HAMILTONIAN CALCULATION
420 28 : DO i = 1, rtbse_env%n_spin
421 28 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
422 : ! G0W0 Hamiltonian
423 14 : CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(i), rtbse_env%real_workspace(1))
424 : ! NOTE : Gamma point is not always the zero k-point
425 : ! C * Lambda
426 14 : CALL cp_fm_column_scale(rtbse_env%real_workspace(1), bs_env%eigenval_G0W0(:, 1, i))
427 : ! C * Lambda * C^T
428 : CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
429 : 1.0_dp, rtbse_env%real_workspace(1), bs_env%fm_mo_coeff_Gamma(i), &
430 14 : 0.0_dp, rtbse_env%real_workspace(2))
431 : ! S * C * Lambda * C^T
432 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
433 : 1.0_dp, rtbse_env%S_fm, rtbse_env%real_workspace(2), &
434 14 : 0.0_dp, rtbse_env%real_workspace(1))
435 : ! S * C * Lambda * C^T * S = H
436 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
437 : 1.0_dp, rtbse_env%real_workspace(1), rtbse_env%S_fm, &
438 14 : 0.0_dp, rtbse_env%real_workspace(2))
439 14 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_reference(i))
440 : ELSE
441 : ! KS Hamiltonian
442 0 : CALL cp_fm_to_cfm(msourcer=bs_env%fm_ks_Gamma(i), mtarget=rtbse_env%ham_reference(i))
443 : END IF
444 : END DO
445 : ! ****** END SINGLE PARTICLE HAMILTONIAN CALCULATION
446 :
447 14 : CALL timestop(handle)
448 14 : END SUBROUTINE initialize_singleparticle_hamiltonian
449 :
450 : ! **************************************************************************************************
451 : !> \brief Calculates the Hartree potential
452 : !> \param rtbse_env RT-BSE environment
453 : !> \author Stepan Marek (09.24)
454 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
455 : ! **************************************************************************************************
456 14 : SUBROUTINE initialize_hartree_potential(rtbse_env)
457 : TYPE(rtbse_env_type), POINTER :: rtbse_env
458 : CHARACTER(len=*), PARAMETER :: routineN = "initialize_hartree_potential"
459 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
460 : INTEGER :: i, handle
461 :
462 14 : CALL timeset(routineN, handle)
463 : ! Get pointers to parameters from qs_env
464 14 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
465 :
466 : ! ****** START HARTREE POTENTIAL REFERENCE CALCULATION
467 : ! Calculate Coulomb RI elements, necessary for Hartree calculation
468 14 : CALL init_hartree(rtbse_env, rtbse_env%v_dbcsr)
469 : ! Calculate the original Hartree potential
470 : ! Uses rho_orig - same as rho for initial run but different for continued run
471 28 : DO i = 1, rtbse_env%n_spin
472 14 : CALL get_hartree(rtbse_env, rtbse_env%rho_orig(i), rtbse_env%hartree_curr(i))
473 : ! Scaling by spin degeneracy
474 14 : CALL cp_fm_scale(rtbse_env%spin_degeneracy, rtbse_env%hartree_curr(i))
475 : ! Subtract the reference from the reference Hamiltonian
476 14 : CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(i), mtarget=rtbse_env%ham_workspace(1))
477 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
478 28 : CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
479 : END DO
480 : ! ****** END HARTREE POTENTIAL REFERENCE CALCULATION
481 :
482 14 : CALL timestop(handle)
483 14 : END SUBROUTINE initialize_hartree_potential
484 :
485 : ! **************************************************************************************************
486 : !> \brief Calculates the COHSEX reference self-energy
487 : !> \param rtbse_env RT-BSE environment
488 : !> \author Stepan Marek (09.24)
489 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
490 : ! **************************************************************************************************
491 14 : SUBROUTINE initialize_cohsex_selfenergy(rtbse_env)
492 : TYPE(rtbse_env_type), POINTER :: rtbse_env
493 : CHARACTER(len=*), PARAMETER :: routineN = "initialize_cohsex_selfenergy"
494 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
495 : INTEGER :: i, handle
496 :
497 14 : CALL timeset(routineN, handle)
498 : ! Get pointers to parameters from qs_env
499 14 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
500 :
501 : ! ****** START COHSEX REFERENCE CALCULATION
502 14 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
503 : ! In a non-HF calculation, copy the actual correlation part of the interaction
504 14 : CALL copy_fm_to_dbcsr(bs_env%fm_W_MIC_freq_zero, rtbse_env%w_dbcsr)
505 : ELSE
506 : ! In HF, correlation is set to zero
507 0 : CALL dbcsr_set(rtbse_env%w_dbcsr, 0.0_dp)
508 : END IF
509 : ! Add the Hartree to the screened_dbt tensor - now W = V + W^c
510 14 : CALL dbcsr_add(rtbse_env%w_dbcsr, rtbse_env%v_dbcsr, 1.0_dp, 1.0_dp)
511 14 : CALL dbt_copy_matrix_to_tensor(rtbse_env%w_dbcsr, rtbse_env%screened_dbt)
512 : ! Calculate the COHSEX starting energies
513 28 : DO i = 1, rtbse_env%n_spin
514 14 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
515 : ! Subtract the v_xc from COH part of the self-energy, as V_xc is also not updated during the timestepping
516 : ! TODO : Allow no COH calculation for static screening
517 14 : CALL get_sigma(rtbse_env, rtbse_env%sigma_COH(i), -0.5_dp, rtbse_env%S_inv_fm)
518 : ! Copy and subtract from the complex reference hamiltonian
519 14 : CALL cp_fm_to_cfm(msourcer=rtbse_env%sigma_COH(i), mtarget=rtbse_env%ham_workspace(1))
520 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
521 14 : CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
522 : END IF
523 : ! Calculate exchange part - TODO : should this be applied for different spins? - TEST with O2 HF propagation?
524 : ! So far only closed shell tested
525 : ! Uses rho_orig - same as rho for initial run but different for continued run
526 : ! For KS reference this is the time-dependent Fock exchange (w_dbcsr = v only).
527 14 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(i), -1.0_dp, rtbse_env%rho_orig(i))
528 : ! Subtract from the complex reference Hamiltonian
529 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
530 28 : CMPLX(-1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(i))
531 : END DO
532 : ! ****** END COHSEX REFERENCE CALCULATION
533 :
534 14 : CALL timestop(handle)
535 14 : END SUBROUTINE initialize_cohsex_selfenergy
536 :
537 : ! **************************************************************************************************
538 : !> \brief Custom reimplementation of the delta pulse routines
539 : !> \param rtbse_env RT-BSE environment
540 : !> \author Stepan Marek (09.24)
541 : ! **************************************************************************************************
542 6 : SUBROUTINE apply_delta_pulse(rtbse_env)
543 : TYPE(rtbse_env_type), POINTER :: rtbse_env
544 : CHARACTER(len=*), PARAMETER :: routineN = "apply_delta_pulse"
545 : REAL(kind=dp) :: intensity, metric
546 : REAL(kind=dp), DIMENSION(3) :: kvec
547 : INTEGER :: i, k, handle
548 :
549 6 : CALL timeset(routineN, handle)
550 :
551 : ! Report application
552 6 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A28)') ' RTBSE| Applying delta pulse'
553 : ! Extra minus for the propagation of density
554 6 : intensity = -rtbse_env%dft_control%rtp_control%delta_pulse_scale
555 : metric = 0.0_dp
556 24 : kvec(:) = rtbse_env%dft_control%rtp_control%delta_pulse_direction(:)
557 6 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A38,E14.4E3,E14.4E3,E14.4E3)') &
558 12 : " RTBSE| Delta pulse elements (a.u.) : ", intensity*kvec(:)
559 : ! So far no spin dependence, but can be added by different structure of delta pulse
560 6 : CALL cp_fm_set_all(rtbse_env%real_workspace(1), 0.0_dp)
561 24 : DO k = 1, 3
562 : CALL cp_fm_scale_and_add(1.0_dp, rtbse_env%real_workspace(1), &
563 24 : kvec(k), rtbse_env%moments_field(k, 1))
564 : END DO
565 : ! enforce hermiticity of the effective Hamiltonian
566 6 : CALL cp_fm_transpose(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2))
567 : CALL cp_fm_scale_and_add(0.5_dp, rtbse_env%real_workspace(1), &
568 6 : 0.5_dp, rtbse_env%real_workspace(2))
569 : ! Prepare the exponential/exponent for propagation
570 6 : IF (rtbse_env%mat_exp_method == do_bch) THEN
571 : ! Multiply by the S_inv matrix - in the classic ordering
572 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
573 : intensity, rtbse_env%S_inv_fm, rtbse_env%real_workspace(1), &
574 6 : 0.0_dp, rtbse_env%real_workspace(2))
575 12 : DO i = 1, rtbse_env%n_spin
576 : ! Sets real part to zero
577 12 : CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_workspace(i))
578 : END DO
579 0 : ELSE IF (rtbse_env%mat_exp_method == do_exact) THEN
580 0 : DO i = 1, rtbse_env%n_spin
581 0 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(1), mtarget=rtbse_env%ham_effective(i))
582 : CALL cp_cfm_gexp(rtbse_env%ham_effective(i), rtbse_env%S_cfm, rtbse_env%ham_workspace(i), &
583 0 : CMPLX(0.0, intensity, kind=dp), rtbse_env%rho_workspace)
584 : END DO
585 : END IF
586 : ! Propagate the density by the effect of the delta pulse
587 6 : CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rtbse_env%rho, rtbse_env%rho_new)
588 6 : metric = rho_metric(rtbse_env%rho_new, rtbse_env%rho, rtbse_env%n_spin)
589 6 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, ('(A42,E38.8E3)')) " RTBSE| Metric difference after delta kick", metric
590 : ! Copy the new density to the old density
591 12 : DO i = 1, rtbse_env%n_spin
592 12 : CALL cp_cfm_to_cfm(rtbse_env%rho_new(i), rtbse_env%rho(i))
593 : END DO
594 :
595 6 : CALL timestop(handle)
596 6 : END SUBROUTINE apply_delta_pulse
597 : ! **************************************************************************************************
598 : !> \brief Determines the metric for the density matrix, used for convergence criterion
599 : !> \param rho_new Array of new density matrices (one for each spin index)
600 : !> \param rho_old Array of old density matrices (one for each spin index)
601 : !> \param nspin Number of spin indices
602 : !> \param workspace_opt Optionally provide external workspace to save some allocation time
603 : ! **************************************************************************************************
604 22438 : FUNCTION rho_metric(rho_new, rho_old, nspin, workspace_opt) RESULT(metric)
605 : TYPE(cp_cfm_type), DIMENSION(:), POINTER, INTENT(IN):: rho_new, &
606 : rho_old
607 : INTEGER, INTENT(IN) :: nspin
608 : TYPE(cp_cfm_type), POINTER, OPTIONAL :: workspace_opt
609 : TYPE(cp_cfm_type) :: workspace
610 : REAL(kind=dp) :: metric
611 22438 : REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: partial_metric
612 : INTEGER :: j
613 : COMPLEX(kind=dp) :: scale_factor
614 :
615 67314 : ALLOCATE (partial_metric(nspin))
616 :
617 : ! Only allocate/deallocate storage if required
618 22438 : IF (PRESENT(workspace_opt)) THEN
619 0 : workspace = workspace_opt
620 : ELSE
621 22438 : CALL cp_cfm_create(workspace, rho_new(1)%matrix_struct)
622 : END IF
623 22438 : scale_factor = 1.0
624 44900 : DO j = 1, nspin
625 22462 : CALL cp_cfm_to_cfm(rho_new(j), workspace)
626 : ! Get the difference in the resulting matrix
627 22462 : CALL cp_cfm_scale_and_add(scale_factor, workspace, -scale_factor, rho_old(j))
628 : ! Now, get the relevant number
629 44900 : partial_metric(j) = cp_cfm_norm(workspace, 'M')
630 : END DO
631 : metric = 0.0_dp
632 : ! For more than one spin, do Cartesian sum of the different spin norms
633 44900 : DO j = 1, nspin
634 44900 : metric = metric + partial_metric(j)*partial_metric(j)
635 : END DO
636 22438 : metric = SQRT(metric)
637 : ! Deallocate workspace
638 22438 : IF (.NOT. PRESENT(workspace_opt)) CALL cp_cfm_release(workspace)
639 22438 : DEALLOCATE (partial_metric)
640 22438 : END FUNCTION rho_metric
641 :
642 : ! **************************************************************************************************
643 : !> \brief Determines the metric of the antihermitian part of the matrix
644 : !> \param real_fm Real part of the full matrix
645 : !> \param imag_fm Imaginary part of the full matrix
646 : ! **************************************************************************************************
647 0 : SUBROUTINE antiherm_metric(real_fm, imag_fm, workspace, metric)
648 : TYPE(cp_fm_type), INTENT(IN) :: real_fm
649 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: imag_fm
650 : REAL(kind=dp), INTENT(OUT) :: metric
651 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: workspace
652 : COMPLEX(kind=dp) :: complex_one
653 :
654 : ! Get the complex and complex conjugate matrix
655 0 : IF (PRESENT(imag_fm)) THEN
656 0 : CALL cp_fm_to_cfm(real_fm, imag_fm, workspace(1))
657 : ELSE
658 0 : CALL cp_fm_to_cfm(msourcer=real_fm, mtarget=workspace(1))
659 : END IF
660 0 : CALL cp_cfm_transpose(workspace(1), "C", workspace(2))
661 : ! Subtract these, and get the metric
662 0 : complex_one = CMPLX(1.0, 0.0, kind=dp)
663 0 : CALL cp_cfm_scale_and_add(complex_one, workspace(1), -complex_one, workspace(2))
664 0 : metric = cp_cfm_norm(workspace(1), "M")
665 0 : END SUBROUTINE antiherm_metric
666 :
667 : ! **************************************************************************************************
668 : !> \brief For Taylor and Exact exp_method, calculates the matrix exponential of the
669 : !> effective Hamiltonian. For BCH, calculates just the effective Hamiltonian. For other methods,
670 : !> aborts the execution, as they are not implemented yet.
671 : !> \param rtbse_env Entry point of the calculation. Uses rho_workspace for Taylor and BCH. For exact,
672 : !> uses complex_workspace, complex_ham, complex_s, real_eigvals and exp_eigvals.
673 : !> Results are stored in ham_workspace.
674 : ! **************************************************************************************************
675 2410 : SUBROUTINE ham_to_exp(rtbse_env, ham, ham_exp)
676 : TYPE(rtbse_env_type), POINTER :: rtbse_env
677 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham, &
678 : ham_exp
679 : CHARACTER(len=*), PARAMETER :: routineN = "ham_to_exp"
680 : INTEGER :: j, handle
681 2410 : CALL timeset(routineN, handle)
682 4820 : DO j = 1, rtbse_env%n_spin
683 4820 : IF (rtbse_env%mat_exp_method == do_bch) THEN
684 : ! In Taylor and BCH, we first evaluate the entire exponent and then evaluate exponential in series
685 : ! In order to produce correct result, need to remultiply by inverse overlap matrix
686 : CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
687 : 1.0_dp, rtbse_env%S_inv_fm, ham(j), &
688 2410 : 0.0_dp, rtbse_env%rho_workspace(1))
689 :
690 : ! The evolution of density matrix is derived from the right multiplying term
691 : ! Imaginary part of the exponent = -real part of the matrix
692 2410 : CALL cp_cfm_scale(CMPLX(0.0, -rtbse_env%sim_dt/2, kind=dp), rtbse_env%rho_workspace(1))
693 : ! In BCH, exponential is not calculated explicitly, but the propagation is solved in series
694 2410 : CALL cp_cfm_to_cfm(rtbse_env%rho_workspace(1), ham_exp(j))
695 0 : ELSE IF (rtbse_env%mat_exp_method == do_exact) THEN
696 : CALL cp_cfm_gexp(ham(j), rtbse_env%S_cfm, ham_exp(j), &
697 0 : CMPLX(0.0, -rtbse_env%sim_dt/2, kind=dp), rtbse_env%rho_workspace)
698 : ELSE
699 0 : CPABORT("Only BCH and Taylor matrix exponentiation implemented")
700 : END IF
701 : END DO
702 :
703 2410 : CALL timestop(handle)
704 2410 : END SUBROUTINE ham_to_exp
705 : ! **************************************************************************************************
706 : !> \brief Updates the effective Hamiltonian, given a density matrix rho
707 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
708 : !> \param qs_env QS env
709 : !> \param rho Real and imaginary parts ( + spin) of the density at current time
710 : ! **************************************************************************************************
711 2410 : SUBROUTINE update_effective_ham(rtbse_env, rho)
712 : TYPE(rtbse_env_type), POINTER :: rtbse_env
713 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
714 : CHARACTER(len=*), PARAMETER :: routineN = "update_effective_ham"
715 : INTEGER :: k, j, nspin, handle
716 :
717 2410 : CALL timeset(routineN, handle)
718 : ! Shorthand
719 2410 : nspin = rtbse_env%n_spin
720 : ! Reset the effective Hamiltonian to KS Hamiltonian + G0W0 - reference COHSEX - reference Hartree
721 4820 : DO j = 1, nspin
722 : ! Sets the imaginary part to zero
723 4820 : CALL cp_cfm_to_cfm(rtbse_env%ham_reference(j), rtbse_env%ham_effective(j))
724 : END DO
725 : ! Determine the field at current time
726 2410 : IF (rtbse_env%dft_control%apply_efield_field) THEN
727 632 : CALL make_field(rtbse_env%dft_control, rtbse_env%field, rtbse_env%sim_step, rtbse_env%sim_time)
728 : ELSE
729 : ! No field
730 7112 : rtbse_env%field(:) = 0.0_dp
731 : END IF
732 4820 : DO j = 1, nspin
733 9640 : DO k = 1, 3
734 : ! Minus sign due to charge of electrons
735 7230 : CALL cp_fm_to_cfm(msourcer=rtbse_env%moments_field(k, 1), mtarget=rtbse_env%ham_workspace(1))
736 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
737 9640 : CMPLX(rtbse_env%field(k), 0.0, kind=dp), rtbse_env%ham_workspace(1))
738 : END DO
739 2410 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN
740 : ! Add the COH part - so far static but can be dynamic in principle through the W updates
741 2410 : CALL get_sigma(rtbse_env, rtbse_env%sigma_COH(j), -0.5_dp, rtbse_env%S_inv_fm)
742 2410 : CALL cp_fm_to_cfm(msourcer=rtbse_env%sigma_COH(j), mtarget=rtbse_env%ham_workspace(1))
743 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
744 2410 : CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
745 : END IF
746 : ! Calculate the (S)EX part - based on provided rho
747 : ! iGW = - rho W
748 2410 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(j), -1.0_dp, rho(j))
749 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
750 2410 : CMPLX(1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(j))
751 : ! Calculate Hartree potential
752 : ! Hartree potential is scaled by number of electrons in each MO - spin degeneracy
753 : CALL get_hartree(rtbse_env, rho(j), &
754 2410 : rtbse_env%hartree_curr(j))
755 2410 : CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(j), mtarget=rtbse_env%ham_workspace(1))
756 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), &
757 2410 : CMPLX(rtbse_env%spin_degeneracy, 0.0, kind=dp), rtbse_env%ham_workspace(1))
758 : ! Enforce hermiticity of the effective Hamiltonian
759 : ! Important components without forced Hermiticity - moments matrix, sigma matrices, Hartree matrix
760 : ! single particle Ham
761 2410 : CALL cp_cfm_transpose(rtbse_env%ham_effective(j), 'C', rtbse_env%ham_workspace(1))
762 : CALL cp_cfm_scale_and_add(CMPLX(0.5, 0.0, kind=dp), rtbse_env%ham_effective(j), &
763 4820 : CMPLX(0.5, 0.0, kind=dp), rtbse_env%ham_workspace(1))
764 : END DO
765 2410 : CALL timestop(handle)
766 2410 : END SUBROUTINE update_effective_ham
767 : ! **************************************************************************************************
768 : !> \brief Self-consistently (ETRS) propagates the density to the next timestep
769 : !> \note Uses rtbse_env%rho_new_last, assumes correct timestep information is given in rtbse_env
770 : !> \param rho_start Initial density matrix
771 : !> \param rho_mid Midpoint density (propagated to by the initial Hamiltonian)
772 : !> \param rho_end Endpoint density (propagated to by endpoint Hamiltonian)
773 : !> \param converged Whether the resulting rho_end is self-consistent
774 : !> \param k How many SC iterations were done
775 : !> \param metric The difference metric from the last self-consistent iteration (for printing/evaluation)
776 : ! **************************************************************************************************
777 678 : SUBROUTINE etrs_scf_loop(rtbse_env, rho_start, rho_mid, rho_end, converged, k, metric)
778 : TYPE(rtbse_env_type), POINTER :: rtbse_env
779 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_start, &
780 : rho_mid, &
781 : rho_end
782 : LOGICAL :: converged
783 : INTEGER :: k
784 : REAL(kind=dp) :: metric
785 : CHARACTER(len=*), PARAMETER :: routineN = "etrs_scf_loop"
786 : INTEGER :: j, handle
787 :
788 678 : CALL timeset(routineN, handle)
789 :
790 : ! This method determines the density matrix at time (t+dt) by guessing the effective Hamiltonian at (t + dt)
791 : ! and using the Hamiltonian at time (t), it propagates density from time (t) while ensuring that the density
792 : ! at (t + dt/2) is the same for both forward and backwards propagation. Then, density at (t + dt) is
793 : ! used to calculate the new Hamiltonian at (t+dt), which is then used to get the new propagator, and so on
794 : ! until the density matrix does not change within certain limit
795 : ! Pseudocode of the algorithm
796 : ! rho_M = exp(-i S^(-1) H[rho(t)] dt/2) rho(t) exp(i H[rho(t)] S^(-1) dt/2)
797 : ! rho(t+dt, 0) = rho_M
798 : ! for j in 1,max_self_iter
799 : ! rho(t+dt,j) = exp(- i S^(-1) H[rho(t+dt,j-1)] dt/2) rho_M exp(i H [rho(t+dt,j-1)] S^(-1) dt/2)
800 : ! if ||rho(t+dt,j) - rho(t+dt,j-1)|| < epsilon
801 : ! break
802 :
803 : ! Initial setup - calculate the Hamiltonian
804 678 : CALL update_effective_ham(rtbse_env, rho_start)
805 : ! Create the exponential
806 678 : CALL ham_to_exp(rtbse_env, rtbse_env%ham_effective, rtbse_env%ham_workspace)
807 : ! Propagate to rho_mid
808 678 : CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rho_start, rho_mid)
809 : ! Propagate to initial guess
810 678 : CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rho_mid, rtbse_env%rho_new_last)
811 : ! Update bookkeeping to the next timestep - Hamiltonians are now evaluated at the next timestep
812 678 : rtbse_env%sim_step = rtbse_env%sim_step + 1
813 678 : rtbse_env%sim_time = rtbse_env%sim_time + rtbse_env%sim_dt
814 678 : converged = .FALSE.
815 678 : CALL print_etrs_info_header(rtbse_env)
816 1732 : DO k = 1, rtbse_env%etrs_max_iter
817 : ! Get the Hamiltonian following from the last timestep
818 1732 : CALL update_effective_ham(rtbse_env, rtbse_env%rho_new_last)
819 1732 : CALL ham_to_exp(rtbse_env, rtbse_env%ham_effective, rtbse_env%ham_workspace)
820 : ! Propagate to new guess
821 1732 : CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rho_mid, rho_end)
822 : ! Check for self-consistency
823 1732 : metric = rho_metric(rho_end, rtbse_env%rho_new_last, rtbse_env%n_spin)
824 : ! ETRS info - only for log level > medium
825 1732 : CALL print_etrs_info(rtbse_env, k, metric)
826 1732 : IF (metric < rtbse_env%etrs_threshold) THEN
827 678 : converged = .TRUE.
828 678 : EXIT
829 : ELSE
830 : ! Copy rho_new to rho_new_last
831 2108 : DO j = 1, rtbse_env%n_spin
832 : ! Leaving for free convergence
833 2108 : CALL cp_cfm_to_cfm(rho_end(j), rtbse_env%rho_new_last(j))
834 : END DO
835 : END IF
836 : END DO
837 : ! Error handling in the case where the propagation did not converge is left to the main routine
838 678 : CALL timestop(handle)
839 678 : END SUBROUTINE etrs_scf_loop
840 :
841 : ! **************************************************************************************************
842 : !> \brief Does the BCH iterative determination of the exponential
843 : !> \param propagator_matrix Matrix X which is to be exponentiated
844 : !> \param target_matrix Matrix Y which the exponential acts upon
845 : !> \param result_matrix Propagated matrix
846 : !> \param workspace Matrices dedicated for work, 4 fm matrices with dimensions of X required
847 : !> \param threshold_opt Optionally, a threshold under which the iteration is considered converged (default 1e-10)
848 : !> \param max_iter_opt Optionally, maximum number of BCH iterations (default 20)
849 : ! **************************************************************************************************
850 6264 : SUBROUTINE bch_propagate(propagator_matrix, target_matrix, result_matrix, workspace, threshold_opt, max_iter_opt)
851 : ! Array of complex propagator matrix X, such that
852 : ! the propagated matrix will follow Y' = e^X Y e^(-X), for each spin
853 : ! effect of e^(-X) is calculated - provide the X on the left hand side
854 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: propagator_matrix
855 : ! Matrix Y to be propagated into matrix Y'
856 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: target_matrix
857 : ! Matrix Y' is stored here on exit
858 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: result_matrix, workspace
859 : ! Threshold for the metric which decides when to truncate the BCH expansion
860 : REAL(kind=dp), OPTIONAL :: threshold_opt
861 : INTEGER, OPTIONAL :: max_iter_opt
862 : CHARACTER(len=*), PARAMETER :: routineN = "bch_propagate"
863 : REAL(kind=dp) :: threshold, prefactor, metric
864 : INTEGER :: max_iter, i, n_spin, n_ao, k, &
865 : w_stride, handle
866 : LOGICAL :: converged
867 : CHARACTER(len=77) :: error
868 :
869 3132 : CALL timeset(routineN, handle)
870 :
871 3132 : converged = .FALSE.
872 :
873 3132 : IF (PRESENT(threshold_opt)) THEN
874 3132 : threshold = threshold_opt
875 : ELSE
876 0 : threshold = 1.0e-10
877 : END IF
878 :
879 3132 : IF (PRESENT(max_iter_opt)) THEN
880 3132 : max_iter = max_iter_opt
881 : ELSE
882 : max_iter = 20
883 : END IF
884 :
885 3132 : n_spin = SIZE(target_matrix)
886 : n_ao = 0
887 3132 : CALL cp_cfm_get_info(target_matrix(1), nrow_global=n_ao)
888 3132 : w_stride = n_spin
889 :
890 : ! Initiate
891 6270 : DO i = 1, n_spin
892 3138 : CALL cp_cfm_to_cfm(target_matrix(i), result_matrix(i))
893 6270 : CALL cp_cfm_to_cfm(target_matrix(i), workspace(i))
894 : END DO
895 :
896 : ! Start the BCH iterations
897 : ! So far, no spin mixing terms
898 20662 : DO k = 1, max_iter
899 20662 : prefactor = 1.0_dp/REAL(k, kind=dp)
900 41342 : DO i = 1, n_spin
901 : CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, &
902 : CMPLX(prefactor, 0.0, kind=dp), propagator_matrix(i), workspace(i), &
903 20680 : CMPLX(0.0, 0.0, kind=dp), workspace(i + w_stride))
904 : CALL parallel_gemm("N", "C", n_ao, n_ao, n_ao, &
905 : CMPLX(prefactor, 0.0, kind=dp), workspace(i), propagator_matrix(i), &
906 20680 : CMPLX(1.0, 0.0, kind=dp), workspace(i + w_stride))
907 : ! Add to the result
908 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), result_matrix(i), &
909 41342 : CMPLX(1.0, 0.0, kind=dp), workspace(i + w_stride))
910 : END DO
911 20662 : metric = rho_metric(workspace(w_stride + 1:), workspace(1:w_stride), n_spin)
912 20662 : IF (metric <= threshold) THEN
913 : converged = .TRUE.
914 : EXIT
915 : ELSE
916 35072 : DO i = 1, n_spin
917 35072 : CALL cp_cfm_to_cfm(workspace(i + w_stride), workspace(i))
918 : END DO
919 : END IF
920 : END DO
921 3132 : IF (.NOT. converged) THEN
922 0 : WRITE (error, '(A35,E13.4E3,A16,E13.4E3)') "BCH did not converge, BCH Metric : ", &
923 0 : metric, "BCH Threshold : ", threshold
924 0 : CPABORT(error)
925 : END IF
926 :
927 3132 : CALL timestop(handle)
928 3132 : END SUBROUTINE bch_propagate
929 :
930 : ! **************************************************************************************************
931 : !> \brief Updates the density in rtbse_env, using the provided exponential
932 : !> The new density is saved to a different matrix, which enables for comparison of matrices
933 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
934 : !> \param exponential Real and imaginary parts ( + spin) of the exponential propagator
935 : ! **************************************************************************************************
936 3132 : SUBROUTINE propagate_density(rtbse_env, exponential, rho_old, rho_new)
937 : TYPE(rtbse_env_type) :: rtbse_env
938 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: exponential, &
939 : rho_old, &
940 : rho_new
941 : CHARACTER(len=*), PARAMETER :: routineN = "propagate_density"
942 : INTEGER :: j, handle
943 :
944 3132 : CALL timeset(routineN, handle)
945 3132 : IF (rtbse_env%mat_exp_method == do_exact) THEN
946 : ! For these methods, exponential is explicitly constructed
947 0 : DO j = 1, rtbse_env%n_spin
948 : ! rho * (exp^dagger)
949 : CALL parallel_gemm("N", "C", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
950 : CMPLX(1.0, 0.0, kind=dp), rho_old(j), exponential(j), &
951 0 : CMPLX(0.0, 0.0, kind=dp), rtbse_env%rho_workspace(1))
952 : ! exp * rho * (exp^dagger)
953 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
954 : CMPLX(1.0, 0.0, kind=dp), exponential(j), rtbse_env%rho_workspace(1), &
955 0 : CMPLX(0.0, 0.0, kind=dp), rho_new(j))
956 : END DO
957 3132 : ELSE IF (rtbse_env%mat_exp_method == do_bch .OR. rtbse_env%linearized) THEN
958 : ! Same number of iterations as ETRS
959 : CALL bch_propagate(exponential, rho_old, rho_new, rtbse_env%rho_workspace, threshold_opt=rtbse_env%exp_accuracy, &
960 3132 : max_iter_opt=rtbse_env%etrs_max_iter)
961 : ELSE
962 0 : CPABORT("Only BCH and exact matrix exponentiation implemented.")
963 : END IF
964 :
965 3132 : CALL timestop(handle)
966 3132 : END SUBROUTINE propagate_density
967 :
968 : ! **************************************************************************************************
969 : !> \brief Outputs the number of electrons in the system from the density matrix
970 : !> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
971 : !> \param rtbse_env Entry point - rtbse environment
972 : !> \param rho Density matrix in AO basis
973 : !> \param electron_n_re Real number of electrons
974 : !> \param electron_n_im Imaginary number of electrons, which can arise from numerical non-hermiticity
975 : ! **************************************************************************************************
976 678 : SUBROUTINE get_electron_number(rtbse_env, rho, electron_n_re, electron_n_im)
977 : TYPE(rtbse_env_type) :: rtbse_env
978 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
979 : REAL(kind=dp), INTENT(OUT) :: electron_n_re, electron_n_im
980 : COMPLEX(kind=dp) :: electron_n_buffer
981 : INTEGER :: j
982 :
983 678 : electron_n_re = 0.0_dp
984 678 : electron_n_im = 0.0_dp
985 678 : CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%rho_workspace(1))
986 1356 : DO j = 1, rtbse_env%n_spin
987 678 : CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rho(j), electron_n_buffer)
988 678 : electron_n_re = electron_n_re + REAL(electron_n_buffer, kind=dp)
989 1356 : electron_n_im = electron_n_im + REAL(AIMAG(electron_n_buffer), kind=dp)
990 : END DO
991 : ! Scale by spin degeneracy
992 678 : electron_n_re = electron_n_re*rtbse_env%spin_degeneracy
993 678 : electron_n_im = electron_n_im*rtbse_env%spin_degeneracy
994 678 : END SUBROUTINE get_electron_number
995 : ! **************************************************************************************************
996 : !> \brief Outputs the deviation from idempotence of density matrix
997 : !> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3)
998 : !> \param rtbse_env Entry point - rtbse environment
999 : !> \param rho Density matrix in AO basis
1000 : !> \param electron_n_re Real number of electrons
1001 : !> \param electron_n_im Imaginary number of electrons, which can arise from numerical non-hermiticity
1002 : ! **************************************************************************************************
1003 0 : SUBROUTINE get_idempotence_deviation(rtbse_env, rho, deviation_metric)
1004 : TYPE(rtbse_env_type) :: rtbse_env
1005 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
1006 : REAL(kind=dp), INTENT(OUT) :: deviation_metric
1007 : COMPLEX(kind=dp) :: buffer_1, buffer_2
1008 : REAL(kind=dp) :: buffer_dev
1009 : INTEGER :: j
1010 :
1011 0 : deviation_metric = 0.0_dp
1012 0 : buffer_dev = 0.0_dp
1013 : ! First, determine Tr(S * rho_re) + i Tr (S * rho_im)
1014 0 : CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%rho_workspace(1))
1015 0 : DO j = 1, rtbse_env%n_spin
1016 0 : CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rho(j), buffer_1)
1017 0 : buffer_dev = buffer_dev + REAL(ABS(buffer_1)*ABS(buffer_1), kind=dp)
1018 : END DO
1019 : ! Now, determine Tr(S * rho_re * S * rho_re) - Tr(S * rho_im * S * rho_im) + 2i Tr(S * rho_re * S * rho_im)
1020 0 : DO j = 1, rtbse_env%n_spin
1021 : ! S * rho
1022 : CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
1023 : 1.0_dp, rtbse_env%S_fm, rho(j), &
1024 0 : 0.0_dp, rtbse_env%rho_workspace(2))
1025 : ! rho * S * rho
1026 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, &
1027 : CMPLX(1.0, 0.0, kind=dp), rho(j), rtbse_env%rho_workspace(2), &
1028 0 : CMPLX(0.0, 0.0, kind=dp), rtbse_env%rho_workspace(3))
1029 : ! Tr (S * rho * S * rho)
1030 0 : CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rtbse_env%rho_workspace(3), buffer_2)
1031 0 : deviation_metric = deviation_metric + REAL(ABS(buffer_2)*ABS(buffer_2), kind=dp)
1032 : END DO
1033 0 : deviation_metric = SQRT(deviation_metric) - SQRT(buffer_dev)
1034 0 : END SUBROUTINE get_idempotence_deviation
1035 :
1036 : ! **************************************************************************************************
1037 : !> \brief Calculates the self-energy by contraction of screened potential, for complex density
1038 : !> \note Can be used for both the Coulomb hole part and screened exchange part
1039 : !> \param rtbse_env Quickstep environment data, entry point of the calculation
1040 : !> \param sigma_cfm Pointer to the self-energy full matrix, which is overwritten by this routine
1041 : !> \param prefactor_opt Optional scaling factor applied to the contraction, defaults to 1.0
1042 : !> \param greens_cfm Pointer to the Green's function matrix, which is used as input data
1043 : !> \param grid_diag_re_accum Optional accumulator for the real part of the RI-RS grid diagonal
1044 : !> \param grid_diag_im_accum Optional accumulator for the imaginary part of the RI-RS grid diagonal
1045 : !> \author Stepan Marek
1046 : !> \date 09.2024
1047 : ! **************************************************************************************************
1048 6216 : SUBROUTINE get_sigma_complex(rtbse_env, sigma_cfm, prefactor_opt, greens_cfm, &
1049 6216 : grid_diag_re_accum, grid_diag_im_accum)
1050 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1051 : TYPE(cp_cfm_type) :: sigma_cfm ! resulting self energy
1052 : REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt
1053 : TYPE(cp_cfm_type), INTENT(IN) :: greens_cfm ! matrix to contract with RI_W
1054 : REAL(kind=dp), INTENT(INOUT), OPTIONAL :: grid_diag_re_accum(:), grid_diag_im_accum(:)
1055 : REAL(kind=dp) :: prefactor
1056 :
1057 6216 : prefactor = 1.0_dp
1058 6216 : IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt
1059 :
1060 : ! RI-RS screened-exchange backend (linRTBSE only; rirs_kernel is forced .FALSE. for full RTBSE,
1061 : ! so this is inert there). The RI-RS routine does its own Re/Im split, replacing the AO-RI body.
1062 : ! The optional grid_diag_* accumulators harvest diag(φρφ^T) for the Hartree reuse (RI-RS only;
1063 : ! absent on AO-RI calls, which build no grid).
1064 6216 : IF (rtbse_env%rirs_kernel) THEN
1065 : CALL compute_sigma_ri_rs_complex(rtbse_env%bs_env, sigma_cfm, prefactor, greens_cfm, &
1066 : grid_diag_re_accum=grid_diag_re_accum, &
1067 3016 : grid_diag_im_accum=grid_diag_im_accum)
1068 1752 : RETURN
1069 : END IF
1070 :
1071 : ! Carry out the sigma part twice
1072 : ! Real part
1073 4464 : CALL cp_cfm_to_fm(msource=greens_cfm, mtargetr=rtbse_env%real_workspace(1))
1074 4464 : CALL get_sigma(rtbse_env, rtbse_env%real_workspace(2), prefactor, rtbse_env%real_workspace(1))
1075 4464 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%sigma_complex_workspace(1))
1076 : ! Imaginary part
1077 4464 : CALL cp_cfm_to_fm(msource=greens_cfm, mtargeti=rtbse_env%real_workspace(1))
1078 4464 : CALL get_sigma(rtbse_env, rtbse_env%real_workspace(2), prefactor, rtbse_env%real_workspace(1))
1079 4464 : CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace(2), mtarget=sigma_cfm)
1080 : ! Add the real part
1081 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), sigma_cfm, &
1082 4464 : CMPLX(1.0, 0.0, kind=dp), rtbse_env%sigma_complex_workspace(1))
1083 :
1084 : END SUBROUTINE get_sigma_complex
1085 : ! **************************************************************************************************
1086 : !> \brief Calculates the self-energy by contraction of screened potential, for complex density
1087 : !> \note Can be used for both the Coulomb hole part and screened exchange part
1088 : !> \param rtbse_env Quickstep environment data, entry point of the calculation
1089 : !> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine
1090 : !> \param greens_fm Pointer to the Green's function matrix, which is used as input data
1091 : !> \author Stepan Marek
1092 : !> \date 09.2024
1093 : ! **************************************************************************************************
1094 11352 : SUBROUTINE get_sigma_real(rtbse_env, sigma_fm, prefactor_opt, greens_fm)
1095 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1096 : TYPE(cp_fm_type) :: sigma_fm ! resulting self energy
1097 : REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt
1098 : TYPE(cp_fm_type), INTENT(IN) :: greens_fm ! matrix to contract with RI_W
1099 : REAL(kind=dp) :: prefactor
1100 : TYPE(dbcsr_type) :: greens_dbcsr_scratch
1101 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1102 :
1103 11352 : prefactor = 1.0_dp
1104 11352 : IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt
1105 :
1106 : ! Local AO-AO dbcsr scratch for the FM->DBCSR conversion. Previously this
1107 : ! routine used rtbse_env%rho_dbcsr as the workspace, which coupled AO-RI SX
1108 : ! to the AO-RI Hartree allocation path - the historical (RIRS-H + AO-RI-SX)
1109 : ! cross-combo (no longer expressible under the single KERNEL_RI switch)
1110 : ! then segfaulted because rho_dbcsr is skipped when rirs_kernel=T.
1111 11352 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
1112 : CALL dbcsr_create(greens_dbcsr_scratch, name="get_sigma greens scratch", &
1113 11352 : template=bs_env%mat_ao_ao%matrix)
1114 :
1115 11352 : CALL copy_fm_to_dbcsr(greens_fm, greens_dbcsr_scratch)
1116 11352 : CALL get_sigma_dbcsr(rtbse_env, sigma_fm, prefactor, greens_dbcsr_scratch)
1117 :
1118 11352 : CALL dbcsr_release(greens_dbcsr_scratch)
1119 11352 : END SUBROUTINE get_sigma_real
1120 : ! **************************************************************************************************
1121 : !> \brief Calculates the self-energy by contraction of screened potential
1122 : !> \note Can be used for both the Coulomb hole part and screened exchange part
1123 : !> \param greens_fm Pointer to the Green's function matrix, which is used as input data
1124 : !> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine
1125 : !> \author Stepan Marek
1126 : !> \date 01.2024
1127 : ! **************************************************************************************************
1128 11352 : SUBROUTINE get_sigma_dbcsr(rtbse_env, sigma_fm, prefactor_opt, greens_dbcsr)
1129 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1130 : TYPE(cp_fm_type) :: sigma_fm ! resulting self energy
1131 : REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt
1132 : TYPE(dbcsr_type) :: greens_dbcsr
1133 : REAL(kind=dp) :: prefactor
1134 :
1135 11352 : prefactor = 1.0_dp
1136 11352 : IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt
1137 :
1138 : CALL get_sigma_noenv(sigma_fm, prefactor_opt, greens_dbcsr, &
1139 : rtbse_env%screened_dbt, rtbse_env%t_3c_w, &
1140 : rtbse_env%t_3c_work_RI_AO__AO, rtbse_env%t_3c_work2_RI_AO__AO, &
1141 11352 : rtbse_env%greens_dbt)
1142 11352 : END SUBROUTINE get_sigma_dbcsr
1143 : ! **************************************************************************************************
1144 : !> \brief Calculates the self-energy by contraction of screened potential
1145 : !> \note Can be used for both the Coulomb hole part and screened exchange part
1146 : !> \note Separated from the rtbse_env - can be in principle called outside of the RTBSE code
1147 : !> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine
1148 : !> \param prefactor_opt Optional argument for the prefactor (used for Coulomb hole calculation)
1149 : !> \param greens_dbcsr Matrix storing the lesser Green's function elements
1150 : !> \param screened_dbt Tensor storing the W_PQ screened Coulomb interaction RI matrix elements
1151 : !> \param int_3c_dbt Tensor storing the 3c integrals (RI| ORB ORB )
1152 : !> \param work_dbt_3c_1 Tensor workspace optimised for RI_AO__AO contractions
1153 : !> \param work_dbt_3c_2 Tensor workspace optimised for RI_AO__AO contractions
1154 : !> \param work_dbt_2c Tensor workspace for 2c integrals (Green's function and self-energy)
1155 : !> \author Stepan Marek
1156 : !> \date 01.2025
1157 : ! **************************************************************************************************
1158 11352 : SUBROUTINE get_sigma_noenv(sigma_fm, prefactor_opt, greens_dbcsr, screened_dbt, &
1159 : int_3c_dbt, work_dbt_3c_1, work_dbt_3c_2, work_dbt_2c)
1160 : TYPE(cp_fm_type) :: sigma_fm ! resulting self energy
1161 : REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt
1162 : TYPE(dbcsr_type) :: greens_dbcsr
1163 : TYPE(dbt_type) :: screened_dbt, &
1164 : int_3c_dbt, &
1165 : work_dbt_3c_1, &
1166 : work_dbt_3c_2, &
1167 : work_dbt_2c
1168 : CHARACTER(len=*), PARAMETER :: routineN = 'get_sigma'
1169 : REAL(kind=dp) :: prefactor
1170 : TYPE(dbcsr_type) :: sigma_dbcsr
1171 : INTEGER :: handle
1172 :
1173 11352 : CALL timeset(routineN, handle)
1174 :
1175 11352 : IF (PRESENT(prefactor_opt)) THEN
1176 11352 : prefactor = prefactor_opt
1177 : ELSE
1178 0 : prefactor = 1.0_dp
1179 : END IF
1180 :
1181 : ! Three-centre integrals are obtained from build_3c_integrals, from qs_tensors
1182 : ! These should use sparcity, while W and Sigma can be full matrices
1183 : ! The summation is carried out by dbt library - dbt_contract in dbt_api
1184 : ! The building of the tensors might be a bit hard, because it requires a lot of parallel information
1185 : ! Probably just use the tensors already present in bs_env? They seem to be mostly work tensors
1186 : ! Create by template
1187 : CALL dbt_contract(alpha=1.0_dp, &
1188 : tensor_1=screened_dbt, &
1189 : tensor_2=int_3c_dbt, &
1190 : beta=0.0_dp, &
1191 : tensor_3=work_dbt_3c_1, &
1192 : contract_1=[2], notcontract_1=[1], map_1=[1], &
1193 11352 : contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3])!,&
1194 : !filter_eps=bs_env%eps_filter)
1195 : ! t_work1 now contains B^P_(nu beta) = sum _ Q W _ (PQ) (iomega = 0) (Q| nu beta)
1196 : ! Next step is to convert the greens full matrix to dbcsr matrix
1197 11352 : CALL dbt_copy_matrix_to_tensor(greens_dbcsr, work_dbt_2c)
1198 : ! Then contract it
1199 : ! no scaling applied - this has to be applied externally
1200 : CALL dbt_contract(alpha=1.0_dp, &
1201 : tensor_1=work_dbt_3c_1, &
1202 : tensor_2=work_dbt_2c, &
1203 : beta=0.0_dp, &
1204 : tensor_3=work_dbt_3c_2, &
1205 : contract_1=[2], notcontract_1=[1, 3], map_1=[1, 3], &
1206 11352 : contract_2=[2], notcontract_2=[1], map_2=[2])
1207 : ! workspace 2 now contains C ^ P _ (mu beta) sum _ nu B ^ P _ (nu beta) g _ (mu nu)
1208 : CALL dbt_contract(alpha=prefactor, &
1209 : tensor_1=int_3c_dbt, &
1210 : tensor_2=work_dbt_3c_2, &
1211 : beta=0.0_dp, &
1212 : tensor_3=work_dbt_2c, &
1213 : contract_1=[1, 3], notcontract_1=[2], map_1=[1], &
1214 11352 : contract_2=[1, 2], notcontract_2=[3], map_2=[2])!,&
1215 : !filter_eps=bs_env%eps_filter)
1216 : ! Finally, convert the COH tensor to matrix and then to fm matrix
1217 : ! TODO : extra workspace?
1218 11352 : CALL dbcsr_create(sigma_dbcsr, name="sigma", template=greens_dbcsr)
1219 11352 : CALL dbt_copy_tensor_to_matrix(work_dbt_2c, sigma_dbcsr)
1220 11352 : CALL copy_dbcsr_to_fm(sigma_dbcsr, sigma_fm)
1221 11352 : CALL dbcsr_release(sigma_dbcsr)
1222 : ! Clear workspaces - saves memory?
1223 11352 : CALL dbt_clear(work_dbt_3c_1)
1224 11352 : CALL dbt_clear(work_dbt_3c_2)
1225 11352 : CALL dbt_clear(work_dbt_2c)
1226 11352 : CALL timestop(handle)
1227 :
1228 11352 : END SUBROUTINE get_sigma_noenv
1229 : ! **************************************************************************************************
1230 : !> \brief Creates the RI matrix and populates it with correct values
1231 : !> \note Tensor contains Hartree elements in the auxiliary basis
1232 : !> \param qs_env Quickstep environment - entry point of calculation
1233 : !> \author Stepan Marek
1234 : !> \date 01.2024
1235 : ! **************************************************************************************************
1236 50 : SUBROUTINE init_hartree(rtbse_env, v_dbcsr)
1237 : TYPE(rtbse_env_type), POINTER, INTENT(IN) :: rtbse_env
1238 : TYPE(dbcsr_type) :: v_dbcsr
1239 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1240 : TYPE(libint_potential_type) :: coulomb_op
1241 : TYPE(cp_fm_type) :: V_fm
1242 : TYPE(cp_fm_type) :: metric_fm
1243 : TYPE(cp_fm_type) :: metric_inv_fm, &
1244 : work_fm
1245 : TYPE(dbcsr_type), DIMENSION(:), ALLOCATABLE :: V_dbcsr_a, &
1246 50 : metric_dbcsr
1247 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1248 50 : POINTER :: nl_2c
1249 :
1250 50 : bs_env => rtbse_env%bs_env
1251 :
1252 : ! Allocate for bare Hartree term
1253 100 : ALLOCATE (V_dbcsr_a(1))
1254 100 : ALLOCATE (metric_dbcsr(1))
1255 50 : CALL dbcsr_create(V_dbcsr_a(1), name="Hartree_dbcsr", template=bs_env%mat_RI_RI%matrix)
1256 50 : CALL dbcsr_create(metric_dbcsr(1), name="RI_metric_dbcsr", template=bs_env%mat_RI_RI%matrix)
1257 :
1258 : ! Calculate full coulomb RI basis elements - V _ (PQ) matrix
1259 50 : NULLIFY (nl_2c)
1260 : CALL build_2c_neighbor_lists(nl_2c, bs_env%basis_set_RI, bs_env%basis_set_RI, &
1261 : coulomb_op, "Coulomb_neighbor_2c_list", rtbse_env%qs_env, &
1262 50 : sym_ij=.FALSE., molecular=.TRUE.)
1263 : CALL build_2c_integrals(V_dbcsr_a, bs_env%eps_filter, rtbse_env%qs_env, nl_2c, &
1264 : bs_env%basis_set_RI, bs_env%basis_set_RI, coulomb_op, &
1265 50 : do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI)
1266 : ! Calculate the RI metric elements
1267 : ! nl_2c is automatically rewritten (even reallocated) in this routine
1268 : CALL build_2c_neighbor_lists(nl_2c, bs_env%basis_set_RI, bs_env%basis_set_RI, &
1269 : bs_env%ri_metric, "Metric_neighbor_2c_list", rtbse_env%qs_env, &
1270 50 : sym_ij=.FALSE., molecular=.TRUE.)
1271 : CALL build_2c_integrals(metric_dbcsr, bs_env%eps_filter, rtbse_env%qs_env, nl_2c, &
1272 : bs_env%basis_set_RI, bs_env%basis_set_RI, bs_env%ri_metric, &
1273 50 : do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI)
1274 : ! nl_2c no longer needed
1275 50 : CALL release_neighbor_list_sets(nl_2c)
1276 50 : CALL cp_fm_create(metric_fm, bs_env%fm_RI_RI%matrix_struct)
1277 50 : CALL cp_fm_set_all(metric_fm, 0.0_dp)
1278 50 : CALL cp_fm_create(metric_inv_fm, bs_env%fm_RI_RI%matrix_struct)
1279 50 : CALL cp_fm_set_all(metric_inv_fm, 0.0_dp)
1280 50 : CALL cp_fm_create(work_fm, bs_env%fm_RI_RI%matrix_struct)
1281 50 : CALL cp_fm_set_all(work_fm, 0.0_dp)
1282 50 : CALL copy_dbcsr_to_fm(metric_dbcsr(1), metric_fm)
1283 50 : CALL cp_fm_invert(metric_fm, metric_inv_fm)
1284 50 : CALL cp_fm_create(V_fm, bs_env%fm_RI_RI%matrix_struct)
1285 50 : CALL cp_fm_set_all(V_fm, 0.0_dp)
1286 : ! Multiply by the inverse from each side (M^-1 is symmetric)
1287 : CALL cp_dbcsr_sm_fm_multiply(V_dbcsr_a(1), metric_inv_fm, &
1288 50 : work_fm, bs_env%n_RI)
1289 : CALL parallel_gemm("N", "N", bs_env%n_RI, bs_env%n_RI, bs_env%n_RI, &
1290 50 : 1.0_dp, metric_inv_fm, work_fm, 0.0_dp, V_fm)
1291 : ! Now, create the tensor from the matrix
1292 : ! First, convert full matrix to dbcsr
1293 50 : CALL dbcsr_clear(V_dbcsr_a(1))
1294 50 : CALL copy_fm_to_dbcsr(V_fm, V_dbcsr_a(1))
1295 50 : CALL dbcsr_create(v_dbcsr, "Hartree ri", V_dbcsr_a(1))
1296 50 : CALL dbcsr_copy(v_dbcsr, V_dbcsr_a(1))
1297 : ! Create and copy distinctly, so that unnecessary objects can be destroyed
1298 : ! Destroy all unnecessary matrices
1299 50 : CALL dbcsr_release(V_dbcsr_a(1))
1300 50 : CALL dbcsr_release(metric_dbcsr(1))
1301 50 : DEALLOCATE (V_dbcsr_a)
1302 50 : DEALLOCATE (metric_dbcsr)
1303 50 : CALL cp_fm_release(V_fm)
1304 : ! CALL cp_fm_release(metric_fm(1,1))
1305 50 : CALL cp_fm_release(metric_fm)
1306 : ! DEALLOCATE(metric_fm)
1307 50 : CALL cp_fm_release(work_fm)
1308 50 : CALL cp_fm_release(metric_inv_fm)
1309 300 : END SUBROUTINE init_hartree
1310 : ! **************************************************************************************************
1311 : !> \brief Calculates the Hartree matrix in the atomic orbital basis, given a density matrix, in local arrays
1312 : !> Calculates the values for single spin species present in given rho
1313 : !> \param qs_env Entry point
1314 : !> \param rtbse_env Entry point of GWBSE - uses rho_dbcsr and some complex_workspace
1315 : !> \param rho_ao Density matrix in ao basis
1316 : !> \param v_ao Overwritten by the Hartree matrix in the atomic orbital basis
1317 : !> \author Stepan Marek
1318 : !> \date 01.2025
1319 : ! **************************************************************************************************
1320 5520 : SUBROUTINE get_hartree_env(rtbse_env, rho_fm, v_fm)
1321 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1322 : TYPE(cp_cfm_type) :: rho_fm
1323 : TYPE(cp_fm_type) :: v_fm
1324 : TYPE(mp_para_env_type), POINTER :: para_env
1325 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1326 :
1327 5520 : CALL get_qs_env(rtbse_env%qs_env, para_env=para_env, bs_env=bs_env)
1328 :
1329 : CALL get_hartree_noenv(v_fm, rho_fm, rtbse_env%int_3c_array, rtbse_env%v_dbcsr, &
1330 : rtbse_env%n_RI, bs_env%sizes_RI, &
1331 5520 : para_env, rtbse_env%rho_dbcsr, rtbse_env%v_ao_dbcsr)
1332 5520 : END SUBROUTINE get_hartree_env
1333 : ! **************************************************************************************************
1334 : !> \brief Calculates the Hartree matrix in the atomic orbital basis, given a density matrix, in local arrays
1335 : !> Calculates the values for single spin species present in given rho
1336 : !> \param v_fm Hartree potential in atomic orbital basis - is overwritten by the updated potential
1337 : !> \param rho_fm Density matrix corresponding to single spin species, in atomic orbital basis
1338 : !> \param int_3c Previously allocated array (best to use create_hartree_ri_3c) for 3c integrals
1339 : !> \param v_dbcsr Previously calculated 2c Coulomb repulsion between RI orbitals
1340 : !> \param n_RI Number of RI basis orbitals
1341 : !> \param sizes_RI Number of RI basis orbitals per atom
1342 : !> \param para_env MPI Parallel environment (used for summation across ranks)
1343 : !> \param rho_dbcsr Previously created dbcsr matrix, used as workspace
1344 : !> \param v_ao_dbcsr Previously created dbcsr matrix, used as workspace
1345 : !> \author Stepan Marek
1346 : !> \date 01.2025
1347 : ! **************************************************************************************************
1348 5520 : SUBROUTINE get_hartree_noenv(v_fm, rho_fm, int_3c, v_dbcsr, n_RI, sizes_RI, para_env, rho_dbcsr, v_ao_dbcsr)
1349 : TYPE(cp_fm_type) :: v_fm
1350 : TYPE(cp_cfm_type), INTENT(IN) :: rho_fm
1351 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c
1352 : TYPE(dbcsr_type) :: v_dbcsr
1353 : INTEGER :: n_RI
1354 : INTEGER, DIMENSION(:) :: sizes_RI
1355 : TYPE(mp_para_env_type), POINTER :: para_env
1356 : TYPE(dbcsr_type) :: rho_dbcsr, v_ao_dbcsr
1357 : CHARACTER(len=*), PARAMETER :: routineN = "get_hartree"
1358 : TYPE(dbcsr_iterator_type) :: iterator_matrix
1359 : INTEGER :: i, j, k, n, nblocks, ind_1, ind_2, row_offset, col_offset, &
1360 : row_size, col_size, j_n_AO, k_n_AO, i_n_RI, &
1361 : ri_offset, ind_i, handle
1362 5520 : REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: Pvector, Qvector
1363 5520 : REAL(kind=dp), DIMENSION(:, :), POINTER :: block_matrix
1364 : INTEGER :: nblkrows_local, nblkcols_local, j_blk, k_blk, j_offset, k_offset
1365 5520 : INTEGER, DIMENSION(:), POINTER :: local_blk_rows, local_blk_cols
1366 : LOGICAL :: found
1367 :
1368 : MARK_USED(i_n_RI)
1369 : MARK_USED(ri_offset)
1370 : MARK_USED(ind_i)
1371 :
1372 : ! No memory optimisation so far - calculate all 3cs an all ranks
1373 : ! Importantly - dbcsr blocks are ordered by atoms - i.e. ethene with 6 atoms will have 6x6 block structure
1374 : ! Number of basis states on each basis set is known is post_scf_bandstructure env
1375 :
1376 5520 : CALL timeset(routineN, handle)
1377 :
1378 : ! Allocate the Q and Pvector on each rank
1379 16560 : ALLOCATE (Qvector(n_RI), source=0.0_dp)
1380 11040 : ALLOCATE (Pvector(n_RI), source=0.0_dp)
1381 :
1382 : ! First step - analyze the structure of copied dbcsr matrix on all ranks
1383 5520 : CALL dbcsr_clear(rho_dbcsr)
1384 : ! Only the real part of the density matrix contributes
1385 : ! Use v_fm as workspace
1386 5520 : CALL cp_cfm_to_fm(msource=rho_fm, mtargetr=v_fm)
1387 5520 : CALL copy_fm_to_dbcsr(v_fm, rho_dbcsr)
1388 5520 : j_offset = 0
1389 : CALL dbcsr_get_info(rho_dbcsr, nblkrows_local=nblkrows_local, nblkcols_local=nblkcols_local, &
1390 5520 : local_rows=local_blk_rows, local_cols=local_blk_cols)
1391 11040 : DO j_blk = 1, nblkrows_local
1392 5520 : k_offset = 0
1393 16560 : DO k_blk = 1, nblkcols_local
1394 : ! Check whether we can retrieve the rho block
1395 : ! TODO : Handle transposed case?
1396 : CALL dbcsr_get_block_p(rho_dbcsr, local_blk_rows(j_blk), local_blk_cols(k_blk), &
1397 11040 : block=block_matrix, found=found, row_size=row_size, col_size=col_size)
1398 : ! If the block is not found, then the density matrix here has below threshold values
1399 11040 : IF (.NOT. found) CYCLE
1400 : ! With the block retrieved, add its contributions to the Q-vector
1401 : !$OMP PARALLEL DO DEFAULT(none) PRIVATE(i,j,k) &
1402 11040 : !$OMP SHARED(n_RI, row_size, col_size, Qvector, int_3c, j_offset, k_offset, block_matrix)
1403 : DO i = 1, n_RI
1404 : DO j = 1, row_size
1405 : DO k = 1, col_size
1406 : Qvector(i) = Qvector(i) + int_3c(j_offset + j, k_offset + k, i)*block_matrix(j, k)
1407 : END DO
1408 : END DO
1409 : END DO
1410 : !$OMP END PARALLEL DO
1411 : ! Increment k-offset - setup for the next block
1412 27600 : k_offset = k_offset + col_size
1413 : END DO
1414 : ! Increments the j-offset - row_size is carried over from the last iteration
1415 11040 : j_offset = j_offset + row_size
1416 : END DO
1417 : ! Now, each rank has contributions from D_jk within its scope
1418 : ! Need to sum over different ranks to get the total vector on all ranks
1419 5520 : CALL para_env%sum(Qvector)
1420 : ! Once this is done, Pvector is current on all ranks
1421 : ! Continue with V_PQ summation
1422 5520 : nblocks = dbcsr_get_num_blocks(v_dbcsr)
1423 5520 : CALL dbcsr_iterator_start(iterator_matrix, v_dbcsr)
1424 16560 : DO n = 1, nblocks
1425 : ! TODO : Try OMP parallelisation over different blocks - expect many more available speedup for large systems
1426 : CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, &
1427 11040 : row_offset=row_offset, col_offset=col_offset, row_size=row_size, col_size=col_size)
1428 : ! TODO : Better names for RI
1429 11040 : j_n_AO = sizes_RI(ind_1)
1430 11040 : k_n_AO = sizes_RI(ind_2)
1431 : ! The allocations are as follows
1432 : !$OMP PARALLEL DO DEFAULT(none) PRIVATE(j,k) &
1433 16560 : !$OMP SHARED(block_matrix, Pvector, Qvector,j_n_AO,k_n_AO,row_offset,col_offset)
1434 : DO j = 1, j_n_AO
1435 : DO k = 1, k_n_AO
1436 : Pvector(j + row_offset - 1) = Pvector(j + row_offset - 1) + block_matrix(j, k)*Qvector(k + col_offset - 1)
1437 : END DO
1438 : END DO
1439 : !$OMP END PARALLEL DO
1440 : END DO
1441 5520 : CALL dbcsr_iterator_stop(iterator_matrix)
1442 : ! Again, make sure that the P vector is present on all ranks
1443 5520 : CALL para_env%sum(Pvector)
1444 : ! Now, for the final trick, iterate over local blocks of v_ao_dbcsr to get the Hartree as dbcsr, then convert to fm
1445 : ! TODO : Clear or set blocks to zero
1446 : ! CALL dbcsr_clear(v_ao_dbcsr)
1447 5520 : j_offset = 0
1448 11040 : DO j_blk = 1, nblkrows_local
1449 5520 : k_offset = 0
1450 16560 : DO k_blk = 1, nblkcols_local
1451 : ! Check whether we can retrieve the rho block
1452 : ! TODO : Handle transposed case?
1453 : CALL dbcsr_get_block_p(v_ao_dbcsr, local_blk_rows(j_blk), local_blk_cols(k_blk), &
1454 11040 : block=block_matrix, found=found, row_size=row_size, col_size=col_size)
1455 : ! If the block is not found, reserve it
1456 11040 : IF (.NOT. found) THEN
1457 : ! Reservations
1458 92 : CALL dbcsr_reserve_blocks(v_ao_dbcsr, local_blk_rows(j_blk:j_blk), local_blk_cols(k_blk:k_blk))
1459 : ! Rerun the getter to get the new block
1460 : CALL dbcsr_get_block_p(v_ao_dbcsr, local_blk_rows(j_blk), local_blk_cols(k_blk), &
1461 92 : block=block_matrix, found=found, row_size=row_size, col_size=col_size)
1462 : END IF
1463 : ! With the block retrieved, contract with the P vector
1464 : !$OMP PARALLEL DO DEFAULT(none) PRIVATE(i,j,k) &
1465 11040 : !$OMP SHARED(row_size, col_size, n_RI, block_matrix, Pvector, int_3c, j_offset, k_offset)
1466 : DO j = 1, row_size
1467 : DO k = 1, col_size
1468 : block_matrix(j, k) = 0.0_dp
1469 : DO i = 1, n_RI
1470 : block_matrix(j, k) = block_matrix(j, k) + Pvector(i)*int_3c(j_offset + j, k_offset + k, i)
1471 : END DO
1472 : END DO
1473 : END DO
1474 : !$OMP END PARALLEL DO
1475 : ! Increment k-offset - setup for the next block
1476 27600 : k_offset = k_offset + col_size
1477 : END DO
1478 : ! Increments the j-offset - row_size is carried over from the last iteration
1479 11040 : j_offset = j_offset + row_size
1480 : END DO
1481 : ! Since P vector was present on all the ranks, v_dbcsr_ao has the complete Hartree result
1482 : ! copy_dbcsr_to_fm should set all values in v_fm to zero
1483 5520 : CALL copy_dbcsr_to_fm(v_ao_dbcsr, v_fm)
1484 5520 : DEALLOCATE (Qvector)
1485 5520 : DEALLOCATE (Pvector)
1486 :
1487 5520 : CALL timestop(handle)
1488 22080 : END SUBROUTINE get_hartree_noenv
1489 : ! **************************************************************************************************
1490 : !> \brief Calculates the exponential of a matrix in a generalized eigenvalue problem. Specifically,
1491 : !> it assumes we have a Hermitian matrix A in the eigenvalue problem AX = BXE, where B is some overlap
1492 : !> matrix and E is a diagonal matrix of real eigenvalues. Then, it calculates
1493 : !> exp(B^(-1) A) = X exp(E) X^C B
1494 : !> \param amatrix Matrix to exponentiate
1495 : !> \param bmatrix Overlap matrix
1496 : !> \param exponential Exponential exp(B^(-1) A) is stored here after the routine is finished
1497 : !> \param eig_scale_opt Optionally scale eigenvalues by a complex number before exponentiating them
1498 : !> \param work_opt Optionally provide workspace (of size at least 4) that is used in the calculation
1499 : !> \author Stepan Marek
1500 : !> \date 09.2024
1501 : ! **************************************************************************************************
1502 0 : SUBROUTINE cp_cfm_gexp(amatrix, bmatrix, exponential, eig_scale_opt, work_opt)
1503 : ! TODO : Do interface for real matrices
1504 : TYPE(cp_cfm_type), INTENT(IN) :: amatrix
1505 : TYPE(cp_cfm_type), INTENT(IN) :: bmatrix
1506 : TYPE(cp_cfm_type) :: exponential
1507 : COMPLEX(kind=dp), INTENT(IN), OPTIONAL :: eig_scale_opt
1508 : TYPE(cp_cfm_type), DIMENSION(:), POINTER, OPTIONAL :: work_opt
1509 : CHARACTER(len=*), PARAMETER :: routineN = "cp_cfm_gexp"
1510 : COMPLEX(kind=dp) :: eig_scale
1511 0 : REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues
1512 0 : COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: expvalues
1513 0 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: work
1514 : LOGICAL :: deallocate_work
1515 : INTEGER :: nrow, i, handle
1516 :
1517 0 : CALL timeset(routineN, handle)
1518 :
1519 : ! Argument parsing and sanity checks
1520 0 : IF (PRESENT(eig_scale_opt)) THEN
1521 0 : eig_scale = eig_scale_opt
1522 : ELSE
1523 : eig_scale = CMPLX(1.0, 0.0, kind=dp)
1524 : END IF
1525 :
1526 0 : NULLIFY (work)
1527 0 : deallocate_work = .TRUE.
1528 0 : IF (PRESENT(work_opt)) THEN
1529 0 : deallocate_work = SIZE(work_opt) < 4
1530 : END IF
1531 0 : IF (.NOT. deallocate_work) THEN
1532 0 : work => work_opt
1533 : ELSE
1534 0 : ALLOCATE (work(4))
1535 : ! Allocate the work storage on the fly
1536 0 : DO i = 1, 4
1537 0 : CALL cp_cfm_create(work(i), amatrix%matrix_struct)
1538 : END DO
1539 : END IF
1540 :
1541 0 : nrow = amatrix%matrix_struct%nrow_global
1542 :
1543 0 : ALLOCATE (eigenvalues(nrow))
1544 0 : ALLOCATE (expvalues(nrow))
1545 :
1546 : ! Do not change the amatrix and bmatrix - need to copy them first
1547 0 : CALL cp_cfm_to_cfm(amatrix, work(1))
1548 0 : CALL cp_cfm_to_cfm(bmatrix, work(2))
1549 :
1550 : ! Solve the generalized eigenvalue equation
1551 0 : CALL cp_cfm_geeig(work(1), work(2), work(3), eigenvalues, work(4))
1552 :
1553 : ! Scale and exponentiate the eigenvalues
1554 0 : expvalues(:) = EXP(eigenvalues(:)*eig_scale)
1555 :
1556 : ! Copy eigenvectors to column scale them
1557 0 : CALL cp_cfm_to_cfm(work(3), work(1))
1558 : ! X * exp(E)
1559 0 : CALL cp_cfm_column_scale(work(1), expvalues)
1560 :
1561 : ! Carry out the remaining operations
1562 : ! X * exp(E) * X^C
1563 : CALL parallel_gemm("N", "C", nrow, nrow, nrow, &
1564 : CMPLX(1.0, 0.0, kind=dp), work(1), work(3), &
1565 0 : CMPLX(0.0, 0.0, kind=dp), work(2))
1566 : ! X * exp(E) * X^C * B
1567 : CALL parallel_gemm("N", "N", nrow, nrow, nrow, &
1568 : CMPLX(1.0, 0.0, kind=dp), work(2), bmatrix, &
1569 0 : CMPLX(0.0, 0.0, kind=dp), exponential)
1570 :
1571 : ! Deallocate work storage if necessary
1572 0 : IF (deallocate_work) THEN
1573 0 : DO i = 1, 4
1574 0 : CALL cp_cfm_release(work(i))
1575 : END DO
1576 0 : DEALLOCATE (work)
1577 : END IF
1578 :
1579 0 : DEALLOCATE (eigenvalues)
1580 0 : DEALLOCATE (expvalues)
1581 :
1582 0 : CALL timestop(handle)
1583 0 : END SUBROUTINE cp_cfm_gexp
1584 : END MODULE rt_bse
|
