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