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