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