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