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