diff --git a/docs/methods/properties/optical/index.md b/docs/methods/properties/optical/index.md index 3aed11b1af..657a4cbf13 100644 --- a/docs/methods/properties/optical/index.md +++ b/docs/methods/properties/optical/index.md @@ -7,6 +7,7 @@ maxdepth: 1 --- tddft bethe-salpeter +rtbse ``` Optical spectroscopy is a technique used to study the interaction between light and matter. It diff --git a/docs/methods/properties/optical/rtbse.md b/docs/methods/properties/optical/rtbse.md new file mode 100644 index 0000000000..e12c30b9c6 --- /dev/null +++ b/docs/methods/properties/optical/rtbse.md @@ -0,0 +1,171 @@ +# Real Time Bethe-Salpeter Propagation + +Instead of using TDDFT functionals, an approximation to the self-energy is employed to calculate the +time dependent behaviour of the density matrix +\[[Attaccalite2011](http://dx.doi.org/10.1103/PhysRevB.84.245110)\]. This requires a previous +determination of the screened Coulomb potential, done via the bandstructure +[GW](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.BANDSTRUCTURE.GW) calculation. + +The single particle density matrix $\hat{\rho}$ is propagated following the equation of motion + +$$ \frac{\mathrm{d} \hat{\rho}}{\mathrm{d} t} = -i [\hat{H}(t), \hat{\rho}(t)] $$ + +This equation is solved by steps + +$$ \hat{\rho} (t + \Delta t) = \mathrm{e} ^ {- i \hat{H} (t+\Delta t) \Delta t/2} \mathrm{e} ^ {-i \hat{H}(t) \Delta t/2} +\hat{\rho} (t) \mathrm{e} ^ {i \hat{H}(t) \Delta t/2} \mathrm{e} ^ {i \hat{H} (t + \Delta t) \Delta t/2}$$ + +which is called the _enforced time reversal +scheme_\[[Castro2004](https://doi.org/10.1063/1.1774980)\]. The effective Hamiltonian is given as + +$$ \hat{H}(t) = \hat{h}^{G0W0} + \hat{U} (t) + +\hat{V}^{\mathrm{Hartree}} [\hat{\rho}(t)] - \hat{V}^{\mathrm{Hartree}} [\hat{\rho}_0] + +\hat{\Sigma}^{\mathrm{COHSEX}}[\hat{\rho}(t)] - \hat{\Sigma}^{\mathrm{COHSEX}}[\hat{\rho}_0] +$$ + +where $\hat{\rho}_0$ is the density matrix determined from the molecular orbitals used in +[GW](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.BANDSTRUCTURE.GW) and $\hat{U}(t)$ is the external applied +field. + +## Excitation scheme + +Without the external field $\hat{U}(t)$, the density matrix only rotates in phase but does not +produce any measurable dynamics. The excitation of the dynamics can be done either by a real time +pulse (i.e. at each point, $\hat{U}(t)$ follows form due to some finite time dependent field +$\vec{E}(t)$) or by an infinitely sharp delta pulse, which we can understand as the limit of +$\vec{E}(t) \to I \vec{e} \delta(t)$, where $I$ is the delta pulse intensity and $\vec{e}$ its +direction. + +## Observables + +The dynamics can be traced through time with electric dipole moment associated with the density +matrix + +$$ \mu_i(t) = \mathrm{Tr} (\hat{\rho}(t) \hat{x}_i) +$$ + +The electric polarizability (which is related to the photon absorption spectrum) is then determined +as + +$$ \alpha_{ij} (\omega) = \frac{\mu_i(\omega)}{E_j(\omega)} +$$ + +where we Fourier transformed to the frequency domain. In order to model the Fourier transform of +infinitely oscillating dipole moments, we introduce a damping factor +$\gamma$\[[Müller2020](https://doi.org/10.1002/jcc.26412)\] + +$$ \mu_i(\omega) = \int _ 0 ^ T \mathrm{d}t \mathrm{e}^{-\gamma t} \mathrm{e} ^ {i \omega t} \mu_i(t) +$$ + +One can easily verify that for real FT of the applied field, this leads Lorentzian peaks at the +frequencies of the oscillations of the moments present in the imaginary part of the corresponding +polarizability element. + +## Running the Propagation + +To run the RTBSE propagation, include the +[REAL_TIME_PROPAGATION](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION) section in the input file +and set the [RTP_METHOD](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.RTP_METHOD) to `RTBSE` - +otherwise, the standard TDDFT propagation is employed. + +### ETRS Precision + +The precision of the self-consistency in the ETRS loop is controlled by the +[EPS_ITER](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.EPS_ITER) keyword. Smaller threshold +(larger precision) lead to more stable propagation, but might require smaller timestep/more +self-consistent iterations. + +[MAX_ITER](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.MAX_ITER) keyword is used to determine +the maximum number of self-consistent iterations for a single time step before the cycle is broken +and non-convergence is reported. + +### Exponential Method + +The method used for the exponentiation of the Hamiltonian is set in the +[MAT_EXP](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.MAT_EXP) keyword, with the following +methods implemented for both TDDFT and RTBSE + +- `BCH` - calculates the effect of matrix exponential by series of commutators using + Baker-Campbell-Hausdorff expansion + +and the following methods implemented only for RTBSE + +- `EXACT` - Diagonalizes the instantaneous Hamiltonian to determine the exponential exactly + +For inexact methods, a threshold for the cutoff of exponential series is provided by the +[EXP_ACCURACY](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.EXP_ACCURACY) keyword. For these, +the [MAX_ITER](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.MAX_ITER) keyword also sets the +maximum number of iterations before the program is stopped and non-convergence is reported. + +Use [RTP_METHOD](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.RTP_METHOD) to start the +calculation by setting it to TDAGW. + +### Excitation Method + +The real time pulse can be specified in the [EFIELD](#CP2K_INPUT.FORCE_EVAL.DFT.EFIELD) section. + +If delta pulse is required instead, use +[APPLY_DELTA_PULSE](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.APPLY_DELTA_PULSE) with +[DELTA_PULSE_DIRECTION](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.DELTA_PULSE_DIRECTION) used +for defining the $\vec{e}$ vector and +[DELTA_PULSE_SCALE](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.DELTA_PULSE_SCALE) setting the +$I$ scale of the delta pulse (in atomic units). Note that the definition of the vector is different +from the definition used in the TDDFT method. + +The actual value of $I \vec{e}$ is printed out in atomic units. + +### Printing observables + +The code is so far optimised for printing the polarizability elements, which are linked to the +absorption spectrum. The printing of all available properties is controlled in the +[PRINT](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT) section - the ones relevant for +RTBSE propagation are listed here + +- [DENSITY_MATRIX](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.DENSITY_MATRIX) - Prints + the elements of the density matrix in the MO basis into a file at every timestep +- [FIELD](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.FIELD) - Prints the elements of the + electric field applied at every time step +- [MOMENTS](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.MOMENTS) - Prints the electric + dipole moment elements reported at every time step +- [MOMENTS_FT](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.MOMENTS_FT) - Prints the + Fourier transform of the dipole moment elements time series +- [POLARIZABILITY](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.POLARIZABILITY) - Prints + an element of the Fourier transform of polarizability. +- [RESTART](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.RESTART) - Controls the name of + the restart file + +When [RESTART](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.RESTART), +[MOMENTS](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.MOMENTS) and +[FIELD](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PRINT.FIELD) are saved into files, one can +continue running the calculation in the same directory for longer time without rerunning the already +calculated time steps. Note that total length of the propagation time controls the energy/frequency +precision, while timestep size controls the energy/frequency range. + +### Example Input + +A typical input file which runs the RTBSE propagation will have the +[REAL_TIME_PROPAGATION](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION) section similar to this +one + +``` +&REAL_TIME_PROPAGATION + RTP_METHOD RTBSE + ! ETRS self-consistent threshold + EPS_ITER 1.0E-8 + MAT_EXP BCH + EXP_ACCURACY 1.0E-12 + INITIAL_WFN RT_RESTART + APPLY_DELTA_PULSE + DELTA_PULSE_SCALE 0.01 + DELTA_PULSE_DIRECTION 1 0 0 + &PRINT + ! No need to print the field for delta pulse + &MOMENTS + FILENAME MOMENTS + &END MOMENTS + &POLARIZABILITY + FILENAME POLARIZABILITY + &END POLARIZABILITY + &END PRINT +&END REAL_TIME_PROPAGATION +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e3803ef295..f6cfc73875 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -984,7 +984,10 @@ list( emd/rt_propagation_methods.F emd/rt_propagation_output.F emd/rt_propagation_utils.F - emd/rt_propagator_init.F) + emd/rt_propagator_init.F + emd/rt_bse.F + emd/rt_bse_io.F + emd/rt_bse_types.F) list( APPEND diff --git a/src/cp_control_types.F b/src/cp_control_types.F index 7953a6ca28..6b245016be 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -12,7 +12,9 @@ MODULE cp_control_types USE cp_fm_types, ONLY: cp_fm_release,& cp_fm_type - USE input_constants, ONLY: do_full_density + USE input_constants, ONLY: do_full_density,& + rtp_bse_ham_G0W0,& + rtp_method_tddft USE kinds, ONLY: default_path_length,& default_string_length,& dp @@ -69,6 +71,8 @@ MODULE cp_control_types INTEGER :: mat_exp = 0 INTEGER :: propagator = 0 LOGICAL :: fixed_ions = .FALSE. + INTEGER :: rtp_method = rtp_method_tddft + INTEGER :: rtbse_ham = rtp_bse_ham_G0W0 INTEGER :: initial_wfn = 0 REAL(dp) :: eps_exp = 0.0_dp LOGICAL :: initial_step = .FALSE. diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index dc88438662..8c2b9f0e3b 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -2425,6 +2425,10 @@ SUBROUTINE read_rtp_section(dft_control, rtp_section) i_val=dft_control%rtp_control%aspc_order) CALL section_vals_val_get(rtp_section, "EXP_ACCURACY", & r_val=dft_control%rtp_control%eps_exp) + CALL section_vals_val_get(rtp_section, "RTP_METHOD", & + i_val=dft_control%rtp_control%rtp_method) + CALL section_vals_val_get(rtp_section, "RTBSE_HAMILTONIAN", & + i_val=dft_control%rtp_control%rtbse_ham) CALL section_vals_val_get(rtp_section, "PROPAGATOR", & i_val=dft_control%rtp_control%propagator) CALL section_vals_val_get(rtp_section, "EPS_ITER", & diff --git a/src/emd/PACKAGE b/src/emd/PACKAGE index e9af768351..f7cf8c6f30 100644 --- a/src/emd/PACKAGE +++ b/src/emd/PACKAGE @@ -13,5 +13,6 @@ "../pw_env", "../arnoldi", "../dbx", + "../dbt", ], } diff --git a/src/emd/rt_bse.F b/src/emd/rt_bse.F new file mode 100644 index 0000000000..a446f76d4f --- /dev/null +++ b/src/emd/rt_bse.F @@ -0,0 +1,1416 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Routines for the propagation via RT-BSE method. +!> \note The control is handed directly from cp2k_runs +!> \author Stepan Marek (12.23) +! ************************************************************************************************** + +MODULE rt_bse + USE cp_control_types, ONLY: dft_control_type, & + rtp_control_type + USE qs_environment_types, ONLY: get_qs_env, & + qs_environment_type + USE force_env_types, ONLY: force_env_type + USE qs_mo_types, ONLY: mo_set_type, & + init_mo_set + USE atomic_kind_types, ONLY: atomic_kind_type + USE qs_kind_types, ONLY: qs_kind_type + USE particle_types, ONLY: particle_type + USE rt_propagation_types, ONLY: get_rtp, & + rt_prop_type + USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type + USE cp_fm_types, ONLY: cp_fm_type, & + cp_fm_p_type, & + cp_fm_to_fm, & + cp_fm_create, & + cp_fm_set_all, & + cp_fm_release, & + cp_fm_get_element, & + cp_fm_set_element, & + cp_fm_get_info, & + cp_fm_write_unformatted, & + cp_fm_read_unformatted, & + cp_fm_write_formatted + USE cp_cfm_types, ONLY: cp_cfm_type, & + cp_cfm_p_type, & + cp_fm_to_cfm, & + cp_cfm_to_cfm, & + cp_cfm_to_fm, & + cp_cfm_create, & + cp_cfm_get_info, & + cp_cfm_release + USE kinds, ONLY: dp, & + default_path_length + USE cp_dbcsr_api, ONLY: dbcsr_p_type, & + dbcsr_type, & + dbcsr_print, & + dbcsr_has_symmetry, & + dbcsr_desymmetrize, & + dbcsr_create, & + dbcsr_release, & + dbcsr_copy, & + dbcsr_scale, & + dbcsr_add, & + dbcsr_set, & + dbcsr_clear, & + dbcsr_setname, & + dbcsr_iterator_type, & + dbcsr_iterator_start, & + dbcsr_iterator_stop, & + dbcsr_iterator_blocks_left, & + dbcsr_iterator_next_block, & + dbcsr_put_block, & + dbcsr_reserve_blocks, & + dbcsr_get_num_blocks, & + dbcsr_get_block_p, & + dbcsr_get_info + USE OMP_LIB, ONLY: omp_get_thread_num, & + omp_get_num_threads, & + omp_set_num_threads, & + omp_get_max_threads + USE dbt_api, ONLY: dbt_create, & + dbt_clear, & + dbt_contract, & + dbt_copy_matrix_to_tensor, & + dbt_copy_tensor_to_matrix, & + dbt_copy, & + dbt_destroy, & + dbt_type, & + dbt_pgrid_type, & + dbt_pgrid_create, & + dbt_pgrid_destroy, & + dbt_mp_environ_pgrid, & + dbt_default_distvec, & + dbt_distribution_type, & + dbt_distribution_new, & + dbt_distribution_destroy, & + dbt_iterator_type, & + dbt_iterator_start, & + dbt_iterator_stop, & + dbt_iterator_blocks_left, & + dbt_iterator_next_block, & + dbt_put_block, & + dbt_get_block, & + dbt_reserve_blocks, & + dbt_get_num_blocks, & + dbt_get_info + USE libint_2c_3c, ONLY: libint_potential_type + USE mp2_ri_2c, ONLY: RI_2c_integral_mat + USE qs_tensors, ONLY: neighbor_list_3c_destroy, & + build_2c_integrals, & + build_2c_neighbor_lists + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type, & + release_neighbor_list_sets + USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm, & + copy_fm_to_dbcsr, & + dbcsr_allocate_matrix_set, & + dbcsr_deallocate_matrix_set, & + cp_dbcsr_sm_fm_multiply, & + copy_cfm_to_dbcsr, & + copy_dbcsr_to_cfm + USE cp_fm_basic_linalg, ONLY: cp_fm_scale, & + cp_fm_invert, & + cp_fm_trace, & + cp_fm_transpose, & + cp_fm_norm, & + cp_fm_column_scale, & + cp_fm_scale_and_add + USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add, & + cp_cfm_scale, & + cp_cfm_transpose, & + cp_cfm_norm, & + cp_cfm_trace, & + cp_cfm_column_scale + USE cp_fm_diag, ONLY: cp_fm_geeig + USE cp_cfm_diag, ONLY: cp_cfm_geeig + USE parallel_gemm_api, ONLY: parallel_gemm + USE qs_moments, ONLY: build_local_moment_matrix, & + build_berry_moment_matrix + USE moments_utils, ONLY: get_reference_point + USE qs_mo_io, ONLY: read_mo_set_from_restart + USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix + USE force_env_methods, ONLY: force_env_calc_energy_force + USE qs_energy_init, ONLY: qs_energies_init + USE qs_energy_utils, ONLY: qs_energies_properties + USE efield_utils, ONLY: make_field + USE rt_propagator_init, ONLY: rt_initialize_rho_from_mos + USE rt_propagation_methods, ONLY: s_matrices_create + USE gw_large_cell_Gamma, ONLY: compute_3c_integrals + USE gw_integrals, ONLY: build_3c_integral_block + USE matrix_exp, ONLY: taylor_full_complex + USE message_passing, ONLY: mp_para_env_type + USE input_section_types, ONLY: section_vals_get_subs_vals, & + section_vals_type + USE cell_types, ONLY: cell_type + USE input_constants, ONLY: rtp_bse_ham_ks, & + rtp_bse_ham_g0w0, & + use_mom_ref_coac, & + use_mom_ref_zero, & + do_taylor, & + do_bch, & + do_exact, & + use_rt_restart + USE rt_bse_types, ONLY: rtbse_env_type, & + create_rtbse_env, & + release_rtbse_env, & + multiply_cfm_fm, & + multiply_fm_cfm + USE rt_bse_io, ONLY: output_moments, & + read_moments, & + read_field, & + output_moments_ft, & + output_polarizability, & + output_field, & + output_mos_contravariant, & + read_restart, & + output_restart, & + print_timestep_info, & + print_etrs_info_header, & + print_etrs_info, & + print_rtbse_header_info + USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_integrate_potential_product, ONLY: integrate_v_rspace + USE qs_collocate_density, ONLY: calculate_rho_elec + USE cp_files, ONLY: open_file, & + file_exists, & + close_file + USE physcon, ONLY: femtoseconds + USE mathconstants, ONLY: twopi + +#include "../base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse" + + #:include "rt_bse_macros.fypp" + + PUBLIC :: run_propagation_bse + + INTERFACE get_sigma + MODULE PROCEDURE get_sigma_dbcsr, get_sigma_real, get_sigma_complex + END INTERFACE + +CONTAINS + +! ************************************************************************************************** +!> \brief Runs the electron-only real time BSE propagation +!> \param qs_env Quickstep environment data, containing the config from input files +!> \param force_env Force environment data, entry point of the calculation +! ************************************************************************************************** + SUBROUTINE run_propagation_bse(qs_env, force_env) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(force_env_type), POINTER :: force_env + CHARACTER(len=*), PARAMETER :: routineN = 'run_propagation_bse' + TYPE(rtbse_env_type), POINTER :: rtbse_env + INTEGER :: i, j, k, handle + LOGICAL :: converged + REAL(kind=dp) :: metric, enum_re, enum_im, & + idempotence_dev, a_metric_1, a_metric_2 + + CALL timeset(routineN, handle) + + ! Run the initial SCF calculation / read SCF restart information + CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., consistent_energies=.FALSE.) + + ! Allocate all persistant storage and read input that does not need further processing + CALL create_rtbse_env(rtbse_env, qs_env, force_env) + + IF (rtbse_env%unit_nr > 0) CALL print_rtbse_header_info(rtbse_env) + + ! Initialize non-trivial values + ! - calculates the moment operators + ! - populates the initial density matrix + ! - reads the restart density if requested + ! - reads the field and moment trace from previous runs + ! - populates overlap and inverse overlap matrices + ! - calculates/populates the G0W0/KS Hamiltonian, respectively + ! - calculates the Hartree reference potential + ! - calculates the COHSEX reference self-energy + ! - prints some info about loaded files into the output + CALL initialize_rtbse_env(rtbse_env) + + ! Setup the time based on the starting step + ! Assumes identical dt between two runs + rtbse_env%sim_time = REAL(rtbse_env%sim_start, dp)*rtbse_env%sim_dt + ! Output 0 time moments and field + IF (.NOT. rtbse_env%restart_extracted) THEN + CALL output_field(rtbse_env) + CALL output_moments(rtbse_env, rtbse_env%rho) + END IF + + ! Do not apply the delta kick if we are doing a restart calculation + IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. (.NOT. rtbse_env%restart_extracted)) THEN + CALL apply_delta_pulse(rtbse_env) + END IF + + ! ********************** Start the time loop ********************** + DO i = rtbse_env%sim_start, rtbse_env%sim_nsteps + + ! Update the simulation time + rtbse_env%sim_time = REAL(i, dp)*rtbse_env%sim_dt + rtbse_env%sim_step = i + ! Start the enforced time reversal method + ! This method determines the density matrix at time (t+dt) by guessing the effective Hamiltonian at (t + dt) + ! and using the Hamiltonian at time (t), it propagates density from time (t) while ensuring that the density + ! at (t + dt/2) is the same for both forward and backwards propagation. Then, density at (t + dt) is + ! used to calculate the new Hamiltonian at (t+dt), which is then used to get the new propagator, and so on + ! until the density matrix does not change within certain limit + ! Pseudocode of the algorithm + ! rho_M = exp(-i S^(-1) H[rho(t)] dt/2) rho(t) exp(i H[rho(t)] S^(-1) dt/2) + ! rho(t+dt, 0) = rho_M + ! for j in 1,max_self_iter + ! 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) + ! if ||rho(t+dt,j) - rho(t+dt,j-1)|| < epsilon + ! break + ! *************** Determine rho_M *************** + ! Update the effective Hamiltonian + ! - sets the correct external field, updates Hartree and COHSEX terms to current density + CALL update_effective_ham(rtbse_env, rtbse_env%rho) + ! Convert the effective hamiltonian into the exponential + ! TODO : Store the result more explicitly - explicit dependence on input/output matrices + CALL ham_to_exp(rtbse_env) + ! Propagate the density to mid-point + CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rtbse_env%rho, rtbse_env%rho_M) + ! In initial iteration, copy rho_M to rho_new_last - that is our initial guess + DO j = 1, rtbse_env%n_spin + CALL cp_cfm_to_cfm(rtbse_env%rho_M(j), rtbse_env%rho_new_last(j)) + END DO + ! *********** Start the self-consistent loop ************************ + ! The guess Hamiltonian uses the guess density and field from the next timestep + rtbse_env%sim_time = REAL(i + 1, dp)*rtbse_env%sim_dt + rtbse_env%sim_step = i + 1 + converged = .FALSE. + ! Prints the grid for etrs info dumping - only for log level > medium + CALL print_etrs_info_header(rtbse_env) + DO k = 1, rtbse_env%etrs_max_iter + CALL update_effective_ham(rtbse_env, rtbse_env%rho_new_last) + CALL ham_to_exp(rtbse_env) + CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rtbse_env%rho_M, rtbse_env%rho_new) + ! *** Self-consistency check *** + metric = rho_metric(rtbse_env%rho_new, rtbse_env%rho_new_last, rtbse_env%n_spin) + ! ETRS info - only for log level > medium + CALL print_etrs_info(rtbse_env, k, metric) + IF (metric < rtbse_env%etrs_threshold) THEN + converged = .TRUE. + EXIT + ELSE + ! Copy rho_new to rho_new_last + DO j = 1, rtbse_env%n_spin + ! Leaving for free convergence + CALL cp_cfm_to_cfm(rtbse_env%rho_new(j), rtbse_env%rho_new_last(j)) + END DO + END IF + END DO + CALL get_electron_number(rtbse_env, rtbse_env%rho_new, enum_re, enum_im) + IF (.FALSE.) THEN + ! Not all of these are used, but they are all good metrics to check the convergence in problematic cases + ! TODO : Allow for conditional warning + CALL get_idempotence_deviation(rtbse_env, rtbse_env%rho_new, idempotence_dev) + DO j = 1, rtbse_env%n_spin + CALL cp_cfm_to_fm(rtbse_env%sigma_SEX(j), rtbse_env%real_workspace(1), rtbse_env%real_workspace(2)) + CALL antiherm_metric(real_fm=rtbse_env%real_workspace(1), imag_fm=rtbse_env%real_workspace(2), & + workspace=rtbse_env%rho_workspace, metric=a_metric_1) + CALL antiherm_metric(real_fm=rtbse_env%hartree_curr(j), & + workspace=rtbse_env%rho_workspace, metric=a_metric_2) + END DO + END IF + CALL print_timestep_info(rtbse_env, i, metric, enum_re, k) + CPASSERT(converged) + DO j = 1, rtbse_env%n_spin + CALL cp_cfm_to_cfm(rtbse_env%rho_new(j), rtbse_env%rho(j)) + END DO + ! Print the updated field + CALL output_field(rtbse_env) + ! If needed, print out the density matrix in MO basis + CALL output_mos_contravariant(rtbse_env, rtbse_env%rho, rtbse_env%rho_section) + ! Also handles outputting to memory + CALL output_moments(rtbse_env, rtbse_env%rho) + ! Output restart files, so that the restart starts at the following time index + CALL output_restart(rtbse_env, rtbse_env%rho, i + 1) + END DO + ! ********************** End the time loop ********************** + + ! Carry out the FT + CALL output_moments_ft(rtbse_env) + CALL output_polarizability(rtbse_env) + + ! Deallocate everything + CALL release_rtbse_env(rtbse_env) + + CALL timestop(handle) + END SUBROUTINE run_propagation_bse + +! ************************************************************************************************** +!> \brief Calculates the initial values, based on restart/scf density, and other non-trivial values +!> \param rtbse_env RT-BSE environment +!> \param qs_env Quickstep environment (needed for reference to previous calculations) +!> \author Stepan Marek (09.24) +! ************************************************************************************************** + SUBROUTINE initialize_rtbse_env(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + CHARACTER(len=*), PARAMETER :: routineN = "initialize_rtbse_env" + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments_dbcsr_p + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s + REAL(kind=dp), DIMENSION(:), POINTER :: custom_ref_point, & + occupations + REAL(kind=dp), DIMENSION(3) :: rpoint + INTEGER :: i, k, handle + + CALL timeset(routineN, handle) + + ! Get pointers to parameters from qs_env + CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, matrix_s=matrix_s) + + ! ****** START MOMENTS OPERATOR CALCULATION + ! Construct moments from dbcsr + NULLIFY (moments_dbcsr_p) + ALLOCATE (moments_dbcsr_p(3)) + DO k = 1, 3 + ! Make sure the pointer is empty + NULLIFY (moments_dbcsr_p(k)%matrix) + ! Allocate a new matrix that the pointer points to + ALLOCATE (moments_dbcsr_p(k)%matrix) + ! Create the matrix storage - matrix copies the structure of overlap matrix + CALL dbcsr_copy(moments_dbcsr_p(k)%matrix, matrix_s(1)%matrix) + END DO + ! Run the moment calculation + ! check for presence to prevent memory errors + NULLIFY (custom_ref_point) + ALLOCATE (custom_ref_point(3), source=0.0_dp) + rpoint(:) = 0.0_dp + CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, reference=use_mom_ref_coac, ref_point=custom_ref_point) + CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint) + ! Copy to full matrix + DO k = 1, 3 + ! Again, matrices are created from overlap template + CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, rtbse_env%moments(k)) + END DO + ! Now, repeat without reference point to get the moments for field + CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, reference=use_mom_ref_zero, ref_point=custom_ref_point) + DEALLOCATE (custom_ref_point) + CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint) + DO k = 1, 3 + CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, rtbse_env%moments_field(k)) + END DO + + ! Now can deallocate dbcsr matrices + DO k = 1, 3 + CALL dbcsr_release(moments_dbcsr_p(k)%matrix) + DEALLOCATE (moments_dbcsr_p(k)%matrix) + END DO + DEALLOCATE (moments_dbcsr_p) + ! ****** END MOMENTS OPERATOR CALCULATION + + ! ****** START INITIAL DENSITY MATRIX CALCULATION + ! Get the rho from fm_MOS + ! Uses real orbitals only - no kpoints + ALLOCATE (occupations(rtbse_env%n_ao)) + ! Iterate over both spins + DO i = 1, rtbse_env%n_spin + occupations(:) = 0.0_dp + occupations(1:rtbse_env%n_occ(i)) = 1.0_dp + ! Create real part + CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(i), rtbse_env%real_workspace(1)) + CALL cp_fm_column_scale(rtbse_env%real_workspace(1), occupations) + CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%real_workspace(1), bs_env%fm_mo_coeff_Gamma(i), & + 0.0_dp, rtbse_env%real_workspace(2)) + ! Sets imaginary part to zero + CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%rho(i)) + ! Save the reference value for the case of delta kick + CALL cp_cfm_to_cfm(rtbse_env%rho(i), rtbse_env%rho_orig(i)) + END DO + DEALLOCATE (occupations) + ! If the restart field is provided, overwrite rho from restart + IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN + CALL read_restart(rtbse_env) + END IF + ! ****** END INITIAL DENSITY MATRIX CALCULATION + ! ****** START MOMENTS TRACE LOAD + ! The moments are only loaded if the consistent save files exist + CALL read_moments(rtbse_env) + ! ****** END MOMENTS TRACE LOAD + + ! ****** START FIELD TRACE LOAD + ! The moments are only loaded if the consistent save files exist + CALL read_field(rtbse_env) + ! ****** END FIELD TRACE LOAD + + ! ****** START OVERLAP + INVERSE OVERLAP CALCULATION + CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, rtbse_env%S_fm) + CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%S_cfm) + CALL cp_fm_invert(rtbse_env%S_fm, rtbse_env%S_inv_fm) + ! ****** END OVERLAP + INVERSE OVERLAP CALCULATION + + ! ****** START SINGLE PARTICLE HAMILTONIAN CALCULATION + DO i = 1, rtbse_env%n_spin + IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN + ! G0W0 Hamiltonian + CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(i), rtbse_env%real_workspace(1)) + ! NOTE : Gamma point is not always the zero k-point + ! C * Lambda + CALL cp_fm_column_scale(rtbse_env%real_workspace(1), bs_env%eigenval_G0W0(:, 1, i)) + ! C * Lambda * C^T + CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%real_workspace(1), bs_env%fm_mo_coeff_Gamma(i), & + 0.0_dp, rtbse_env%real_workspace(2)) + ! S * C * Lambda * C^T + CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%S_fm, rtbse_env%real_workspace(2), & + 0.0_dp, rtbse_env%real_workspace(1)) + ! S * C * Lambda * C^T * S = H + CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%real_workspace(1), rtbse_env%S_fm, & + 0.0_dp, rtbse_env%real_workspace(2)) + CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_reference(i)) + ELSE + ! KS Hamiltonian + CALL cp_fm_to_cfm(msourcer=bs_env%fm_ks_Gamma(i), mtarget=rtbse_env%ham_reference(i)) + END IF + END DO + ! ****** END SINGLE PARTICLE HAMILTONIAN CALCULATION + + ! ****** START HARTREE POTENTIAL REFERENCE CALCULATION + ! Calculate Coulomb RI elements, necessary for Hartree calculation + CALL init_hartree(rtbse_env, rtbse_env%v_dbcsr) + IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN + ! In a non-HF calculation, copy the actual correlation part of the interaction + CALL copy_fm_to_dbcsr(bs_env%fm_W_MIC_freq_zero, rtbse_env%w_dbcsr) + ELSE + ! In HF, correlation is set to zero + CALL dbcsr_set(rtbse_env%w_dbcsr, 0.0_dp) + END IF + ! Add the Hartree to the screened_dbt tensor - now W = V + W^c + CALL dbcsr_add(rtbse_env%w_dbcsr, rtbse_env%v_dbcsr, 1.0_dp, 1.0_dp) + CALL dbt_copy_matrix_to_tensor(rtbse_env%w_dbcsr, rtbse_env%screened_dbt) + ! Calculate the original Hartree potential + ! Uses rho_orig - same as rho for initial run but different for continued run + DO i = 1, rtbse_env%n_spin + CALL get_hartree_local(rtbse_env, rtbse_env%rho_orig(i), rtbse_env%hartree_curr(i)) + ! Scaling by spin degeneracy + CALL cp_fm_scale(rtbse_env%spin_degeneracy, rtbse_env%hartree_curr(i)) + ! Subtract the reference from the reference Hamiltonian + CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(i), mtarget=rtbse_env%ham_workspace(1)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), & + CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1)) + END DO + ! ****** END HARTREE POTENTIAL REFERENCE CALCULATION + + ! ****** START COHSEX REFERENCE CALCULATION + ! Calculate the COHSEX starting energies + IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN + ! Subtract the v_xc from COH part of the self-energy, as V_xc is also not updated during the timestepping + DO i = 1, rtbse_env%n_spin + ! TODO : Allow no COH calculation for static screening + CALL get_sigma(rtbse_env, rtbse_env%sigma_COH(i), -0.5_dp, rtbse_env%S_inv_fm) + ! Copy and subtract from the complex reference hamiltonian + CALL cp_fm_to_cfm(msourcer=rtbse_env%sigma_COH(i), mtarget=rtbse_env%ham_workspace(1)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), & + CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1)) + ! Calculate exchange part - TODO : should this be applied for different spins? - TEST with O2 HF propagation? + ! So far only closed shell tested + ! Uses rho_orig - same as rho for initial run but different for continued run + CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(i), -1.0_dp, rtbse_env%rho_orig(i)) + ! Subtract from the complex reference Hamiltonian + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), & + CMPLX(-1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(i)) + END DO + ELSE + ! KS Hamiltonian - use time-dependent Fock exchange + DO i = 1, rtbse_env%n_spin + CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(i), -1.0_dp, rtbse_env%rho_orig(i)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), & + CMPLX(-1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(i)) + END DO + END IF + ! ****** END COHSEX REFERENCE CALCULATION + + CALL timestop(handle) + END SUBROUTINE initialize_rtbse_env + +! ************************************************************************************************** +!> \brief Custom reimplementation of the delta pulse routines +!> \param rtbse_env RT-BSE environment +!> \author Stepan Marek (09.24) +! ************************************************************************************************** + SUBROUTINE apply_delta_pulse(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + CHARACTER(len=*), PARAMETER :: routineN = "apply_delta_pulse" + REAL(kind=dp) :: intensity, metric + REAL(kind=dp), DIMENSION(3) :: kvec + INTEGER :: i, k, handle + + CALL timeset(routineN, handle) + + ! Report application + IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A27)') ' RTBSE| Applying delta pulse' + ! Extra minus for the propagation of density + intensity = -rtbse_env%dft_control%rtp_control%delta_pulse_scale + metric = 0.0_dp + kvec(:) = rtbse_env%dft_control%rtp_control%delta_pulse_direction(:) + IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A38,E14.4E3,E14.4E3,E14.4E3)') & + " RTBSE| Delta pulse elements (a.u.) : ", intensity*kvec(:) + ! So far no spin dependence, but can be added by different structure of delta pulse + CALL cp_fm_set_all(rtbse_env%real_workspace(1), 0.0_dp) + DO k = 1, 3 + CALL cp_fm_scale_and_add(1.0_dp, rtbse_env%real_workspace(1), & + kvec(k), rtbse_env%moments_field(k)) + END DO + ! enforce hermiticity of the effective Hamiltonian + CALL cp_fm_transpose(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2)) + CALL cp_fm_scale_and_add(0.5_dp, rtbse_env%real_workspace(1), & + 0.5_dp, rtbse_env%real_workspace(2)) + ! Prepare the exponential/exponent for propagation + IF (rtbse_env%mat_exp_method == do_bch) THEN + ! Multiply by the S_inv matrix - in the classic ordering + CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + intensity, rtbse_env%S_inv_fm, rtbse_env%real_workspace(1), & + 0.0_dp, rtbse_env%real_workspace(2)) + DO i = 1, rtbse_env%n_spin + ! Sets real part to zero + CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_effective(i)) + END DO + ELSE IF (rtbse_env%mat_exp_method == do_exact) THEN + DO i = 1, rtbse_env%n_spin + CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(1), mtarget=rtbse_env%ham_workspace(i)) + CALL cp_cfm_gexp(rtbse_env%ham_workspace(i), rtbse_env%S_cfm, rtbse_env%ham_effective(i), & + CMPLX(0.0, intensity, kind=dp), rtbse_env%rho_workspace) + END DO + END IF + ! Propagate the density by the effect of the delta pulse + CALL propagate_density(rtbse_env, rtbse_env%ham_effective, rtbse_env%rho, rtbse_env%rho_new) + metric = rho_metric(rtbse_env%rho_new, rtbse_env%rho, rtbse_env%n_spin) + IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, ('(A42,E38.8E3)')) " RTBSE| Metric difference after delta kick", metric + ! Copy the new density to the old density + DO i = 1, rtbse_env%n_spin + CALL cp_cfm_to_cfm(rtbse_env%rho_new(i), rtbse_env%rho(i)) + END DO + + CALL timestop(handle) + END SUBROUTINE apply_delta_pulse +! ************************************************************************************************** +!> \brief Determines the metric for the density matrix, used for convergence criterion +!> \param rho_new Array of new density matrices (one for each spin index) +!> \param rho_old Array of old density matrices (one for each spin index) +!> \param nspin Number of spin indices +!> \param workspace_opt Optionally provide external workspace to save some allocation time +! ************************************************************************************************** + FUNCTION rho_metric(rho_new, rho_old, nspin, workspace_opt) RESULT(metric) + TYPE(cp_cfm_type), DIMENSION(:), POINTER, INTENT(IN):: rho_new, & + rho_old + INTEGER, INTENT(IN) :: nspin + TYPE(cp_cfm_type), POINTER, OPTIONAL :: workspace_opt + TYPE(cp_cfm_type) :: workspace + REAL(kind=dp) :: metric + REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: partial_metric + INTEGER :: j + COMPLEX(kind=dp) :: scale_factor + + ALLOCATE (partial_metric(nspin)) + + ! Only allocate/deallocate storage if required + IF (PRESENT(workspace_opt)) THEN + workspace = workspace_opt + ELSE + CALL cp_cfm_create(workspace, rho_new(1)%matrix_struct) + END IF + scale_factor = 1.0 + DO j = 1, nspin + CALL cp_cfm_to_cfm(rho_new(j), workspace) + ! Get the difference in the resulting matrix + CALL cp_cfm_scale_and_add(scale_factor, workspace, -scale_factor, rho_old(j)) + ! Now, get the relevant number + partial_metric(j) = cp_cfm_norm(workspace, 'M') + END DO + metric = 0.0_dp + ! For more than one spin, do Cartesian sum of the different spin norms + DO j = 1, nspin + metric = metric + partial_metric(j)*partial_metric(j) + END DO + metric = SQRT(metric) + ! Deallocate workspace + IF (.NOT. PRESENT(workspace_opt)) CALL cp_cfm_release(workspace) + DEALLOCATE (partial_metric) + END FUNCTION + +! ************************************************************************************************** +!> \brief Determines the metric of the antihermitian part of the matrix +!> \param real_fm Real part of the full matrix +!> \param imag_fm Imaginary part of the full matrix +! ************************************************************************************************** + SUBROUTINE antiherm_metric(real_fm, imag_fm, workspace, metric) + TYPE(cp_fm_type), INTENT(IN) :: real_fm + TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: imag_fm + REAL(kind=dp), INTENT(OUT) :: metric + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: workspace + COMPLEX(kind=dp) :: complex_one + + ! Get the complex and complex conjugate matrix + IF (PRESENT(imag_fm)) THEN + CALL cp_fm_to_cfm(real_fm, imag_fm, workspace(1)) + ELSE + CALL cp_fm_to_cfm(msourcer=real_fm, mtarget=workspace(1)) + END IF + CALL cp_cfm_transpose(workspace(1), "C", workspace(2)) + ! Subtract these, and get the metric + complex_one = CMPLX(1.0, 0.0, kind=dp) + CALL cp_cfm_scale_and_add(complex_one, workspace(1), -complex_one, workspace(2)) + metric = cp_cfm_norm(workspace(1), "M") + END SUBROUTINE + +! ************************************************************************************************** +!> \brief For Taylor and Exact exp_method, calculates the matrix exponential of the +!> effective Hamiltonian. For BCH, calculates just the effective Hamiltonian. For other methods, +!> aborts the execution, as they are not implemented yet. +!> \param rtbse_env Entry point of the calculation. Uses rho_workspace for Taylor and BCH. For exact, +!> uses complex_workspace, complex_ham, complex_s, real_eigvals and exp_eigvals. +!> Results are stored in ham_workspace. +! ************************************************************************************************** + SUBROUTINE ham_to_exp(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + CHARACTER(len=*), PARAMETER :: routineN = "ham_to_exp" + INTEGER :: j, handle + CALL timeset(routineN, handle) + DO j = 1, rtbse_env%n_spin + IF (rtbse_env%mat_exp_method == do_bch) THEN + ! In Taylor and BCH, we first evaluate the entire exponent and then evaluate exponential in series + ! In order to produce correct result, need to remultiply by inverse overlap matrix + CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%S_inv_fm, rtbse_env%ham_effective(j), & + 0.0_dp, rtbse_env%rho_workspace(1)) + + ! The evolution of density matrix is derived from the right multiplying term + ! Imaginary part of the exponent = -real part of the matrix + CALL cp_cfm_scale(CMPLX(0.0, -rtbse_env%sim_dt/2, kind=dp), rtbse_env%rho_workspace(1)) + ! In BCH, exponential is not calculated explicitly, but the propagation is solved in series + CALL cp_cfm_to_cfm(rtbse_env%rho_workspace(1), rtbse_env%ham_workspace(j)) + ELSE IF (rtbse_env%mat_exp_method == do_exact) THEN + CALL cp_cfm_gexp(rtbse_env%ham_effective(j), rtbse_env%S_cfm, rtbse_env%ham_workspace(j), & + CMPLX(0.0, -rtbse_env%sim_dt/2, kind=dp), rtbse_env%rho_workspace) + ELSE + CPABORT("Only BCH and Taylor matrix exponentiation implemented") + END IF + END DO + + CALL timestop(handle) + END SUBROUTINE +! ************************************************************************************************** +!> \brief Updates the effective Hamiltonian, given a density matrix rho +!> \param rtbse_env Entry point of the calculation - contains current state of variables +!> \param qs_env QS env +!> \param rho Real and imaginary parts ( + spin) of the density at current time +! ************************************************************************************************** + SUBROUTINE update_effective_ham(rtbse_env, rho) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho + CHARACTER(len=*), PARAMETER :: routineN = "update_effective_ham" + INTEGER :: k, j, nspin, handle + + CALL timeset(routineN, handle) + ! Shorthand + nspin = rtbse_env%n_spin + ! Reset the effective Hamiltonian to KS Hamiltonian + G0W0 - reference COHSEX - reference Hartree + DO j = 1, nspin + ! Sets the imaginary part to zero + CALL cp_cfm_to_cfm(rtbse_env%ham_reference(j), rtbse_env%ham_effective(j)) + END DO + ! Determine the field at current time + IF (rtbse_env%dft_control%apply_efield_field) THEN + CALL make_field(rtbse_env%dft_control, rtbse_env%field, rtbse_env%sim_step, rtbse_env%sim_time) + ELSE + ! No field + rtbse_env%field(:) = 0.0_dp + END IF + DO j = 1, nspin + DO k = 1, 3 + ! Minus sign due to charge of electrons + CALL cp_fm_to_cfm(msourcer=rtbse_env%moments_field(k), mtarget=rtbse_env%ham_workspace(1)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), & + CMPLX(rtbse_env%field(k), 0.0, kind=dp), rtbse_env%ham_workspace(1)) + END DO + IF (rtbse_env%ham_reference_type == rtp_bse_ham_g0w0) THEN + ! Add the COH part - so far static but can be dynamic in principle throught the W updates + CALL get_sigma(rtbse_env, rtbse_env%sigma_COH(j), -0.5_dp, rtbse_env%S_inv_fm) + CALL cp_fm_to_cfm(msourcer=rtbse_env%sigma_COH(j), mtarget=rtbse_env%ham_workspace(1)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), & + CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1)) + END IF + ! Calculate the (S)EX part - based on provided rho + ! iGW = - rho W + CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX(j), -1.0_dp, rho(j)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), & + CMPLX(1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(j)) + ! Calculate Hartree potential + ! Hartree potential is scaled by number of electrons in each MO - spin degeneracy + CALL get_hartree_local(rtbse_env, rho(j), & + rtbse_env%hartree_curr(j)) + CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(j), mtarget=rtbse_env%ham_workspace(1)) + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_effective(j), & + CMPLX(rtbse_env%spin_degeneracy, 0.0, kind=dp), rtbse_env%ham_workspace(1)) + ! Enforce hermiticity of the effective Hamiltonian + ! Important components without forced Hermiticity - moments matrix, sigma matrices, Hartree matrix + ! single particle Ham + CALL cp_cfm_transpose(rtbse_env%ham_effective(j), 'C', rtbse_env%ham_workspace(1)) + CALL cp_cfm_scale_and_add(CMPLX(0.5, 0.0, kind=dp), rtbse_env%ham_effective(j), & + CMPLX(0.5, 0.0, kind=dp), rtbse_env%ham_workspace(1)) + END DO + CALL timestop(handle) + END SUBROUTINE update_effective_ham + +! ************************************************************************************************** +!> \brief Does the BCH iterative determination of the exponential +!> \param propagator_matrix Matrix X which is to be exponentiated +!> \param target_matrix Matrix Y which the exponential acts upon +!> \param result_matrix Propagated matrix +!> \param workspace Matrices dedicated for work, 4 fm matrices with dimensions of X required +!> \param threshold_opt Optionally, a threshold under which the iteration is considered converged (default 1e-10) +!> \param max_iter_opt Optionally, maximum number of BCH iterations (default 20) +! ************************************************************************************************** + SUBROUTINE bch_propagate(propagator_matrix, target_matrix, result_matrix, workspace, threshold_opt, max_iter_opt) + ! Array of complex propagator matrix X, such that + ! the propagated matrix will follow Y' = e^X Y e^(-X), for each spin + ! effect of e^(-X) is calculated - provide the X on the left hand side + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: propagator_matrix + ! Matrix Y to be propagated into matrix Y' + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: target_matrix + ! Matrix Y' is stored here on exit + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: result_matrix, workspace + ! Threshold for the metric which decides when to truncate the BCH expansion + REAL(kind=dp), OPTIONAL :: threshold_opt + INTEGER, OPTIONAL :: max_iter_opt + CHARACTER(len=*), PARAMETER :: routineN = "bch_propagate" + REAL(kind=dp) :: threshold, prefactor, metric + INTEGER :: max_iter, i, n_spin, n_ao, k, & + w_stride, handle + LOGICAL :: converged + CHARACTER(len=77) :: error + + CALL timeset(routineN, handle) + + converged = .FALSE. + + IF (PRESENT(threshold_opt)) THEN + threshold = threshold_opt + ELSE + threshold = 1.0e-10 + END IF + + IF (PRESENT(max_iter_opt)) THEN + max_iter = max_iter_opt + ELSE + max_iter = 20 + END IF + + n_spin = SIZE(target_matrix) + n_ao = 0 + CALL cp_cfm_get_info(target_matrix(1), nrow_global=n_ao) + w_stride = n_spin + + ! Initiate + DO i = 1, n_spin + CALL cp_cfm_to_cfm(target_matrix(i), result_matrix(i)) + CALL cp_cfm_to_cfm(target_matrix(i), workspace(i)) + END DO + + ! Start the BCH iterations + ! So far, no spin mixing terms + DO k = 1, max_iter + prefactor = 1.0_dp/REAL(k, kind=dp) + DO i = 1, n_spin + CALL parallel_gemm("N", "N", n_ao, n_ao, n_ao, & + CMPLX(prefactor, 0.0, kind=dp), propagator_matrix(i), workspace(i), & + CMPLX(0.0, 0.0, kind=dp), workspace(i + w_stride)) + CALL parallel_gemm("N", "C", n_ao, n_ao, n_ao, & + CMPLX(prefactor, 0.0, kind=dp), workspace(i), propagator_matrix(i), & + CMPLX(1.0, 0.0, kind=dp), workspace(i + w_stride)) + ! Add to the result + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), result_matrix(i), & + CMPLX(1.0, 0.0, kind=dp), workspace(i + w_stride)) + END DO + metric = rho_metric(workspace(w_stride + 1:), workspace(1:w_stride), n_spin) + IF (metric <= threshold) THEN + converged = .TRUE. + EXIT + ELSE + DO i = 1, n_spin + CALL cp_cfm_to_cfm(workspace(i + w_stride), workspace(i)) + END DO + END IF + END DO + IF (.NOT. converged) THEN + WRITE (error, '(A35,E13.4E3,A16,E13.4E3)') "BCH did not converge, BCH Metric : ", & + metric, "BCH Threshold : ", threshold + CPABORT(error) + END IF + + CALL timestop(handle) + END SUBROUTINE bch_propagate + +! ************************************************************************************************** +!> \brief Updates the density in rtbse_env, using the provided exponential +!> The new density is saved to a different matrix, which enables for comparison of matrices +!> \param rtbse_env Entry point of the calculation - contains current state of variables +!> \param exponential Real and imaginary parts ( + spin) of the exponential propagator +! ************************************************************************************************** + SUBROUTINE propagate_density(rtbse_env, exponential, rho_old, rho_new) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: exponential, & + rho_old, & + rho_new + CHARACTER(len=*), PARAMETER :: routineN = "propagate_density" + INTEGER :: j, handle + + CALL timeset(routineN, handle) + IF (rtbse_env%mat_exp_method == do_exact) THEN + ! For these methods, exponential is explicitly constructed + DO j = 1, rtbse_env%n_spin + ! rho * (exp^dagger) + CALL parallel_gemm("N", "C", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + CMPLX(1.0, 0.0, kind=dp), rho_old(j), exponential(j), & + CMPLX(0.0, 0.0, kind=dp), rtbse_env%rho_workspace(1)) + ! exp * rho * (exp^dagger) + CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + CMPLX(1.0, 0.0, kind=dp), exponential(j), rtbse_env%rho_workspace(1), & + CMPLX(0.0, 0.0, kind=dp), rho_new(j)) + END DO + ELSE IF (rtbse_env%mat_exp_method == do_bch) THEN + ! Same number of iterations as ETRS + CALL bch_propagate(exponential, rho_old, rho_new, rtbse_env%rho_workspace, threshold_opt=rtbse_env%exp_accuracy, & + max_iter_opt=rtbse_env%etrs_max_iter) + ELSE + CPABORT("Only BCH and exact matrix exponentiation implemented.") + END IF + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Outputs the number of electrons in the system from the density matrix +!> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3) +!> \param rtbse_env Entry point - rtbse environment +!> \param rho Density matrix in AO basis +!> \param electron_n_re Real number of electrons +!> \param electron_n_im Imaginary number of electrons, which can arise from numerical non-hermiticity +! ************************************************************************************************** + SUBROUTINE get_electron_number(rtbse_env, rho, electron_n_re, electron_n_im) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho + REAL(kind=dp), INTENT(OUT) :: electron_n_re, electron_n_im + COMPLEX(kind=dp) :: electron_n_buffer + INTEGER :: j + + electron_n_re = 0.0_dp + electron_n_im = 0.0_dp + CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%rho_workspace(1)) + DO j = 1, rtbse_env%n_spin + CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rho(j), electron_n_buffer) + electron_n_re = electron_n_re + REAL(electron_n_buffer, kind=dp) + electron_n_im = electron_n_im + REAL(AIMAG(electron_n_buffer), kind=dp) + END DO + ! Scale by spin degeneracy + electron_n_re = electron_n_re*rtbse_env%spin_degeneracy + electron_n_im = electron_n_im*rtbse_env%spin_degeneracy + END SUBROUTINE get_electron_number +! ************************************************************************************************** +!> \brief Outputs the deviation from idempotence of density matrix +!> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3) +!> \param rtbse_env Entry point - rtbse environment +!> \param rho Density matrix in AO basis +!> \param electron_n_re Real number of electrons +!> \param electron_n_im Imaginary number of electrons, which can arise from numerical non-hermiticity +! ************************************************************************************************** + SUBROUTINE get_idempotence_deviation(rtbse_env, rho, deviation_metric) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho + REAL(kind=dp), INTENT(OUT) :: deviation_metric + COMPLEX(kind=dp) :: buffer_1, buffer_2 + REAL(kind=dp) :: buffer_dev + INTEGER :: j + + deviation_metric = 0.0_dp + buffer_dev = 0.0_dp + ! First, determine Tr(S * rho_re) + i Tr (S * rho_im) + CALL cp_fm_to_cfm(msourcer=rtbse_env%S_fm, mtarget=rtbse_env%rho_workspace(1)) + DO j = 1, rtbse_env%n_spin + CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rho(j), buffer_1) + buffer_dev = buffer_dev + REAL(ABS(buffer_1)*ABS(buffer_1)) + END DO + ! Now, determine Tr(S * rho_re * S * rho_re) - Tr(S * rho_im * S * rho_im) + 2i Tr(S * rho_re * S * rho_im) + DO j = 1, rtbse_env%n_spin + ! S * rho + CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%S_fm, rho(j), & + 0.0_dp, rtbse_env%rho_workspace(2)) + ! rho * S * rho + CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + CMPLX(1.0, 0.0, kind=dp), rho(j), rtbse_env%rho_workspace(2), & + CMPLX(0.0, 0.0, kind=dp), rtbse_env%rho_workspace(3)) + ! Tr (S * rho * S * rho) + CALL cp_cfm_trace(rtbse_env%rho_workspace(1), rtbse_env%rho_workspace(3), buffer_2) + deviation_metric = deviation_metric + REAL(ABS(buffer_2)*ABS(buffer_2)) + END DO + deviation_metric = SQRT(deviation_metric) - SQRT(buffer_dev) + END SUBROUTINE get_idempotence_deviation + +! ************************************************************************************************** +!> \brief Calculates the self-energy by contraction of screened potential, for complex density +!> \note Can be used for both the Coulomb hole part and screened exchange part +!> \param rtbse_env Quickstep environment data, entry point of the calculation +!> \param sigma_cfm Pointer to the self-energy full matrix, which is overwritten by this routine +!> \param greens_cfm Pointer to the Green's function matrix, which is used as input data +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + SUBROUTINE get_sigma_complex(rtbse_env, sigma_cfm, prefactor_opt, greens_cfm) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type) :: sigma_cfm ! resulting self energy + REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt + TYPE(cp_cfm_type), INTENT(IN) :: greens_cfm ! matrix to contract with RI_W + REAL(kind=dp) :: prefactor + + prefactor = 1.0_dp + IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt + + ! Carry out the sigma part twice + ! Real part + CALL cp_cfm_to_fm(msource=greens_cfm, mtargetr=rtbse_env%real_workspace(1)) + CALL get_sigma(rtbse_env, rtbse_env%real_workspace(2), prefactor, rtbse_env%real_workspace(1)) + CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), mtarget=rtbse_env%ham_workspace(1)) + ! Imaginary part + CALL cp_cfm_to_fm(msource=greens_cfm, mtargeti=rtbse_env%real_workspace(1)) + CALL get_sigma(rtbse_env, rtbse_env%real_workspace(2), prefactor, rtbse_env%real_workspace(1)) + CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace(2), mtarget=sigma_cfm) + ! Add the real part + CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), sigma_cfm, & + CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1)) + + END SUBROUTINE get_sigma_complex +! ************************************************************************************************** +!> \brief Calculates the self-energy by contraction of screened potential, for complex density +!> \note Can be used for both the Coulomb hole part and screened exchange part +!> \param rtbse_env Quickstep environment data, entry point of the calculation +!> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine +!> \param greens_fm Pointer to the Green's function matrix, which is used as input data +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + SUBROUTINE get_sigma_real(rtbse_env, sigma_fm, prefactor_opt, greens_fm) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_fm_type) :: sigma_fm ! resulting self energy + REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt + TYPE(cp_fm_type), INTENT(IN) :: greens_fm ! matrix to contract with RI_W + REAL(kind=dp) :: prefactor + + prefactor = 1.0_dp + IF (PRESENT(prefactor_opt)) prefactor = prefactor_opt + + ! Carry out the sigma part twice + ! Convert to dbcsr + CALL copy_fm_to_dbcsr(greens_fm, rtbse_env%rho_dbcsr) + CALL get_sigma(rtbse_env, sigma_fm, prefactor, rtbse_env%rho_dbcsr) + + END SUBROUTINE get_sigma_real +! ************************************************************************************************** +!> \brief Calculates the self-energy by contraction of screened potential +!> \note Can be used for both the Coulomb hole part and screened exchange part +!> \param qs_env Quickstep environment data, entry point of the calculation +!> \param greens_fm Pointer to the Green's function matrix, which is used as input data +!> \param sigma_fm Pointer to the self-energy full matrix, which is overwritten by this routine +!> \author Stepan Marek +!> \date 01.2024 +! ************************************************************************************************** + SUBROUTINE get_sigma_dbcsr(rtbse_env, sigma_fm, prefactor_opt, greens_dbcsr) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_fm_type) :: sigma_fm ! resulting self energy + REAL(kind=dp), INTENT(IN), OPTIONAL :: prefactor_opt + TYPE(dbcsr_type) :: greens_dbcsr + CHARACTER(len=*), PARAMETER :: routineN = 'get_sigma' + REAL(kind=dp) :: prefactor + TYPE(dbcsr_type) :: sigma_dbcsr + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (PRESENT(prefactor_opt)) THEN + prefactor = prefactor_opt + ELSE + prefactor = 1.0_dp + END IF + + ! Three-centre integrals are obtained from build_3c_integrals, from qs_tensors + ! These should use sparcity, while W and Sigma can be full matrices + ! The summation is carried out by dbt library - dbt_contract in dbt_api + ! The building of the tensors might be a bit hard, because it requires a lot of parallel information + ! Probably just use the tensors already present in bs_env? They seem to be mostly work tensors + ! Create by template + CALL dbt_contract(alpha=1.0_dp, & + tensor_1=rtbse_env%screened_dbt, & + tensor_2=rtbse_env%t_3c_w, & + beta=0.0_dp, & + tensor_3=rtbse_env%t_3c_work_RI_AO__AO, & + contract_1=[2], notcontract_1=[1], map_1=[1], & + contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3])!,& + !filter_eps=bs_env%eps_filter) + ! t_work1 now contains B^P_(nu beta) = sum _ Q W _ (PQ) (iomega = 0) (Q| nu beta) + ! Next step is to convert the greens full matrix to dbcsr matrix + CALL dbt_copy_matrix_to_tensor(greens_dbcsr, rtbse_env%greens_dbt) + ! Then contract it + ! no scaling applied - this has to be applied externally + CALL dbt_contract(alpha=1.0_dp, & + tensor_1=rtbse_env%t_3c_work_RI_AO__AO, & + tensor_2=rtbse_env%greens_dbt, & + beta=0.0_dp, & + tensor_3=rtbse_env%t_3c_work2_RI_AO__AO, & + contract_1=[2], notcontract_1=[1, 3], map_1=[1, 3], & + contract_2=[2], notcontract_2=[1], map_2=[2]) + ! workspace 2 now contains C ^ P _ (mu beta) sum _ nu B ^ P _ (nu beta) g _ (mu nu) + CALL dbt_contract(alpha=prefactor, & + tensor_1=rtbse_env%t_3c_w, & + tensor_2=rtbse_env%t_3c_work2_RI_AO__AO, & + beta=0.0_dp, & + tensor_3=rtbse_env%sigma_dbt, & + contract_1=[1, 3], notcontract_1=[2], map_1=[1], & + contract_2=[1, 2], notcontract_2=[3], map_2=[2])!,& + !filter_eps=bs_env%eps_filter) + ! Finally, convert the COH tensor to matrix and then to fm matrix + CALL dbcsr_create(sigma_dbcsr, name="sigma", template=rtbse_env%rho_dbcsr) + CALL dbt_copy_tensor_to_matrix(rtbse_env%sigma_dbt, sigma_dbcsr) + CALL copy_dbcsr_to_fm(sigma_dbcsr, sigma_fm) + CALL dbcsr_release(sigma_dbcsr) + ! Clear workspaces - saves memory? + CALL dbt_clear(rtbse_env%t_3c_work_RI_AO__AO) + CALL dbt_clear(rtbse_env%t_3c_work2_RI_AO__AO) + CALL dbt_clear(rtbse_env%sigma_dbt) + CALL dbt_clear(rtbse_env%greens_dbt) + CALL timestop(handle) + + END SUBROUTINE get_sigma_dbcsr +! ************************************************************************************************** +!> \brief Creates the RI matrix and populates it with correct values +!> \note Tensor contains Hartree elements in the auxiliary basis +!> \param qs_env Quickstep environment - entry point of calculation +!> \author Stepan Marek +!> \date 01.2024 +! ************************************************************************************************** + SUBROUTINE init_hartree(rtbse_env, v_dbcsr) + TYPE(rtbse_env_type), POINTER, INTENT(IN) :: rtbse_env + TYPE(dbcsr_type) :: v_dbcsr + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(libint_potential_type) :: coulomb_op + TYPE(cp_fm_type) :: V_fm + TYPE(cp_fm_type) :: metric_fm + TYPE(cp_fm_type) :: metric_inv_fm, & + work_fm + TYPE(dbcsr_type), DIMENSION(:), ALLOCATABLE :: V_dbcsr_a, & + metric_dbcsr + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: nl_2c + + bs_env => rtbse_env%bs_env + + ! Allocate for bare Hartree term + ALLOCATE (V_dbcsr_a(1)) + ALLOCATE (metric_dbcsr(1)) + CALL dbcsr_create(V_dbcsr_a(1), name="Hartree_dbcsr", template=bs_env%mat_RI_RI%matrix) + CALL dbcsr_create(metric_dbcsr(1), name="RI_metric_dbcsr", template=bs_env%mat_RI_RI%matrix) + + ! Calculate full coulomb RI basis elements - V _ (PQ) matrix + NULLIFY (nl_2c) + CALL build_2c_neighbor_lists(nl_2c, bs_env%basis_set_RI, bs_env%basis_set_RI, & + coulomb_op, "Coulomb_neighbor_2c_list", rtbse_env%qs_env, & + sym_ij=.FALSE., molecular=.TRUE.) + CALL build_2c_integrals(V_dbcsr_a, bs_env%eps_filter, rtbse_env%qs_env, nl_2c, & + bs_env%basis_set_RI, bs_env%basis_set_RI, coulomb_op, & + do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI) + ! Calculate the RI metric elements + ! nl_2c is automatically rewritten (even reallocated) in this routine + CALL build_2c_neighbor_lists(nl_2c, bs_env%basis_set_RI, bs_env%basis_set_RI, & + bs_env%ri_metric, "Metric_neighbor_2c_list", rtbse_env%qs_env, & + sym_ij=.FALSE., molecular=.TRUE.) + CALL build_2c_integrals(metric_dbcsr, bs_env%eps_filter, rtbse_env%qs_env, nl_2c, & + bs_env%basis_set_RI, bs_env%basis_set_RI, bs_env%ri_metric, & + do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI) + ! nl_2c no longer needed + CALL release_neighbor_list_sets(nl_2c) + CALL cp_fm_create(metric_fm, bs_env%fm_RI_RI%matrix_struct) + CALL cp_fm_set_all(metric_fm, 0.0_dp) + CALL cp_fm_create(metric_inv_fm, bs_env%fm_RI_RI%matrix_struct) + CALL cp_fm_set_all(metric_inv_fm, 0.0_dp) + CALL cp_fm_create(work_fm, bs_env%fm_RI_RI%matrix_struct) + CALL cp_fm_set_all(work_fm, 0.0_dp) + CALL copy_dbcsr_to_fm(metric_dbcsr(1), metric_fm) + CALL cp_fm_invert(metric_fm, metric_inv_fm) + CALL cp_fm_create(V_fm, bs_env%fm_RI_RI%matrix_struct) + CALL cp_fm_set_all(V_fm, 0.0_dp) + ! Multiply by the inverse from each side (M^-1 is symmetric) + CALL cp_dbcsr_sm_fm_multiply(V_dbcsr_a(1), metric_inv_fm, & + work_fm, bs_env%n_RI) + CALL parallel_gemm("N", "N", bs_env%n_RI, bs_env%n_RI, bs_env%n_RI, & + 1.0_dp, metric_inv_fm, work_fm, 0.0_dp, V_fm) + ! Now, create the tensor from the matrix + ! First, convert full matrix to dbcsr + CALL dbcsr_clear(V_dbcsr_a(1)) + CALL copy_fm_to_dbcsr(V_fm, V_dbcsr_a(1)) + CALL dbcsr_create(v_dbcsr, "Hartree ri", V_dbcsr_a(1)) + CALL dbcsr_copy(v_dbcsr, V_dbcsr_a(1)) + ! Create and copy distinctly, so that unnecessary objects can be destroyed + ! Destroy all unnecessary matrices + CALL dbcsr_release(V_dbcsr_a(1)) + CALL dbcsr_release(metric_dbcsr(1)) + DEALLOCATE (V_dbcsr_a) + DEALLOCATE (metric_dbcsr) + CALL cp_fm_release(V_fm) + ! CALL cp_fm_release(metric_fm(1,1)) + CALL cp_fm_release(metric_fm) + ! DEALLOCATE(metric_fm) + CALL cp_fm_release(work_fm) + CALL cp_fm_release(metric_inv_fm) + END SUBROUTINE +! ************************************************************************************************** +!> \brief Calculates the Hartree matrix in the atomic orbital basis, given a density matrix, in local arrays +!> Calculates the values for single spin species present in given rho +!> \param qs_env Entry point +!> \param rtbse_env Entry point of GWBSE - uses rho_dbcsr and some complex_workspace +!> \param rho_ao Density matrix in ao basis +!> \param v_ao Overwritten by the Hartree matrix in the atomic orbital basis +!> \author Stepan Marek +!> \date 02.2024 +! ************************************************************************************************** + SUBROUTINE get_hartree_local(rtbse_env, rho_ao, v_ao) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), INTENT(IN) :: rho_ao + TYPE(cp_fm_type) :: v_ao + CHARACTER(len=*), PARAMETER :: routineN = "get_hartree_local" + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(dbcsr_iterator_type) :: iterator_matrix + INTEGER :: i, j, k, n, nblocks, ind_1, ind_2, row_offset, col_offset, & + row_size, col_size, j_n_AO, k_n_AO, i_n_RI, n_RI, & + ri_offset, ind_i, handle + REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: Pvector, Qvector + REAL(kind=dp), DIMENSION(:, :), POINTER :: block_matrix + REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c + TYPE(dbcsr_type) :: v_dbcsr_ao + + ! No memory optimisation so far - calculate all 3cs an all ranks + ! Importantly - dbcsr blocks are ordered by atoms - i.e. ethene with 6 atoms will have 6x6 block structure + ! Number of basis states on each basis set is known is post_scf_bandstructure env + + CALL timeset(routineN, handle) + ! Get the relevant environments + CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, para_env=para_env) + n_RI = rtbse_env%n_RI + int_3c => rtbse_env%int_3c_array + + ! Allocate the Q and Pvector on each rank + ALLOCATE (Qvector(rtbse_env%n_RI), source=0.0_dp) + ALLOCATE (Pvector(rtbse_env%n_RI), source=0.0_dp) + + ! First step - analyze the structure of copied dbcsr matrix on all ranks + CALL dbcsr_clear(rtbse_env%rho_dbcsr) + ! Only the real part of the density matrix contributes + CALL cp_cfm_to_fm(msource=rho_ao, mtargetr=rtbse_env%real_workspace(1)) + CALL copy_fm_to_dbcsr(rtbse_env%real_workspace(1), rtbse_env%rho_dbcsr) + nblocks = dbcsr_get_num_blocks(rtbse_env%rho_dbcsr) + CALL dbcsr_iterator_start(iterator_matrix, rtbse_env%rho_dbcsr) + DO n = 1, nblocks + CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, & + row_offset=row_offset, col_offset=col_offset, row_size=row_size, col_size=col_size) + ! Now we have a block corresponding to a single atom pair + j_n_AO = bs_env%sizes_AO(ind_1) + k_n_AO = bs_env%sizes_AO(ind_2) + ! For each atom pair, we need to get contributions from all RI atoms + ! The allocations are as follows + ! Now, get the relevant 3c integrals + ri_offset = 0 + DO ind_i = 1, bs_env%n_atom + i_n_RI = bs_env%sizes_RI(ind_i) + !$OMP PARALLEL DO DEFAULT(none) PRIVATE(i,j,k) & + !$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) + DO i = 1, i_n_RI + DO j = 1, j_n_AO + DO k = 1, k_n_AO + ! Different OMP threads write to different places in memory - should be safe from data race + Qvector(ri_offset + i) = Qvector(ri_offset + i) + int_3c(j + bs_env%i_ao_start_from_atom(ind_1) - 1, & + k + bs_env%i_ao_start_from_atom(ind_2) - 1, & + i + bs_env%i_RI_start_from_atom(ind_i) - 1) & + *block_matrix(j, k) + END DO + END DO + END DO + !$OMP END PARALLEL DO + ri_offset = ri_offset + i_n_RI + END DO + END DO + CALL dbcsr_iterator_stop(iterator_matrix) + ! Now, each rank has contributions from D_jk within its scope + ! Need to sum over different ranks to get the total vector on all ranks + CALL para_env%sum(Qvector) + ! Once this is done, Pvector is current on all ranks + ! Continue with V_PQ summation + nblocks = dbcsr_get_num_blocks(rtbse_env%v_dbcsr) + CALL dbcsr_iterator_start(iterator_matrix, rtbse_env%v_dbcsr) + DO n = 1, nblocks + ! TODO : Try OMP parallelisation over different blocks - expect many more available speedup for large systems + CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, & + row_offset=row_offset, col_offset=col_offset, row_size=row_size, col_size=col_size) + ! TODO : Better names for RI + j_n_AO = bs_env%sizes_RI(ind_1) + k_n_AO = bs_env%sizes_RI(ind_2) + ! The allocations are as follows + !$OMP PARALLEL DO DEFAULT(none) PRIVATE(j,k) & + !$OMP SHARED(block_matrix, Pvector, Qvector,j_n_AO,k_n_AO,row_offset,col_offset) + DO j = 1, j_n_AO + DO k = 1, k_n_AO + Pvector(j + row_offset - 1) = Pvector(j + row_offset - 1) + block_matrix(j, k)*Qvector(k + col_offset - 1) + END DO + END DO + !$OMP END PARALLEL DO + END DO + CALL dbcsr_iterator_stop(iterator_matrix) + ! Again, make sure that the P vector is present on all ranks + CALL para_env%sum(Pvector) + ! Now, for the final trick, iterate over local blocks of v_dbcsr_ao to get the Hartree as dbcsr, then convert to fm + CALL dbcsr_create(v_dbcsr_ao, "Hartree ao", rtbse_env%rho_dbcsr) + CALL copy_fm_to_dbcsr(v_ao, v_dbcsr_ao) + nblocks = dbcsr_get_num_blocks(v_dbcsr_ao) + CALL dbcsr_iterator_start(iterator_matrix, v_dbcsr_ao) + DO n = 1, nblocks + CALL dbcsr_iterator_next_block(iterator_matrix, ind_1, ind_2, block_matrix, & + row_offset=row_offset, col_offset=col_offset) + j_n_AO = bs_env%sizes_AO(ind_1) + k_n_AO = bs_env%sizes_AO(ind_2) + ri_offset = 0 + block_matrix(:, :) = 0.0_dp + DO ind_i = 1, bs_env%n_atom + i_n_RI = bs_env%sizes_RI(ind_i) + ! TODO : In principle, can parallelize over both directions here + !$OMP PARALLEL DO DEFAULT(none) PRIVATE(i,j,k) & + !$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) + DO j = 1, j_n_AO + DO k = 1, k_n_AO + ! Add all the summed up values + DO i = 1, i_n_RI + block_matrix(j, k) = block_matrix(j, k) + int_3c(j + bs_env%i_ao_start_from_atom(ind_1) - 1, & + k + bs_env%i_ao_start_from_atom(ind_2) - 1, & + i + bs_env%i_RI_start_from_atom(ind_i) - 1)*Pvector(i + ri_offset) + END DO + END DO + END DO + !$OMP END PARALLEL DO + ri_offset = ri_offset + i_n_RI + END DO + END DO + CALL dbcsr_iterator_stop(iterator_matrix) + ! Since P vector was present on all the ranks, v_dbcsr_ao has the complete Hartree result + CALL copy_dbcsr_to_fm(v_dbcsr_ao, v_ao) + CALL dbcsr_release(v_dbcsr_ao) + DEALLOCATE (Qvector) + DEALLOCATE (Pvector) + + CALL timestop(handle) + END SUBROUTINE get_hartree_local +! ************************************************************************************************** +!> \brief Calculates the exponential of a matrix in a generalized eigenvalue problem. Specifically, +!> it assumes we have a Hermitian matrix A in the eigenvalue problem AX = BXE, where B is some overlap +!> matrix and E is a diagonal matrix of real eigenvalues. Then, it calculates +!> exp(B^(-1) A) = X exp(E) X^C B +!> \param amatrix Matrix to exponentiate +!> \param bmatrix Overlap matrix +!> \param exponential Exponential exp(B^(-1) A) is stored here after the routine is finished +!> \param eig_scale_opt Optionally scale eigenvalues by a complex number before exponentiating them +!> \param work_opt Optionally provide workspace (of size at least 4) that is used in the calculation +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + SUBROUTINE cp_cfm_gexp(amatrix, bmatrix, exponential, eig_scale_opt, work_opt) + ! TODO : Do interface for real matrices + TYPE(cp_cfm_type), INTENT(IN) :: amatrix + TYPE(cp_cfm_type), INTENT(IN) :: bmatrix + TYPE(cp_cfm_type) :: exponential + COMPLEX(kind=dp), INTENT(IN), OPTIONAL :: eig_scale_opt + TYPE(cp_cfm_type), DIMENSION(:), POINTER, OPTIONAL :: work_opt + CHARACTER(len=*), PARAMETER :: routineN = "cp_cfm_gexp" + COMPLEX(kind=dp) :: eig_scale + REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvalues + COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: expvalues + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: work + LOGICAL :: deallocate_work + INTEGER :: nrow, i, handle + + CALL timeset(routineN, handle) + + ! Argument parsing and sanity checks + IF (PRESENT(eig_scale_opt)) THEN + eig_scale = eig_scale_opt + ELSE + eig_scale = CMPLX(1.0, 0.0, kind=dp) + END IF + + NULLIFY (work) + IF (PRESENT(work_opt) .AND. SIZE(work_opt) >= 4) THEN + deallocate_work = .FALSE. + work => work_opt + ELSE + deallocate_work = .TRUE. + ALLOCATE (work(4)) + ! Allocate the work storage on the fly + DO i = 1, 4 + CALL cp_cfm_create(work(i), amatrix%matrix_struct) + END DO + END IF + + nrow = amatrix%matrix_struct%nrow_global + + ALLOCATE (eigenvalues(nrow)) + ALLOCATE (expvalues(nrow)) + + ! Do not change the amatrix and bmatrix - need to copy them first + CALL cp_cfm_to_cfm(amatrix, work(1)) + CALL cp_cfm_to_cfm(bmatrix, work(2)) + + ! Solve the generalized eigenvalue equation + CALL cp_cfm_geeig(work(1), work(2), work(3), eigenvalues, work(4)) + + ! Scale and exponentiate the eigenvalues + expvalues(:) = EXP(eigenvalues(:)*eig_scale) + + ! Copy eigenvectors to column scale them + CALL cp_cfm_to_cfm(work(3), work(1)) + ! X * exp(E) + CALL cp_cfm_column_scale(work(1), expvalues) + + ! Carry out the remaining operations + ! X * exp(E) * X^C + CALL parallel_gemm("N", "C", nrow, nrow, nrow, & + CMPLX(1.0, 0.0, kind=dp), work(1), work(3), & + CMPLX(0.0, 0.0, kind=dp), work(2)) + ! X * exp(E) * X^C * B + CALL parallel_gemm("N", "N", nrow, nrow, nrow, & + CMPLX(1.0, 0.0, kind=dp), work(2), bmatrix, & + CMPLX(0.0, 0.0, kind=dp), exponential) + + ! Deallocate work storage if necessary + IF (deallocate_work) THEN + DO i = 1, 4 + CALL cp_cfm_release(work(i)) + END DO + DEALLOCATE (work) + END IF + + DEALLOCATE (eigenvalues) + DEALLOCATE (expvalues) + + CALL timestop(handle) + END SUBROUTINE cp_cfm_gexp +END MODULE rt_bse diff --git a/src/emd/rt_bse_io.F b/src/emd/rt_bse_io.F new file mode 100644 index 0000000000..20983df73f --- /dev/null +++ b/src/emd/rt_bse_io.F @@ -0,0 +1,803 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Input/output from the propagation via RT-BSE method. +!> \author Stepan Marek (08.24) +! ************************************************************************************************** + +MODULE rt_bse_io + USE cp_fm_types, ONLY: cp_fm_type, & + cp_fm_p_type, & + cp_fm_create, & + cp_fm_release, & + cp_fm_read_unformatted, & + cp_fm_write_unformatted, & + cp_fm_write_formatted + USE cp_cfm_types, ONLY: cp_cfm_type, & + cp_cfm_p_type, & + cp_fm_to_cfm, & + cp_cfm_to_fm + USE kinds, ONLY: dp, & + default_path_length + USE cp_fm_basic_linalg, ONLY: cp_fm_trace, & + cp_fm_transpose, & + cp_fm_norm + USE parallel_gemm_api, ONLY: parallel_gemm + USE cp_log_handling, ONLY: cp_logger_type, & + cp_get_default_logger + USE cp_output_handling, ONLY: cp_print_key_unit_nr, & + cp_print_key_finished_output, & + cp_print_key_generate_filename, & + low_print_level, & + medium_print_level + USE input_section_types, ONLY: section_vals_type + USE rt_bse_types, ONLY: rtbse_env_type, & + multiply_cfm_fm, & + multiply_fm_cfm + USE cp_files, ONLY: open_file, & + file_exists, & + close_file + USE input_constants, ONLY: do_exact, & + do_bch, & + rtp_bse_ham_g0w0, & + rtp_bse_ham_ks, & + use_rt_restart + USE physcon, ONLY: femtoseconds, & + evolt + USE mathconstants, ONLY: twopi + +#include "../base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse_io" + + #:include "rt_bse_macros.fypp" + + PUBLIC :: output_moments, & + read_moments, & + output_moments_ft, & + output_polarizability, & + output_field, & + read_field, & + output_mos_contravariant, & + output_mos_covariant, & + output_restart, & + read_restart, & + print_etrs_info_header, & + print_etrs_info, & + print_timestep_info, & + print_rtbse_header_info + +CONTAINS + +! ************************************************************************************************** +!> \brief Writes the header and basic info to the standard output +!> \param rtbse_env Entry point - rtbse environment +! ************************************************************************************************** + SUBROUTINE print_rtbse_header_info(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + + WRITE (rtbse_env%unit_nr, *) '' + WRITE (rtbse_env%unit_nr, '(A)') ' /-----------------------------------------------'// & + '------------------------------\' + WRITE (rtbse_env%unit_nr, '(A)') ' | '// & + ' |' + WRITE (rtbse_env%unit_nr, '(A)') ' | Real Time Bethe-Salpeter Propagation'// & + ' |' + WRITE (rtbse_env%unit_nr, '(A)') ' | '// & + ' |' + WRITE (rtbse_env%unit_nr, '(A)') ' \-----------------------------------------------'// & + '------------------------------/' + WRITE (rtbse_env%unit_nr, *) '' + + ! Methods used + WRITE (rtbse_env%unit_nr, '(A19)', advance="no") ' Exponential method' + SELECT CASE (rtbse_env%mat_exp_method) + CASE (do_bch) + WRITE (rtbse_env%unit_nr, '(A61)') 'BCH' + CASE (do_exact) + WRITE (rtbse_env%unit_nr, '(A61)') 'EXACT' + END SELECT + + WRITE (rtbse_env%unit_nr, '(A22)', advance="no") ' Reference Hamiltonian' + SELECT CASE (rtbse_env%ham_reference_type) + CASE (rtp_bse_ham_g0w0) + WRITE (rtbse_env%unit_nr, '(A58)') 'G0W0' + CASE (rtp_bse_ham_ks) + WRITE (rtbse_env%unit_nr, '(A58)') 'Kohn-Sham' + END SELECT + + WRITE (rtbse_env%unit_nr, '(A18,L62)') ' Apply delta pulse', & + rtbse_env%dft_control%rtp_control%apply_delta_pulse + + WRITE (rtbse_env%unit_nr, '(A)') '' + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Writes the update after single etrs iteration - only for log level > medium +!> \param rtbse_env Entry point - rtbse environment +! ************************************************************************************************** + SUBROUTINE print_etrs_info(rtbse_env, step, metric) + TYPE(rtbse_env_type) :: rtbse_env + INTEGER :: step + REAL(kind=dp) :: metric + TYPE(cp_logger_type), POINTER :: logger + + logger => cp_get_default_logger() + + IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN + WRITE (rtbse_env%unit_nr, '(A7,I5, E20.8E3)') ' RTBSE|', step, metric + END IF + + END SUBROUTINE +! ************************************************************************************************** +!> \brief Writes the header for the etrs iteration updates - only for log level > medium +!> \param rtbse_env Entry point - rtbse environment +! ************************************************************************************************** + SUBROUTINE print_etrs_info_header(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + + logger => cp_get_default_logger() + + IF (logger%iter_info%print_level > medium_print_level .AND. rtbse_env%unit_nr > 0) THEN + WRITE (rtbse_env%unit_nr, '(A13, A20)') ' RTBSE| Iter.', 'Convergence' + END IF + + END SUBROUTINE +! ************************************************************************************************** +!> \brief Writes the header for the etrs iteration updates - only for log level > medium +!> \param rtbse_env Entry point - rtbse environment +! ************************************************************************************************** + SUBROUTINE print_timestep_info(rtbse_env, step, convergence, electron_num_re, etrs_num) + TYPE(rtbse_env_type) :: rtbse_env + INTEGER :: step + REAL(kind=dp) :: convergence + REAL(kind=dp) :: electron_num_re + INTEGER :: etrs_num + TYPE(cp_logger_type), POINTER :: logger + + logger => cp_get_default_logger() + + IF (logger%iter_info%print_level > low_print_level .AND. rtbse_env%unit_nr > 0) THEN + WRITE (rtbse_env%unit_nr, '(A23,A20,A20,A17)') " RTBSE| Simulation step", "Convergence", & + "Electron number", "ETRS Iterations" + WRITE (rtbse_env%unit_nr, '(A7,I16,E20.8E3,E20.8E3,I17)') ' RTBSE|', step, convergence, & + electron_num_re, etrs_num + END IF + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Outputs the matrix in MO basis for matrix coefficients corresponding to contravariant +!> operator, i.e. density matrix +!> \param rtbse_env Entry point - gwbse environment +!> \param rho Density matrix in AO basis +!> \param rtp_section RTP input section +! ************************************************************************************************** + SUBROUTINE output_mos_contravariant(rtbse_env, rho, print_key_section) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho + TYPE(section_vals_type), POINTER :: print_key_section + TYPE(cp_logger_type), POINTER :: logger + INTEGER :: j, rho_unit_re, rho_unit_im + CHARACTER(len=14), DIMENSION(4) :: file_labels + + file_labels(1) = "_SPIN_A_RE.dat" + file_labels(2) = "_SPIN_A_IM.dat" + file_labels(3) = "_SPIN_B_RE.dat" + file_labels(4) = "_SPIN_B_IM.dat" + logger => cp_get_default_logger() + ! Start by multiplying the current density by MOS + DO j = 1, rtbse_env%n_spin + rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1)) + rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j)) + ! Transform the density matrix into molecular orbitals basis and print it out + ! S * rho + CALL multiply_fm_cfm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%S_fm, rho(j), & + 0.0_dp, rtbse_env%rho_workspace(1)) + ! C^T * S * rho + CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), rtbse_env%rho_workspace(1), & + 0.0_dp, rtbse_env%rho_workspace(2)) + ! C^T * S * rho * S + CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%rho_workspace(2), rtbse_env%S_fm, & + 0.0_dp, rtbse_env%rho_workspace(1)) + ! C^T * S * rho * S * C + CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), & + 0.0_dp, rtbse_env%rho_workspace(2)) + ! Print real and imaginary parts separately + CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), & + rtbse_env%real_workspace(1), rtbse_env%real_workspace(2)) + CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re) + CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im) + CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section) + CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section) + END DO + END SUBROUTINE +! ************************************************************************************************** +!> \brief Outputs the matrix in MO basis for matrix components corresponding to covariant representation, +!> i.e. the Hamiltonian matrix +!> \param rtbse_env Entry point - gwbse environment +!> \param cohsex cohsex matrix in AO basis, covariant representation +!> \param rtp_section RTP input section +! ************************************************************************************************** + SUBROUTINE output_mos_covariant(rtbse_env, ham, print_key_section) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham + TYPE(section_vals_type), POINTER :: print_key_section + TYPE(cp_logger_type), POINTER :: logger + INTEGER :: j, rho_unit_re, rho_unit_im + CHARACTER(len=21), DIMENSION(4) :: file_labels + + file_labels(1) = "_SPIN_A_RE.dat" + file_labels(2) = "_SPIN_A_IM.dat" + file_labels(3) = "_SPIN_B_RE.dat" + file_labels(4) = "_SPIN_B_IM.dat" + logger => cp_get_default_logger() + DO j = 1, rtbse_env%n_spin + rho_unit_re = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j - 1)) + rho_unit_im = cp_print_key_unit_nr(logger, print_key_section, extension=file_labels(2*j)) + ! C^T * cohsex + CALL multiply_fm_cfm("T", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%bs_env%fm_mo_coeff_Gamma(j), ham(j), & + 0.0_dp, rtbse_env%rho_workspace(1)) + ! C^T * cohsex * C + CALL multiply_cfm_fm("N", "N", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_ao, & + 1.0_dp, rtbse_env%rho_workspace(1), rtbse_env%bs_env%fm_mo_coeff_Gamma(j), & + 0.0_dp, rtbse_env%rho_workspace(2)) + ! Print real and imaginary parts separately + CALL cp_cfm_to_fm(rtbse_env%rho_workspace(2), & + rtbse_env%real_workspace(1), rtbse_env%real_workspace(2)) + CALL cp_fm_write_formatted(rtbse_env%real_workspace(1), rho_unit_re) + CALL cp_fm_write_formatted(rtbse_env%real_workspace(2), rho_unit_im) + CALL cp_print_key_finished_output(rho_unit_re, logger, print_key_section) + CALL cp_print_key_finished_output(rho_unit_im, logger, print_key_section) + END DO + END SUBROUTINE +! ************************************************************************************************** +!> \brief Prints the current field components into a file provided by input +!> \param rtbse_env Entry point - gwbse environment +!> \param rtp_section RTP input section +! ************************************************************************************************** + SUBROUTINE output_field(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + INTEGER :: field_unit, n, i + + ! Get logger + logger => cp_get_default_logger() + ! Get file descriptor + field_unit = cp_print_key_unit_nr(logger, rtbse_env%field_section, extension=".dat") + ! If the file descriptor is non-zero, output field + ! TODO : Output also in SI + IF (field_unit /= -1) THEN + WRITE (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%sim_time*femtoseconds, & + rtbse_env%field(1), rtbse_env%field(2), rtbse_env%field(3) + END IF + ! Write the output to memory for FT + ! Need the absolute index + n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1 + DO i = 1, 3 + rtbse_env%field_trace(i)%series(n) = rtbse_env%field(i) + END DO + rtbse_env%time_trace%series(n) = rtbse_env%sim_time + CALL cp_print_key_finished_output(field_unit, logger, rtbse_env%field_section) + + END SUBROUTINE +! ************************************************************************************************** +!> \brief Reads the field from the files provided by input - useful for the continuation run +!> \param rtbse_env Entry point - gwbse environment +!> \param rtp_section RTP input section +! ************************************************************************************************** + SUBROUTINE read_field(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + CHARACTER(len=default_path_length) :: save_name + INTEGER :: k, n, field_unit + + ! Get logger + logger => cp_get_default_logger() + ! Get file name + save_name = cp_print_key_generate_filename(logger, rtbse_env%field_section, extension=".dat", my_local=.FALSE.) + IF (file_exists(save_name)) THEN + CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", & + unit_number=field_unit) + DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start + n = k - rtbse_env%sim_start_orig + 1 + READ (field_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') rtbse_env%time_trace%series(n), & + rtbse_env%field_trace(1)%series(n), rtbse_env%field_trace(2)%series(n), rtbse_env%field_trace(3)%series(n) + ! Set the time units back to atomic units + rtbse_env%time_trace%series(n) = rtbse_env%time_trace%series(n)/femtoseconds + END DO + CALL close_file(field_unit) + ELSE IF (.NOT. rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. & + rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN + CPWARN("Restart without RT field file - unknown field trace set to zero.") + END IF + END SUBROUTINE read_field + +! ************************************************************************************************** +!> \brief Outputs the expectation value of moments from a given density matrix +!> \note Moments matrix is provided by the rtbse_env, uses rho_workspace(1:3) +!> \param rtbse_env Entry point - gwbse environment +!> \param rho Density matrix in AO basis +!> \param rtp_section RTP section of the input parameters, where moments destination may be present +! ************************************************************************************************** + SUBROUTINE output_moments(rtbse_env, rho) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho + TYPE(cp_logger_type), POINTER :: logger + INTEGER :: i, j, n, moments_unit_re, moments_unit_im + CHARACTER(len=14), DIMENSION(4) :: file_labels + REAL(kind=dp), DIMENSION(3) :: moments + + ! Start by getting the relevant file unit + file_labels(1) = "_SPIN_A_RE.dat" + file_labels(2) = "_SPIN_A_IM.dat" + file_labels(3) = "_SPIN_B_RE.dat" + file_labels(4) = "_SPIN_B_IM.dat" + logger => cp_get_default_logger() + DO j = 1, rtbse_env%n_spin + moments_unit_re = cp_print_key_unit_nr(logger, rtbse_env%moments_section, extension=file_labels(2*j - 1)) + moments_unit_im = cp_print_key_unit_nr(logger, rtbse_env%moments_section, extension=file_labels(2*j)) + ! If, for any reason, the file unit is not provided, skip to next cycle immediately + ! TODO : Specify output units in config + ! Need to transpose due to the definition of trace function + CALL cp_cfm_to_fm(msource=rho(j), mtargetr=rtbse_env%real_workspace(2)) + DO i = 1, 3 + ! Moments should be symmetric, test without transopose? + CALL cp_fm_transpose(rtbse_env%moments(i), rtbse_env%real_workspace(1)) + CALL cp_fm_trace(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), moments(i)) + ! Scale by spin degeneracy and electron charge + moments(i) = -moments(i)*rtbse_env%spin_degeneracy + END DO + ! Output to the file + IF (moments_unit_re > 0) THEN + ! If outputting to standard output, also prepend identifying 'MOMENTS_TRACE|' string + IF (rtbse_env%unit_nr == moments_unit_re) THEN + WRITE (moments_unit_re, '(A18,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + " MOMENTS_TRACE_RE|", rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3) + ELSE + WRITE (moments_unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3) + CALL cp_print_key_finished_output(moments_unit_re, logger, rtbse_env%moments_section) + END IF + END IF + ! Save to memory for FT - real part + n = rtbse_env%sim_step - rtbse_env%sim_start_orig + 1 + DO i = 1, 3 + rtbse_env%moments_trace(i)%series(n) = CMPLX(moments(i), 0.0, kind=dp) + END DO + ! Same for imaginary part + CALL cp_cfm_to_fm(msource=rho(j), mtargeti=rtbse_env%real_workspace(2)) + DO i = 1, 3 + CALL cp_fm_transpose(rtbse_env%moments(i), rtbse_env%real_workspace(1)) + CALL cp_fm_trace(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), moments(i)) + ! Scale by spin degeneracy and electron charge + moments(i) = -moments(i)*rtbse_env%spin_degeneracy + END DO + ! Output to the file + IF (moments_unit_im > 0) THEN + ! If outputting to standard output, also prepend identifying 'MOMENTS_TRACE|' string + IF (rtbse_env%unit_nr == moments_unit_im) THEN + WRITE (moments_unit_im, '(A18,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + " MOMENTS_TRACE_IM|", rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3) + ELSE + WRITE (moments_unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + rtbse_env%sim_time*femtoseconds, moments(1), moments(2), moments(3) + ! Close the files + CALL cp_print_key_finished_output(moments_unit_im, logger, rtbse_env%moments_section) + END IF + END IF + ! Save to memory for FT - imag part + DO i = 1, 3 + rtbse_env%moments_trace(i)%series(n) = rtbse_env%moments_trace(i)%series(n) + CMPLX(0.0, moments(i), kind=dp) + END DO + END DO + END SUBROUTINE +! ************************************************************************************************** +!> \brief Outputs the restart info (last finished iteration step) + restard density matrix +!> \param restart_section Print key section for the restart files +!> \param rho Density matrix in AO basis +!> \param time_index Time index to be written into the info file +! ************************************************************************************************** + SUBROUTINE output_restart(rtbse_env, rho, time_index) + TYPE(rtbse_env_type), POINTER :: rtbse_env + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho + INTEGER :: time_index + TYPE(cp_fm_type), DIMENSION(:), POINTER :: workspace + CHARACTER(len=17), DIMENSION(4) :: file_labels + TYPE(cp_logger_type), POINTER :: logger + INTEGER :: rho_unit_nr, i + + ! Default labels distinguishing up to two spin species and real/imaginary parts + file_labels(1) = "_SPIN_A_RE.matrix" + file_labels(2) = "_SPIN_A_IM.matrix" + file_labels(3) = "_SPIN_B_RE.matrix" + file_labels(4) = "_SPIN_B_IM.matrix" + + logger => cp_get_default_logger() + + workspace => rtbse_env%real_workspace + + DO i = 1, rtbse_env%n_spin + CALL cp_cfm_to_fm(rho(i), workspace(1), workspace(2)) + ! Real part + rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i - 1), & + file_form="UNFORMATTED", file_position="REWIND") + CALL cp_fm_write_unformatted(workspace(1), rho_unit_nr) + CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section) + ! Imag part + rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=file_labels(2*i), & + file_form="UNFORMATTED", file_position="REWIND") + CALL cp_fm_write_unformatted(workspace(2), rho_unit_nr) + CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section) + ! Info + rho_unit_nr = cp_print_key_unit_nr(logger, rtbse_env%restart_section, extension=".info", & + file_form="UNFORMATTED", file_position="REWIND") + IF (rho_unit_nr > 0) WRITE (rho_unit_nr) time_index + CALL cp_print_key_finished_output(rho_unit_nr, logger, rtbse_env%restart_section) + END DO + END SUBROUTINE output_restart +! ************************************************************************************************** +!> \brief Reads the density matrix from restart files and updates the starting time +!> \param restart_section Print key section for the restart files +!> \param rho Density matrix in AO basis +!> \param time_index Time index to be written into the info file +! ************************************************************************************************** + SUBROUTINE read_restart(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + CHARACTER(len=default_path_length) :: save_name, save_name_2 + INTEGER :: rho_unit_nr, j + CHARACTER(len=17), DIMENSION(4) :: file_labels + + ! This allows the delta kick and output of moment at time 0 in all cases + ! except the case when both imaginary and real parts of the density are read + rtbse_env%restart_extracted = .FALSE. + logger => cp_get_default_logger() + ! Start by probing/loading info file + save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, extension=".info", my_local=.FALSE.) + IF (file_exists(save_name)) THEN + CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", & + unit_number=rho_unit_nr) + READ (rho_unit_nr) rtbse_env%sim_start + CALL close_file(rho_unit_nr) + IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A31,I25,A24)') " RTBSE| Starting from timestep ", & + rtbse_env%sim_start, ", delta kick NOT applied" + ELSE + CPWARN("Restart required but no info file found - starting from sim_step given in input") + END IF + + ! Default labels distinguishing up to two spin species and real/imaginary parts + file_labels(1) = "_SPIN_A_RE.matrix" + file_labels(2) = "_SPIN_A_IM.matrix" + file_labels(3) = "_SPIN_B_RE.matrix" + file_labels(4) = "_SPIN_B_IM.matrix" + DO j = 1, rtbse_env%n_spin + save_name = cp_print_key_generate_filename(logger, rtbse_env%restart_section, & + extension=file_labels(2*j - 1), my_local=.FALSE.) + save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%restart_section, & + extension=file_labels(2*j), my_local=.FALSE.) + IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN + CALL open_file(save_name, file_status="OLD", file_form="UNFORMATTED", file_action="READ", & + unit_number=rho_unit_nr) + CALL cp_fm_read_unformatted(rtbse_env%real_workspace(1), rho_unit_nr) + CALL close_file(rho_unit_nr) + CALL open_file(save_name_2, file_status="OLD", file_form="UNFORMATTED", file_action="READ", & + unit_number=rho_unit_nr) + CALL cp_fm_read_unformatted(rtbse_env%real_workspace(2), rho_unit_nr) + CALL close_file(rho_unit_nr) + CALL cp_fm_to_cfm(rtbse_env%real_workspace(1), rtbse_env%real_workspace(2), & + rtbse_env%rho(j)) + rtbse_env%restart_extracted = .TRUE. + ELSE + CPWARN("Restart without some restart matrices - starting from SCF density.") + END IF + END DO + END SUBROUTINE read_restart +! ************************************************************************************************** +!> \brief Reads the moments and time traces from the save files +!> \param rtbse_env GW-BSE environment (assumes consistent setup, i.e. a continuation run). +!> Otherwise, the traces are set at zero +! ************************************************************************************************** + SUBROUTINE read_moments(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + CHARACTER(len=default_path_length) :: save_name, save_name_2 + INTEGER :: i, j, k, moments_unit_re, moments_unit_im, n + CHARACTER(len=17), DIMENSION(4) :: file_labels + REAL(kind=dp), DIMENSION(3) :: moments_re, moments_im + REAL(kind=dp) :: timestamp + + logger => cp_get_default_logger() + + ! Read moments from the previous run + ! Default labels distinguishing up to two spin species and real/imaginary parts + file_labels(1) = "_SPIN_A_RE.dat" + file_labels(2) = "_SPIN_A_IM.dat" + file_labels(3) = "_SPIN_B_RE.dat" + file_labels(4) = "_SPIN_B_IM.dat" + DO j = 1, rtbse_env%n_spin + save_name = cp_print_key_generate_filename(logger, rtbse_env%moments_section, & + extension=file_labels(2*j - 1), my_local=.FALSE.) + save_name_2 = cp_print_key_generate_filename(logger, rtbse_env%moments_section, & + extension=file_labels(2*j), my_local=.FALSE.) + IF (file_exists(save_name) .AND. file_exists(save_name_2)) THEN + CALL open_file(save_name, file_status="OLD", file_form="FORMATTED", file_action="READ", & + unit_number=moments_unit_re) + CALL open_file(save_name_2, file_status="OLD", file_form="FORMATTED", file_action="READ", & + unit_number=moments_unit_im) + ! Extra time step for the initial one + DO k = rtbse_env%sim_start_orig, rtbse_env%sim_start + ! Determine the absolute time step - offset in memory + n = k - rtbse_env%sim_start_orig + 1 + READ (moments_unit_re, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') timestamp, & + moments_re(1), moments_re(2), moments_re(3) + READ (moments_unit_im, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') timestamp, & + moments_im(1), moments_im(2), moments_im(3) + DO i = 1, 3 + rtbse_env%moments_trace(i)%series(n) = CMPLX(moments_re(i), moments_im(i), kind=dp) + END DO + END DO + ! Change back to atomic units in the trace + rtbse_env%time_trace%series(:) = rtbse_env%time_trace%series(:)/femtoseconds + CALL close_file(moments_unit_re) + CALL close_file(moments_unit_im) + ELSE IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN + CPWARN("Restart without previous moments - unknown moments trace set to zero.") + END IF + END DO + END SUBROUTINE read_moments +! ************************************************************************************************** +!> \brief Outputs the Fourier transform of moments stored in the environment memory to the configured file +!> \param rtbse_env GW-BSE environment +! ************************************************************************************************** + SUBROUTINE output_moments_ft(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + REAL(kind=dp), DIMENSION(:), POINTER :: omega_series, & + ft_real_series, & + ft_imag_series, & + value_series + ! Stores the data in ready for output format + ! - first dimension is 6 - 1 - real part along x, 2 - imag part along x, 3 - real part along y, ... + REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: ft_full_series + INTEGER :: i, n, ft_unit + + logger => cp_get_default_logger() + + n = rtbse_env%sim_nsteps + 2 + NULLIFY (omega_series) + ALLOCATE (omega_series(n), source=0.0_dp) + NULLIFY (ft_real_series) + ALLOCATE (ft_real_series(n), source=0.0_dp) + NULLIFY (ft_imag_series) + ALLOCATE (ft_imag_series(n), source=0.0_dp) + NULLIFY (value_series) + ALLOCATE (value_series(n), source=0.0_dp) + ALLOCATE (ft_full_series(6, n)) + ! Carry out for each direction independently and real and imaginary parts also independently + DO i = 1, 3 + ! Real part of the value first + value_series(:) = REAL(rtbse_env%moments_trace(i)%series(:)) + CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, & + rtbse_env%ft_damping, rtbse_env%ft_start) + ft_full_series(2*i - 1, :) = ft_real_series(:) + ft_full_series(2*i, :) = ft_imag_series(:) + ! Now imaginary part + value_series(:) = AIMAG(rtbse_env%moments_trace(i)%series(:)) + CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, & + rtbse_env%ft_damping, rtbse_env%ft_start) + ft_full_series(2*i - 1, :) = ft_full_series(2*i - 1, :) - ft_imag_series + ft_full_series(2*i, :) = ft_full_series(2*i, :) + ft_real_series + END DO + DEALLOCATE (ft_real_series) + DEALLOCATE (ft_imag_series) + DEALLOCATE (value_series) + ! Now, write these to file + ft_unit = cp_print_key_unit_nr(logger, rtbse_env%ft_section, extension=".dat", & + file_form="FORMATTED", file_position="REWIND") + IF (ft_unit > 0) THEN + DO i = 1, n + WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + omega_series(i), ft_full_series(1, i), ft_full_series(2, i), & + ft_full_series(3, i), ft_full_series(4, i), & + ft_full_series(5, i), ft_full_series(6, i) + END DO + END IF + CALL cp_print_key_finished_output(ft_unit, logger, rtbse_env%ft_section) + DEALLOCATE (omega_series) + DEALLOCATE (ft_full_series) + END SUBROUTINE output_moments_ft +! ************************************************************************************************** +!> \brief Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega), +!> where i and j are provided by the configuration. The tensor element is energy dependent and +!> has real and imaginary parts +!> \param rtbse_env GW-BSE environment +! ************************************************************************************************** + SUBROUTINE output_polarizability(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + TYPE(cp_logger_type), POINTER :: logger + REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_series, & + ft_real_series, & + ft_imag_series, & + value_series + COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: moment_series, & + field_series, & + polarizability_series + INTEGER :: pol_unit, & + i, n + + logger => cp_get_default_logger() + + n = rtbse_env%sim_nsteps + 2 + ! All allocations together, although could save some memory, if required by consequent deallocations + ALLOCATE (omega_series(n), source=0.0_dp) + ALLOCATE (ft_real_series(n), source=0.0_dp) + ALLOCATE (ft_imag_series(n), source=0.0_dp) + ALLOCATE (value_series(n), source=0.0_dp) + ALLOCATE (moment_series(n), source=CMPLX(0.0, 0.0, kind=dp)) + ALLOCATE (field_series(n), source=CMPLX(0.0, 0.0, kind=dp)) + ALLOCATE (polarizability_series(n), source=CMPLX(0.0, 0.0, kind=dp)) + + ! The moment ft + ! Real part + value_series(:) = REAL(rtbse_env%moments_trace(rtbse_env%pol_element(1))%series(:)) + CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, & + rtbse_env%ft_damping, rtbse_env%ft_start) + moment_series(:) = CMPLX(ft_real_series(:), ft_imag_series(:), kind=dp) + ! Imaginary part + value_series(:) = AIMAG(rtbse_env%moments_trace(rtbse_env%pol_element(1))%series(:)) + CALL ft_simple(rtbse_env%time_trace%series, value_series, omega_series, ft_real_series, ft_imag_series, & + rtbse_env%ft_damping, rtbse_env%ft_start) + moment_series(:) = moment_series(:) + CMPLX(-ft_imag_series(:), ft_real_series(:), kind=dp) + + ! Calculate the field transform - store it in ft_real_series + IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse) THEN + ! Only divide by constant magnitude in atomic units + field_series(:) = CMPLX(rtbse_env%dft_control%rtp_control%delta_pulse_scale, 0.0, kind=dp) + ELSE + ! Calculate the transform of the field as well and divide by it + ! The field FT is not damped - assume field is localised in time + ! The field is strictly real + CALL ft_simple(rtbse_env%time_trace%series, rtbse_env%field_trace(rtbse_env%pol_element(2))%series, & + omega_series, ft_real_series, ft_imag_series, 0.0_dp, rtbse_env%ft_start) + ! Regularization for the case when ft_series becomes identically zero - TODO : Set in config + field_series(:) = CMPLX(ft_real_series(:), ft_imag_series(:) + 1.0e-20, kind=dp) + END IF + ! Divide to get the polarizability series + ! Real part + polarizability_series(:) = moment_series(:)/field_series(:) + + ! Change units to eV for energy + ! use value_series for energy and moment_series for polarizability + value_series(:) = omega_series(:)*evolt + ! Use below for printing of field FT in case of problems + ! PRINT *, "=== Field FT" + ! DO i=1,n + ! PRINT '(E20.8E3,E20.8E3,E20.8E3)', value_series(i), REAL(field_series(i)), AIMAG(field_series(i)) + ! END DO + ! PRINT *, "=== Field FT" + ! Print out the polarizability to a file + pol_unit = cp_print_key_unit_nr(logger, rtbse_env%pol_section, extension=".dat", & + file_form="FORMATTED", file_position="REWIND") + IF (pol_unit > 0) THEN + IF (pol_unit == rtbse_env%unit_nr) THEN + WRITE (pol_unit, '(A16,A20,A22,A22)') & + " POLARIZABILITY|", "Energy [eV]", "Real polarizability", "Imag polarizability" + DO i = 1, n + WRITE (pol_unit, '(A16,E20.8E3,E22.8E3,E22.8E3)') & + " POLARIZABILITY|", value_series(i), REAL(polarizability_series(i)), AIMAG(polarizability_series(i)) + END DO + ELSE + WRITE (pol_unit, '(A1,A19,A20,A20,A20)') & + "#", "omega [a.u.]", "Energy [eV]", "Real polarizability", "Imag polarizability" + DO i = 1, n + WRITE (pol_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + omega_series(i), value_series(i), REAL(polarizability_series(i)), AIMAG(polarizability_series(i)) + END DO + CALL cp_print_key_finished_output(pol_unit, logger, rtbse_env%ft_section) + END IF + END IF + + DEALLOCATE (value_series) + DEALLOCATE (ft_real_series) + DEALLOCATE (ft_imag_series) + + DEALLOCATE (field_series) + DEALLOCATE (moment_series) + + DEALLOCATE (omega_series) + DEALLOCATE (polarizability_series) + END SUBROUTINE output_polarizability +! ************************************************************************************************** +!> \brief Naively calculates the Fourier transform - it is not the bottleneck of this calculation +!> \param time_series Timestamps in atomic units of time +!> \param value_series Values to be Fourier transformed - moments, field etc. So far only real. +!> \param omega_series Array to be filled by sampled values of frequency +!> \param result_series FT of the value series - real values (cosines) +!> \param iresult_series FT of the value series - imaginary values (sines) +!> \param damping_opt Supply custom exponential damping - default is 4.0/totalTime, i.e. ratio +!> of last and first element in windowed value series is reduced by e^(-4) +!> \param t0_opt Carry the FT only starting from certain time - allows for exclusion of trace before +!> the pulse application etc. +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + ! So far only for defined one dimensional series + SUBROUTINE ft_simple(time_series, value_series, omega_series, result_series, iresult_series, & + damping_opt, t0_opt) + REAL(kind=dp), DIMENSION(:) :: time_series + REAL(kind=dp), DIMENSION(:) :: value_series + REAL(kind=dp), DIMENSION(:) :: omega_series + REAL(kind=dp), DIMENSION(:) :: result_series + REAL(kind=dp), DIMENSION(:) :: iresult_series + REAL(kind=dp), OPTIONAL :: damping_opt + REAL(kind=dp), OPTIONAL :: t0_opt + CHARACTER(len=*), PARAMETER :: routineN = "ft_simple" + INTEGER :: N, i, j, t0_i, j_wrap, handle + REAL(kind=dp) :: t0, delta_t, delta_omega, damping + + CALL timeset(routineN, handle) + + N = SIZE(time_series) + + t0_i = 1 + IF (PRESENT(t0_opt)) THEN + ! Determine the index at which we start applying the damping + ! Also the index around which we fold around + DO i = 1, N + ! Increase until we break or reach the end of the time series + t0_i = i + IF (time_series(i) >= t0_opt) THEN + EXIT + END IF + END DO + END IF + + t0 = time_series(t0_i) + + ! Default damping so that at the end of the time series, divide value by e^-4 + damping = 4.0_dp/(time_series(N) - time_series(t0_i)) + IF (PRESENT(damping_opt)) THEN + IF (damping_opt >= 0.0_dp) damping = damping_opt + END IF + + ! Construct the grid + ! Assume series equidistant + delta_t = time_series(2) - time_series(1) + delta_omega = twopi/(time_series(N) - time_series(1)) + DO i = 1, N + result_series(i) = 0.0_dp + iresult_series(i) = 0.0_dp + omega_series(i) = delta_omega*(i - 1) + DO j = 1, N + j_wrap = MODULO(j + t0_i - 2, N) + 1 + result_series(i) = result_series(i) + COS(twopi*(i - 1)*(j - 1)/N)* & + EXP(-damping*delta_t*(j - 1))*value_series(j_wrap) + iresult_series(i) = iresult_series(i) + SIN(twopi*(i - 1)*(j - 1)/N)* & + EXP(-damping*delta_t*(j - 1))*value_series(j_wrap) + END DO + END DO + delta_omega = twopi/(time_series(N) - time_series(1)) + result_series(:) = delta_t*result_series(:) + iresult_series(:) = delta_t*iresult_series(:) + + CALL timestop(handle) + + END SUBROUTINE +END MODULE rt_bse_io diff --git a/src/emd/rt_bse_macros.fypp b/src/emd/rt_bse_macros.fypp new file mode 100644 index 0000000000..4a02e474d6 --- /dev/null +++ b/src/emd/rt_bse_macros.fypp @@ -0,0 +1,22 @@ +#!-------------------------------------------------------------------------------------------------! +#! CP2K: A general program to perform molecular dynamics simulations ! +#! Copyright 2000-2024 CP2K developers group ! +#! ! +#! SPDX-License-Identifier: GPL-2.0-or-later ! +#!-------------------------------------------------------------------------------------------------! + +#! ************************************************************************************************** +#!> \brief Macros for propagation via RT-BSE method +#!> \note Contains the SPIN_DO construct +#!> \author Stepan Marek (12.23) +#! ************************************************************************************************** + +#:mute + + #:def SPIN_DO(i, re, im, nspin) + DO ${i}$ = 1, ${nspin}$ + ${re}$ = 2*${i}$-1 + ${im}$ = 2*${i}$ + #:enddef + +#:endmute diff --git a/src/emd/rt_bse_types.F b/src/emd/rt_bse_types.F new file mode 100644 index 0000000000..ac8f26bbfb --- /dev/null +++ b/src/emd/rt_bse_types.F @@ -0,0 +1,755 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Data storage and other types for propagation via RT-BSE method. +!> \author Stepan Marek (01.24) +! ************************************************************************************************** + +MODULE rt_bse_types + + USE kinds, ONLY: dp + USE cp_fm_types, ONLY: cp_fm_type, & + cp_fm_p_type, & + cp_fm_release, & + cp_fm_create, & + cp_fm_set_all, & + cp_fm_write_formatted, & + cp_fm_to_fm + USE cp_cfm_types, ONLY: cp_cfm_p_type, & + cp_cfm_type, & + cp_cfm_set_all, & + cp_cfm_create, & + cp_fm_to_cfm, & + cp_cfm_to_fm, & + cp_cfm_to_cfm, & + cp_cfm_release + USE cp_fm_basic_linalg, ONLY: cp_fm_invert, & + cp_fm_transpose, & + cp_fm_column_scale, & + cp_fm_scale_and_add + USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale + USE cp_dbcsr_api, ONLY: dbcsr_type, & + dbcsr_p_type, & + dbcsr_create, & + dbcsr_release, & + dbcsr_print, & + dbcsr_get_info, & + dbcsr_copy, & + dbcsr_set, & + dbcsr_scale, & + dbcsr_type_complex_8 + USE parallel_gemm_api, ONLY: parallel_gemm + USE dbt_api, ONLY: dbt_type, & + dbt_pgrid_type, & + dbt_pgrid_create, & + dbt_pgrid_destroy, & + dbt_mp_environ_pgrid, & + dbt_default_distvec, & + dbt_distribution_type, & + dbt_distribution_new, & + dbt_distribution_destroy, & + dbt_create, & + dbt_copy_matrix_to_tensor, & + dbt_get_num_blocks, & + dbt_destroy + USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm, & + copy_fm_to_dbcsr + USE qs_mo_types, ONLY: mo_set_type + USE cp_control_types, ONLY: dft_control_type + USE qs_environment_types, ONLY: qs_environment_type, & + get_qs_env + USE force_env_types, ONLY: force_env_type + USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type + USE post_scf_bandstructure_utils, ONLY: create_and_init_bs_env + USE rt_propagation_types, ONLY: rt_prop_type, & + get_rtp + USE rt_propagator_init, ONLY: rt_initialize_rho_from_mos + USE rt_propagation_methods, ONLY: s_matrices_create + USE qs_moments, ONLY: build_local_moment_matrix + USE moments_utils, ONLY: get_reference_point + USE matrix_exp, ONLY: get_nsquare_norder + USE gw_integrals, ONLY: build_3c_integral_block + USE gw_large_cell_Gamma, ONLY: compute_3c_integrals + USE qs_tensors, ONLY: neighbor_list_3c_destroy +! USE gw_utils, ONLY: create_and_init_bs_env_for_gw,& +! setup_AO_and_RI_basis_set,& +! get_RI_basis_and_basis_function_indices,& +! set_heuristic_parameters,& +! set_parallelization_parameters,& +! allocate_and_fill_matrices_and_arrays,& +! create_tensors + USE libint_wrapper, ONLY: cp_libint_static_init + USE input_constants, ONLY: use_mom_ref_coac, & + use_mom_ref_zero, & + use_mom_ref_user, & + use_mom_ref_com, & + rtp_bse_ham_g0w0, & + rtp_bse_ham_ks, & + do_taylor, & + do_bch, & + do_exact + USE physcon, ONLY: angstrom + USE mathconstants, ONLY: z_zero + USE input_section_types, ONLY: section_vals_type, & + section_vals_val_get, & + section_vals_get_subs_vals + +#include "../base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse" + + #:include "rt_bse_macros.fypp" + + PUBLIC :: rtbse_env_type, & + create_rtbse_env, & + release_rtbse_env, & + multiply_cfm_fm, & + multiply_fm_cfm + + ! Created so that we can have an array of pointers to arrays + TYPE series_real_type + REAL(kind=dp), DIMENSION(:), POINTER :: series => NULL() + END TYPE series_real_type + TYPE series_complex_type + COMPLEX(kind=dp), DIMENSION(:), POINTER :: series => NULL() + END TYPE series_complex_type + +! ************************************************************************************************** +!> \param n_spin Number of spin channels that are present +!> \param n_ao Number of atomic orbitals +!> \param n_RI Number of RI orbitals +!> \param n_occ Number of occupied orbitals, spin dependent +!> \param spin_degeneracy Number of electrons per orbital +!> \param field Electric field calculated at the given timestep +!> \param moments Moment operators along cartesian directions - centered at zero charge - used for plotting +!> \param moments_field Moment operators along cartesian directions - used to coupling to the field - +!> origin bound to unit cell +!> \param sim_step Current step of the simulation +!> \param sim_start Starting step of the simulation +!> \param sim_nsteps Number of steps of the simulation +!> \param sim_time Current time of the simulation +!> \param sim_dt Timestep of the simulation +!> \param etrs_threshold Self-consistency threshold for enforced time reversal symmetry propagation +!> \param exp_accuracy Threshold for matrix exponential calculation +!> \param dft_control DFT control parameters +!> \param ham_effective Real and imaginary part of the effective Hamiltonian used to propagate +!> the density matrix +!> \param ham_reference Reference Hamiltonian, which does not change in the propagation = DFT+G0W0 - initial Hartree - initial COHSEX +!> \param ham_workspace Workspace matrices for use with the Hamiltonian propagation - storage of +!> exponential propagators etc. +!> \param rho Density matrix at the current time step +!> \param rho_new Density matrix - workspace in ETRS +!> \param rho_last Density matrix - workspace in ETRS +!> \param rho_new_last Density matrix - workspace in ETRS +!> \param rho_M Density matrix - workspace in ETRS +!> \param S_inv_fm Inverse overlap matrix, full matrix +!> \param S_fm Overlap matrix, full matrix +!> \param S_inv Inverse overlap matrix, sparse matrix +!> \param rho_dbcsr Density matrix, sparse matrix +!> \param rho_workspace Matrices for storage of density matrix at different timesteps for +!> interpolation and self-consistency checks etc. +!> \param complex_workspace Workspace for complex density (exact diagonalisation) +!> \param complex_s Complex overlap matrix (exact diagonalisation) +!> \param real_eigvals Eigenvalues of hermitian matrix (exact diagonalisation) +!> \param exp_eigvals Exponentiated eigenvalues (exact diagonalisation) +!> \param v_dbcsr Sparse matrix with bare Coulomb in RI basis +!> \param w_dbcsr Sparse matrix with correlation part of dressed Coulomb in RI basis (without bare Coulomb) +!> \param screened_dbt Tensor for screened Coulomb interaction +!> \param sigma_dbt Tensor for self-energy +!> \param greens_dbt Tensor for greens function/density matrix +!> \param t_3c_w Tensor containing 3c integrals +!> \param t_3c_work_RI__AO_AO Tensor sigma contraction +!> \param t_3c_work_RI_AO__AO Tensor sigma contraction +!> \param t_3c_work2_RI_AO__AO Tensor sigma contraction +!> \param sigma_SEX Screened exchange self-energy +!> \param sigma_COH Coulomb hole self-energy +!> \param hartree_curr Current Hartree matrix +!> \param etrs_max_iter Maximum number of ETRS iterations +!> \param ham_reference_type Which Hamiltonian to use as single particle basis +!> \param mat_exp_method Which method to use for matrix exponentiation +!> \param unit_nr Number of output unit +!> \param int_3c_array Array containing the local 3c integrals +!> \author Stepan Marek (01.24) +! ************************************************************************************************** + TYPE rtbse_env_type + INTEGER :: n_spin = 1, & + n_ao = -1, & + n_RI = -1 + INTEGER, DIMENSION(2) :: n_occ = -1 + REAL(kind=dp) :: spin_degeneracy = 2 + REAL(kind=dp), DIMENSION(3) :: field = 0.0_dp + TYPE(cp_fm_type), DIMENSION(:), POINTER :: moments => NULL(), & + moments_field => NULL() + INTEGER :: sim_step = 0, & + sim_start = 0, & + ! Needed for continuation runs for loading of previous moments trace + sim_start_orig = 0, & + sim_nsteps = -1 + REAL(kind=dp) :: sim_time = 0.0_dp, & + sim_dt = 0.1_dp, & + etrs_threshold = 1.0e-7_dp, & + exp_accuracy = 1.0e-10_dp, & + ft_damping = 0.0_dp, & + ft_start = 0.0_dp + ! Which element of polarizability to print out + INTEGER, DIMENSION(2) :: pol_element = [1, 1] + TYPE(dft_control_type), POINTER :: dft_control => NULL() + ! DEBUG : Trying keeping the reference to previous environments inside this one + TYPE(qs_environment_type), POINTER :: qs_env => NULL() + TYPE(post_scf_bandstructure_type), POINTER :: bs_env => NULL() + ! Stores data needed for reading/writing to the restart files + TYPE(section_vals_type), POINTER :: restart_section => NULL(), & + field_section => NULL(), & + rho_section => NULL(), & + ft_section => NULL(), & + pol_section => NULL(), & + moments_section => NULL() + LOGICAL :: restart_extracted = .FALSE. + + ! Different indices signify different spins + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_effective => NULL() + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_reference => NULL() + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_workspace => NULL() + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: sigma_SEX => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: sigma_COH => NULL(), & + hartree_curr => NULL() + + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho => NULL(), & + rho_new => NULL(), & + rho_new_last => NULL(), & + rho_M => NULL(), & + rho_orig => NULL() + TYPE(cp_fm_type) :: S_inv_fm = cp_fm_type(), & + S_fm = cp_fm_type() + ! Many routines require overlap in the complex format + TYPE(cp_cfm_type) :: S_cfm = cp_cfm_type() + TYPE(dbcsr_type) :: rho_dbcsr = dbcsr_type() + ! Indices only correspond to different workspaces + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_workspace => NULL() + ! Many methods use real and imaginary parts separately - prevent unnecessary reallocation + TYPE(cp_fm_type), DIMENSION(:), POINTER :: real_workspace => NULL() + ! Workspace required for exact matrix exponentiation + REAL(kind=dp), DIMENSION(:), POINTER :: real_eigvals => NULL() + COMPLEX(kind=dp), DIMENSION(:), POINTER :: exp_eigvals => NULL() + ! Workspace for saving the values for FT + TYPE(series_complex_type), DIMENSION(3) :: moments_trace = series_complex_type() + TYPE(series_real_type) :: time_trace = series_real_type() + TYPE(series_real_type), DIMENSION(3) :: field_trace = series_real_type() + ! Workspace required for hartree_pw + TYPE(dbcsr_type) :: v_dbcsr = dbcsr_type(), & + w_dbcsr = dbcsr_type() + TYPE(dbt_type) :: screened_dbt = dbt_type(), & + sigma_dbt = dbt_type(), & + greens_dbt = dbt_type(), & + t_3c_w = dbt_type(), & + t_3c_work_RI__AO_AO = dbt_type(), & + t_3c_work_RI_AO__AO = dbt_type(), & + t_3c_work2_RI_AO__AO = dbt_type() + ! These matrices are always real + INTEGER :: etrs_max_iter = 10 + INTEGER :: ham_reference_type = 2 + INTEGER :: mat_exp_method = 4 + INTEGER :: unit_nr = -1 + REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c_array => NULL() + + END TYPE rtbse_env_type + +CONTAINS + +! ************************************************************************************************** +!> \brief Allocates structures and prepares rtbse_env for run +!> \param rtbse_env rtbse_env_type that is initialised +!> \param qs_env Entry point of the calculation +!> \author Stepan Marek +!> \date 02.2024 +! ************************************************************************************************** + SUBROUTINE create_rtbse_env(rtbse_env, qs_env, force_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(force_env_type), POINTER :: force_env + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(rt_prop_type), POINTER :: rtp + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + INTEGER :: i, k + TYPE(section_vals_type), POINTER :: input, bs_sec, md_sec + INTEGER, DIMENSION(:), POINTER :: pol_tmp + + ! Allocate the storage for the gwbse environment + NULLIFY (rtbse_env) + ALLOCATE (rtbse_env) + ! Extract the other types first + CALL get_qs_env(qs_env, & + bs_env=bs_env, & + rtp=rtp, & + matrix_s=matrix_s, & + mos=mos, & + dft_control=rtbse_env%dft_control, & + input=input) + bs_sec => section_vals_get_subs_vals(input, "PROPERTIES%BANDSTRUCTURE") + IF (.NOT. ASSOCIATED(bs_env)) THEN + CPABORT("Cannot run RT-aGW without running GW calculation (PROPERTIES) before") + END IF + ! Number of spins + rtbse_env%n_spin = bs_env%n_spin + ! Number of atomic orbitals + rtbse_env%n_ao = bs_env%n_ao + ! Number of auxiliary basis orbitals + rtbse_env%n_RI = bs_env%n_RI + ! Number of occupied orbitals - for closed shell equals to half the number of electrons + rtbse_env%n_occ(:) = bs_env%n_occ(:) + ! Spin degeneracy - number of spins per orbital + rtbse_env%spin_degeneracy = bs_env%spin_degeneracy + ! Default field is zero + rtbse_env%field(:) = 0.0_dp + ! Default time is zero + rtbse_env%sim_step = 0 + rtbse_env%sim_time = 0 + ! Time step is taken from rtp + md_sec => section_vals_get_subs_vals(force_env%root_section, "MOTION%MD") + CALL section_vals_val_get(md_sec, "TIMESTEP", r_val=rtbse_env%sim_dt) + ! rtbse_env%sim_dt = rtp%dt + ! Threshold for etrs is taken from the eps_energy from RT propagation + rtbse_env%etrs_threshold = rtbse_env%dft_control%rtp_control%eps_ener + rtbse_env%exp_accuracy = rtbse_env%dft_control%rtp_control%eps_exp + ! Recover custom options + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE_HAMILTONIAN", & + i_val=rtbse_env%ham_reference_type) + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAX_ITER", & + i_val=rtbse_env%etrs_max_iter) + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAT_EXP", & + i_val=rtbse_env%mat_exp_method) + ! Output unit number, recovered from the post_scf_bandstructure_type + rtbse_env%unit_nr = bs_env%unit_nr + ! Sim start index and total number of steps as well + CALL section_vals_val_get(md_sec, "STEP_START_VAL", i_val=rtbse_env%sim_start) + ! Copy this value to sim_start_orig for continuation runs + rtbse_env%sim_start_orig = rtbse_env%sim_start + CALL section_vals_val_get(md_sec, "STEPS", i_val=rtbse_env%sim_nsteps) + ! Get the values for the FT + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS_FT%DAMPING", & + r_val=rtbse_env%ft_damping) + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS_FT%START_TIME", & + r_val=rtbse_env%ft_start) + ! TODO : Units for starting time and damping - so far give atomic units + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%PRINT%POLARIZABILITY%ELEMENT", & + i_vals=pol_tmp) + IF (SIZE(pol_tmp) < 2) CPABORT("Less than two elements provided for polarizability") + rtbse_env%pol_element(:) = pol_tmp(:) + ! Do basic sanity checks for pol_element + DO k = 1, 2 + IF (rtbse_env%pol_element(k) > 3 .OR. rtbse_env%pol_element(k) < 1) & + CPABORT("Polarisation tensor element not 1,2 or 3 in at least one index") + END DO + ! Get the restart section + rtbse_env%restart_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%RESTART") + rtbse_env%restart_extracted = .FALSE. + rtbse_env%field_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%FIELD") + rtbse_env%moments_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS") + rtbse_env%rho_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%DENSITY_MATRIX") + rtbse_env%ft_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%MOMENTS_FT") + rtbse_env%pol_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%POLARIZABILITY") + ! DEBUG : References to previous environments + rtbse_env%qs_env => qs_env + rtbse_env%bs_env => bs_env + + ! Allocate moments matrices + NULLIFY (rtbse_env%moments) + ALLOCATE (rtbse_env%moments(3)) + NULLIFY (rtbse_env%moments_field) + ALLOCATE (rtbse_env%moments_field(3)) + DO k = 1, 3 + ! Matrices are created from overlap template + ! Values are initialized in initialize_rtbse_env + CALL cp_fm_create(rtbse_env%moments(k), bs_env%fm_s_Gamma%matrix_struct) + CALL cp_fm_create(rtbse_env%moments_field(k), bs_env%fm_s_Gamma%matrix_struct) + END DO + + ! Allocate space for density propagation and other operations + NULLIFY (rtbse_env%rho_workspace) + ALLOCATE (rtbse_env%rho_workspace(4)) + DO i = 1, SIZE(rtbse_env%rho_workspace) + CALL cp_cfm_create(rtbse_env%rho_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_set_all(rtbse_env%rho_workspace(i), CMPLX(0.0, 0.0, kind=dp)) + END DO + ! Allocate real workspace + NULLIFY (rtbse_env%real_workspace) + SELECT CASE (rtbse_env%mat_exp_method) + CASE (do_exact) + ALLOCATE (rtbse_env%real_workspace(4)) + CASE (do_bch) + ALLOCATE (rtbse_env%real_workspace(2)) + CASE DEFAULT + CPABORT("Only exact and BCH matrix propagation implemented in RT-BSE") + END SELECT + DO i = 1, SIZE(rtbse_env%real_workspace) + CALL cp_fm_create(rtbse_env%real_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_fm_set_all(rtbse_env%real_workspace(i), 0.0_dp) + END DO + ! Allocate density matrix + NULLIFY (rtbse_env%rho) + ALLOCATE (rtbse_env%rho(rtbse_env%n_spin)) + DO i = 1, rtbse_env%n_spin + CALL cp_cfm_create(rtbse_env%rho(i), matrix_struct=bs_env%fm_s_Gamma%matrix_struct) + END DO + ! Create the inverse overlap matrix, for use in density propagation + ! Start by creating the actual overlap matrix + CALL cp_fm_create(rtbse_env%S_fm, bs_env%fm_s_Gamma%matrix_struct) + CALL cp_fm_create(rtbse_env%S_inv_fm, bs_env%fm_s_Gamma%matrix_struct) + CALL cp_cfm_create(rtbse_env%S_cfm, bs_env%fm_s_Gamma%matrix_struct) + + ! Create the single particle hamiltonian + ! Allocate workspace + NULLIFY (rtbse_env%ham_workspace) + ALLOCATE (rtbse_env%ham_workspace(rtbse_env%n_spin)) + DO i = 1, rtbse_env%n_spin + CALL cp_cfm_create(rtbse_env%ham_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_set_all(rtbse_env%ham_workspace(i), CMPLX(0.0, 0.0, kind=dp)) + END DO + ! Now onto the Hamiltonian itself + NULLIFY (rtbse_env%ham_reference) + ALLOCATE (rtbse_env%ham_reference(rtbse_env%n_spin)) + DO i = 1, rtbse_env%n_spin + CALL cp_cfm_create(rtbse_env%ham_reference(i), bs_env%fm_ks_Gamma(i)%matrix_struct) + END DO + + ! Create the matrices and workspaces for ETRS propagation + NULLIFY (rtbse_env%ham_effective) + NULLIFY (rtbse_env%rho_new) + NULLIFY (rtbse_env%rho_new_last) + NULLIFY (rtbse_env%rho_M) + NULLIFY (rtbse_env%rho_orig) + ALLOCATE (rtbse_env%ham_effective(rtbse_env%n_spin)) + ALLOCATE (rtbse_env%rho_new(rtbse_env%n_spin)) + ALLOCATE (rtbse_env%rho_new_last(rtbse_env%n_spin)) + ALLOCATE (rtbse_env%rho_M(rtbse_env%n_spin)) + ALLOCATE (rtbse_env%rho_orig(rtbse_env%n_spin)) + DO i = 1, rtbse_env%n_spin + CALL cp_cfm_create(rtbse_env%ham_effective(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_set_all(rtbse_env%ham_effective(i), CMPLX(0.0, 0.0, kind=dp)) + CALL cp_cfm_create(rtbse_env%rho_new(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_set_all(rtbse_env%rho_new(i), CMPLX(0.0, 0.0, kind=dp)) + CALL cp_cfm_create(rtbse_env%rho_new_last(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_set_all(rtbse_env%rho_new_last(i), CMPLX(0.0, 0.0, kind=dp)) + CALL cp_cfm_create(rtbse_env%rho_M(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_set_all(rtbse_env%rho_M(i), CMPLX(0.0, 0.0, kind=dp)) + CALL cp_cfm_create(rtbse_env%rho_orig(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + END DO + + ! Fields for exact diagonalisation + NULLIFY (rtbse_env%real_eigvals) + ALLOCATE (rtbse_env%real_eigvals(rtbse_env%n_ao)) + rtbse_env%real_eigvals(:) = 0.0_dp + NULLIFY (rtbse_env%exp_eigvals) + ALLOCATE (rtbse_env%exp_eigvals(rtbse_env%n_ao)) + rtbse_env%exp_eigvals(:) = CMPLX(0.0, 0.0, kind=dp) + + ! Workspace for FT - includes (in principle) the zeroth step and the extra last step + DO i = 1, 3 + NULLIFY (rtbse_env%moments_trace(i)%series) + ALLOCATE (rtbse_env%moments_trace(i)%series(rtbse_env%sim_nsteps + 2), source=z_zero) + NULLIFY (rtbse_env%field_trace(i)%series) + ALLOCATE (rtbse_env%field_trace(i)%series(rtbse_env%sim_nsteps + 2), source=0.0_dp) + END DO + NULLIFY (rtbse_env%time_trace%series) + ALLOCATE (rtbse_env%time_trace%series(rtbse_env%sim_nsteps + 2), source=0.0_dp) + + ! Allocate self-energy parts and dynamic Hartree potential + NULLIFY (rtbse_env%hartree_curr) + NULLIFY (rtbse_env%sigma_SEX) + NULLIFY (rtbse_env%sigma_COH) + ALLOCATE (rtbse_env%hartree_curr(rtbse_env%n_spin)) + ALLOCATE (rtbse_env%sigma_SEX(rtbse_env%n_spin)) + ALLOCATE (rtbse_env%sigma_COH(rtbse_env%n_spin)) + DO i = 1, rtbse_env%n_spin + CALL cp_fm_create(rtbse_env%sigma_COH(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_cfm_create(rtbse_env%sigma_SEX(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_fm_create(rtbse_env%hartree_curr(i), bs_env%fm_ks_Gamma(1)%matrix_struct) + CALL cp_fm_set_all(rtbse_env%sigma_COH(i), 0.0_dp) + CALL cp_cfm_set_all(rtbse_env%sigma_SEX(i), CMPLX(0.0, 0.0, kind=dp)) + CALL cp_fm_set_all(rtbse_env%hartree_curr(i), 0.0_dp) + END DO + + ! Allocate workspaces for get_sigma + CALL create_sigma_workspace(rtbse_env, qs_env) + + ! Depending on the chosen methods, allocate extra workspace + CALL create_hartree_ri_workspace(rtbse_env, qs_env) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Simple reimplementation of cp_fm_release_pp1 for complex matrices +!> \param matrices cp_cfm_p_type(:) +!> \author Stepan Marek +!> \date 02.2024 +! ************************************************************************************************** + SUBROUTINE cp_cfm_release_pa1(matrices) + TYPE(cp_cfm_type), DIMENSION(:), POINTER :: matrices + INTEGER :: i + + DO i = 1, SIZE(matrices) + CALL cp_cfm_release(matrices(i)) + END DO + DEALLOCATE (matrices) + NULLIFY (matrices) + END SUBROUTINE cp_cfm_release_pa1 + +! ************************************************************************************************** +!> \brief Releases the environment allocated structures +!> \param rtbse_env +!> \author Stepan Marek +!> \date 02.2024 +! ************************************************************************************************** + SUBROUTINE release_rtbse_env(rtbse_env) + TYPE(rtbse_env_type), POINTER :: rtbse_env + INTEGER :: i + + CALL cp_cfm_release_pa1(rtbse_env%ham_effective) + CALL cp_cfm_release_pa1(rtbse_env%ham_workspace) + CALL cp_fm_release(rtbse_env%sigma_COH) + CALL cp_cfm_release_pa1(rtbse_env%sigma_SEX) + CALL cp_fm_release(rtbse_env%hartree_curr) + CALL cp_cfm_release_pa1(rtbse_env%ham_reference) + CALL cp_cfm_release_pa1(rtbse_env%rho) + CALL cp_cfm_release_pa1(rtbse_env%rho_workspace) + CALL cp_cfm_release_pa1(rtbse_env%rho_new) + CALL cp_cfm_release_pa1(rtbse_env%rho_new_last) + CALL cp_cfm_release_pa1(rtbse_env%rho_M) + CALL cp_cfm_release_pa1(rtbse_env%rho_orig) + CALL cp_fm_release(rtbse_env%real_workspace) + CALL cp_fm_release(rtbse_env%S_inv_fm) + CALL cp_fm_release(rtbse_env%S_fm) + CALL cp_cfm_release(rtbse_env%S_cfm) + + ! DO i = 1, 3 + ! CALL cp_fm_release(rtbse_env%moments(i)%matrix) + ! CALL cp_fm_release(rtbse_env%moments_field(i)%matrix) + ! END DO + CALL cp_fm_release(rtbse_env%moments) + CALL cp_fm_release(rtbse_env%moments_field) + + CALL release_sigma_workspace(rtbse_env) + + CALL release_hartree_ri_workspace(rtbse_env) + + DEALLOCATE (rtbse_env%real_eigvals) + DEALLOCATE (rtbse_env%exp_eigvals) + DO i = 1, 3 + DEALLOCATE (rtbse_env%moments_trace(i)%series) + DEALLOCATE (rtbse_env%field_trace(i)%series) + END DO + DEALLOCATE (rtbse_env%time_trace%series) + + ! Deallocate the neighbour list that is not deallocated in gw anymore + IF (ASSOCIATED(rtbse_env%bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(rtbse_env%bs_env%nl_3c) + ! Deallocate the storage for the environment itself + DEALLOCATE (rtbse_env) + ! Nullify to make sure it is not used again + NULLIFY (rtbse_env) + + END SUBROUTINE +! ************************************************************************************************** +!> \brief Allocates the workspaces for Hartree RI method +!> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors +!> \param rtbse_env +!> \param qs_env Quickstep environment - entry point of calculation +!> \author Stepan Marek +!> \date 05.2024 +! ************************************************************************************************** + SUBROUTINE create_hartree_ri_workspace(rtbse_env, qs_env) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + REAL(kind=dp) :: size_mb + + CALL get_qs_env(qs_env, bs_env=bs_env) + + CALL dbcsr_create(rtbse_env%rho_dbcsr, name="Sparse density", template=bs_env%mat_ao_ao%matrix) + + ! v_dbcsr is created in init_hartree + + ! TODO : Implement option/decision to not precompute all the 3c integrals + size_mb = REAL(rtbse_env%n_ao*rtbse_env%n_ao*rtbse_env%n_RI*STORAGE_SIZE(size_mb))/(1024_dp*1024_dp) + IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, *) " RTBSE| Approximate size of int_3c in MB", size_mb + + ALLOCATE (rtbse_env%int_3c_array(rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%n_RI)) + CALL build_3c_integral_block(rtbse_env%int_3c_array, qs_env, potential_parameter=bs_env%ri_metric, & + basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, & + basis_i=bs_env%basis_set_RI, & + j_bf_start_from_atom=bs_env%i_ao_start_from_atom, & + k_bf_start_from_atom=bs_env%i_ao_start_from_atom, & + i_bf_start_from_atom=bs_env%i_RI_start_from_atom) + END SUBROUTINE create_hartree_ri_workspace +! ************************************************************************************************** +!> \brief Releases the workspace for the Hartree RI method +!> \param rtbse_env RT-BSE Environment, containing specific RI Hartree storage +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + SUBROUTINE release_hartree_ri_workspace(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + + DEALLOCATE (rtbse_env%int_3c_array) + + CALL dbcsr_release(rtbse_env%rho_dbcsr) + + CALL dbcsr_release(rtbse_env%v_dbcsr) + + END SUBROUTINE release_hartree_ri_workspace +! ************************************************************************************************** +!> \brief Allocates the workspaces for self-energy determination routine +!> \param rtbse_env Structure for holding information and workspace structures +!> \param qs_env Quickstep environment - entry point of calculation +!> \author Stepan Marek +!> \date 02.2024 +! ************************************************************************************************** + SUBROUTINE create_sigma_workspace(rtbse_env, qs_env) + TYPE(rtbse_env_type) :: rtbse_env + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(post_scf_bandstructure_type), POINTER :: bs_env + + CALL get_qs_env(qs_env, bs_env=bs_env) + + ! t_3c_w + CALL dbt_create(bs_env%t_RI__AO_AO, rtbse_env%t_3c_w) + ! TODO : Provide option/decision whether to store the 3c integrals precomputed + CALL compute_3c_integrals(qs_env, bs_env, rtbse_env%t_3c_w) + ! t_3c_work_RI_AO__AO + CALL dbt_create(bs_env%t_RI_AO__AO, rtbse_env%t_3c_work_RI_AO__AO) + ! t_3c_work2_RI_AO__AO + CALL dbt_create(bs_env%t_RI_AO__AO, rtbse_env%t_3c_work2_RI_AO__AO) + ! t_W + ! Populate screened_dbt from gw run + CALL dbcsr_create(rtbse_env%w_dbcsr, name="W", template=bs_env%mat_RI_RI%matrix) + CALL dbt_create(rtbse_env%w_dbcsr, rtbse_env%screened_dbt) + ! sigma_dbt + CALL dbt_create(bs_env%mat_ao_ao%matrix, rtbse_env%sigma_dbt) + ! greens_dbt + CALL dbt_create(bs_env%mat_ao_ao%matrix, rtbse_env%greens_dbt) + END SUBROUTINE +! ************************************************************************************************** +!> \brief Releases the workspaces for self-energy determination +!> \param rtbse_env +!> \author Stepan Marek +!> \date 02.2024 +! ************************************************************************************************** + SUBROUTINE release_sigma_workspace(rtbse_env) + TYPE(rtbse_env_type) :: rtbse_env + + CALL dbt_destroy(rtbse_env%t_3c_w) + CALL dbt_destroy(rtbse_env%t_3c_work_RI_AO__AO) + CALL dbt_destroy(rtbse_env%t_3c_work2_RI_AO__AO) + CALL dbt_destroy(rtbse_env%screened_dbt) + CALL dbt_destroy(rtbse_env%sigma_dbt) + CALL dbt_destroy(rtbse_env%greens_dbt) + CALL dbcsr_release(rtbse_env%w_dbcsr) + END SUBROUTINE +! ************************************************************************************************** +!> \brief Multiplies real matrix by a complex matrix from the right +!> \note So far only converts the real matrix to complex one, potentially doubling the work +!> \param rtbse_env +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + SUBROUTINE multiply_fm_cfm(trans_r, trans_c, na, nb, nc, & + alpha, matrix_r, matrix_c, beta, res) + ! Transposition + CHARACTER(len=1) :: trans_r, trans_c + INTEGER :: na, nb, nc + ! accept real numbers + ! TODO : Just use complex numbers and import z_one, z_zero etc. + REAL(kind=dp) :: alpha, beta + TYPE(cp_fm_type) :: matrix_r + TYPE(cp_cfm_type) :: matrix_c, res + TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im + REAL(kind=dp) :: i_unit + CHARACTER(len=1) :: trans_cr + + CALL cp_fm_create(work_re, matrix_c%matrix_struct) + CALL cp_fm_create(work_im, matrix_c%matrix_struct) + CALL cp_fm_create(res_re, res%matrix_struct) + CALL cp_fm_create(res_im, res%matrix_struct) + CALL cp_cfm_to_fm(matrix_c, work_re, work_im) + SELECT CASE (trans_c) + CASE ("C") + i_unit = -1.0_dp + trans_cr = "T" + CASE ("T") + i_unit = 1.0_dp + trans_cr = "T" + CASE default + i_unit = 1.0_dp + trans_cr = "N" + END SELECT + ! Actual multiplication + CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, & + alpha, matrix_r, work_re, beta, res_re) + CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, & + i_unit*alpha, matrix_r, work_im, beta, res_im) + CALL cp_fm_to_cfm(res_re, res_im, res) + CALL cp_fm_release(work_re) + CALL cp_fm_release(work_im) + CALL cp_fm_release(res_re) + CALL cp_fm_release(res_im) + + END SUBROUTINE multiply_fm_cfm +! ************************************************************************************************** +!> \brief Multiplies complex matrix by a real matrix from the right +!> \note So far only converts the real matrix to complex one, potentially doubling the work +!> \param rtbse_env +!> \author Stepan Marek +!> \date 09.2024 +! ************************************************************************************************** + SUBROUTINE multiply_cfm_fm(trans_c, trans_r, na, nb, nc, & + alpha, matrix_c, matrix_r, beta, res) + ! Transposition + CHARACTER(len=1) :: trans_c, trans_r + INTEGER :: na, nb, nc + ! accept real numbers + ! TODO : complex number support via interface? + REAL(kind=dp) :: alpha, beta + TYPE(cp_cfm_type) :: matrix_c, res + TYPE(cp_fm_type) :: matrix_r + TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im + REAL(kind=dp) :: i_unit + CHARACTER(len=1) :: trans_cr + + CALL cp_fm_create(work_re, matrix_c%matrix_struct) + CALL cp_fm_create(work_im, matrix_c%matrix_struct) + CALL cp_fm_create(res_re, res%matrix_struct) + CALL cp_fm_create(res_im, res%matrix_struct) + CALL cp_cfm_to_fm(matrix_c, work_re, work_im) + SELECT CASE (trans_c) + CASE ("C") + i_unit = -1.0_dp + trans_cr = "T" + CASE ("T") + i_unit = 1.0_dp + trans_cr = "T" + CASE default + i_unit = 1.0_dp + trans_cr = "N" + END SELECT + ! Actual multiplication + CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, & + alpha, work_re, matrix_r, beta, res_re) + CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, & + i_unit*alpha, work_im, matrix_r, beta, res_im) + CALL cp_fm_to_cfm(res_re, res_im, res) + CALL cp_fm_release(work_re) + CALL cp_fm_release(work_im) + CALL cp_fm_release(res_re) + CALL cp_fm_release(res_im) + + END SUBROUTINE multiply_cfm_fm +END MODULE diff --git a/src/emd/rt_propagator_init.F b/src/emd/rt_propagator_init.F index 2231231ebe..4f11e1ccdf 100644 --- a/src/emd/rt_propagator_init.F +++ b/src/emd/rt_propagator_init.F @@ -38,6 +38,7 @@ MODULE rt_propagator_init do_cn,& do_em,& do_etrs,& + do_exact,& do_pade,& do_taylor USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz @@ -160,7 +161,7 @@ SUBROUTINE init_propagators(qs_env) DO imat = 1, SIZE(mos_new) CALL cp_fm_to_fm(mos_new(imat), mos_next(imat)) END DO - ELSEIF (rtp_control%mat_exp == do_bch) THEN + ELSEIF (rtp_control%mat_exp == do_bch .OR. rtp_control%mat_exp == do_exact) THEN ELSE IF (rtp%linear_scaling) THEN CALL compute_exponential_sparse(exp_H_new, propagator_matrix, rtp_control, rtp) diff --git a/src/gw_large_cell_gamma.F b/src/gw_large_cell_gamma.F index d042b4c90e..d230163059 100644 --- a/src/gw_large_cell_gamma.F +++ b/src/gw_large_cell_gamma.F @@ -57,6 +57,7 @@ MODULE gw_large_cell_gamma USE gw_utils, ONLY: analyt_conti_and_print,& de_init_bs_env,& time_to_freq + USE input_constants, ONLY: rtp_method_bse USE input_section_types, ONLY: section_vals_type USE kinds, ONLY: default_string_length,& dp,& @@ -89,7 +90,8 @@ MODULE gw_large_cell_gamma CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_large_cell_gamma' - PUBLIC :: gw_calc_large_cell_Gamma + PUBLIC :: gw_calc_large_cell_Gamma, & + compute_3c_integrals CONTAINS @@ -959,6 +961,19 @@ SUBROUTINE read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time) IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' ' + ! Marek : Reading of the W(w=0) potential for RTP + ! TODO : is the condition bs_env%all_W_exist sufficient for reading? + IF (bs_env%rtp_method == rtp_method_bse) THEN + CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct) + t1 = m_walltime() + CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0) + IF (bs_env%unit_nr > 0) THEN + WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') & + 'Read W^MIC(w=0) from file for frequency point ', 1, ' /', 1, & + ', Execution time', m_walltime() - t1, ' s' + END IF + END IF + CALL timestop(handle) END SUBROUTINE read_W_MIC_time @@ -1052,6 +1067,23 @@ SUBROUTINE compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time) CALL cp_cfm_release(cfm_V_sqrt_ikp) CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau) + ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0 + ! TODO : Should this be activeted here or somewhere else? Maybe directly in the main (gw) routine? + ! Create the matrix TODO: Add release statement in some well understood place - at destruction of bs_env? + IF (bs_env%rtp_method == rtp_method_bse) THEN + CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct) + ! Set to zero + CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp) + ! Sum over all times + DO i_t = 1, bs_env%num_time_freq_points + ! Add the relevant structure with correct weight + CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, & + bs_env%imag_time_weights_freq_zero(i_t), fm_W_MIC_time(i_t)) + END DO + ! Done, save to file + CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env) + END IF + IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' ' CALL timestop(handle) diff --git a/src/gw_utils.F b/src/gw_utils.F index 4e717f71f1..4ce47c34e7 100644 --- a/src/gw_utils.F +++ b/src/gw_utils.F @@ -28,14 +28,9 @@ MODULE gw_utils cp_cfm_to_fm,& cp_cfm_type USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_api, ONLY: dbcsr_create,& - dbcsr_distribution_release,& - dbcsr_distribution_type,& - dbcsr_p_type,& - dbcsr_release,& - dbcsr_set,& - dbcsr_type,& - dbcsr_type_no_symmetry + USE cp_dbcsr_api, ONLY: & + dbcsr_create, dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, & + dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& copy_fm_to_dbcsr,& cp_dbcsr_dist2d_to_dist,& @@ -66,9 +61,11 @@ MODULE gw_utils USE input_constants, ONLY: do_potential_truncated,& large_cell_Gamma,& ri_rpa_g0w0_crossing_newton,& + rtp_method_bse,& small_cell_full_kp,& xc_none - USE input_section_types, ONLY: section_vals_get_subs_vals,& + USE input_section_types, ONLY: section_vals_get,& + section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get,& section_vals_val_set @@ -222,7 +219,8 @@ SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec) ! Recommendation in case of memory issues: first perform GW calculation without calculating ! LDOS (to safe memor). Then, use GW restart files ! in a subsequent calculation to calculate the LDOS - IF (.NOT. bs_env%do_ldos) THEN + ! Marek : TODO - boolean that does not interfere with RTP init but sets this to correct value + IF (.NOT. bs_env%do_ldos .AND. .FALSE.) THEN CALL qs_env_part_release(qs_env) END IF @@ -247,7 +245,11 @@ SUBROUTINE de_init_bs_env(bs_env) ! 2. consume a lot of memory and should not be kept until the quantity is ! deallocated in bs_env_release - IF (ASSOCIATED(bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(bs_env%nl_3c) + IF (ASSOCIATED(bs_env%nl_3c%ij_list) .AND. (bs_env%rtp_method == rtp_method_bse)) THEN + IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, *) "Retaining nl_3c for RTBSE" + ELSE + CALL neighbor_list_3c_destroy(bs_env%nl_3c) + END IF CALL cp_libint_static_cleanup() @@ -1594,12 +1596,14 @@ SUBROUTINE compute_V_xc(qs_env, bs_env) CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_xc' INTEGER :: handle, img, ispin, myfun, nimages - REAL(KIND=dp) :: energy_ex, energy_exc, energy_total + LOGICAL :: hf_present + REAL(KIND=dp) :: energy_ex, energy_exc, energy_total, & + myfraction TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_ks_without_v_xc TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp TYPE(dft_control_type), POINTER :: dft_control TYPE(qs_energy_type), POINTER :: energy - TYPE(section_vals_type), POINTER :: input, xc_section + TYPE(section_vals_type), POINTER :: hf_section, input, xc_section CALL timeset(routineN, handle) @@ -1613,6 +1617,17 @@ SUBROUTINE compute_V_xc(qs_env, bs_env) xc_section => section_vals_get_subs_vals(input, "DFT%XC") CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun) CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=xc_none) + ! IF (ASSOCIATED(section_vals_get_subs_vals(xc_section, "HF", can_return_null=.TRUE.))) THEN + hf_section => section_vals_get_subs_vals(input, "DFT%XC%HF", can_return_null=.TRUE.) + hf_present = .FALSE. + IF (ASSOCIATED(hf_section)) THEN + CALL section_vals_get(hf_section, explicit=hf_present) + END IF + IF (hf_present) THEN + ! Special case for handling hfx + CALL section_vals_val_get(xc_section, "HF%FRACTION", r_val=myfraction) + CALL section_vals_val_set(xc_section, "HF%FRACTION", r_val=0.0_dp) + END IF ! save the energy before the energy gets updated energy_total = energy%total @@ -1627,7 +1642,12 @@ SUBROUTINE compute_V_xc(qs_env, bs_env) DO ispin = 1, bs_env%n_spin ALLOCATE (mat_ks_without_v_xc(ispin)%matrix) - CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix) + IF (hf_present) THEN + CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix, & + matrix_type=dbcsr_type_symmetric) + ELSE + CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix) + END IF END DO ! calculate KS-matrix without XC @@ -1675,7 +1695,12 @@ SUBROUTINE compute_V_xc(qs_env, bs_env) dft_control%nimages = nimages ! set the DFT functional and HF fraction back - CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun) + CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", & + i_val=myfun) + IF (hf_present) THEN + CALL section_vals_val_set(xc_section, "HF%FRACTION", & + r_val=myfraction) + END IF IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN ! calculate KS-matrix again with XC @@ -1717,6 +1742,7 @@ SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env) ALLOCATE (bs_env%imag_freq_points(num_time_freq_points)) ALLOCATE (bs_env%imag_time_points(num_time_freq_points)) + ALLOCATE (bs_env%imag_time_weights_freq_zero(num_time_freq_points)) ALLOCATE (bs_env%weights_cos_t_to_w(num_time_freq_points, num_time_freq_points)) ALLOCATE (bs_env%weights_cos_w_to_t(num_time_freq_points, num_time_freq_points)) ALLOCATE (bs_env%weights_sin_t_to_w(num_time_freq_points, num_time_freq_points)) @@ -1789,6 +1815,7 @@ SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env) END IF bs_env%imag_time_points(:) = points_and_weights(1:num_time_freq_points)/(2.0_dp*E_min) + bs_env%imag_time_weights_freq_zero(:) = points_and_weights(num_time_freq_points + 1:)/(E_min) DEALLOCATE (points_and_weights) diff --git a/src/input_constants.F b/src/input_constants.F index 300b88093a..b9b0f566d4 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -920,7 +920,8 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: do_taylor = 1, & do_pade = 2, & do_arnoldi = 3, & - do_bch = 4 + do_bch = 4, & + do_exact = 5 INTEGER, PARAMETER, PUBLIC :: do_etrs = 1, & do_cn = 2, & @@ -935,6 +936,12 @@ MODULE input_constants ramp_env = 3, & custom_env = 4 + INTEGER, PARAMETER, PUBLIC :: rtp_method_tddft = 1, & + rtp_method_bse = 2 + + INTEGER, PARAMETER, PUBLIC :: rtp_bse_ham_ks = 1, & + rtp_bse_ham_g0w0 = 2 + ! type of reference for RTP%PRINT%PROJECTION_MO INTEGER, PARAMETER, PUBLIC :: proj_mo_ref_scf = 1, & proj_mo_ref_xas_tdp = 2 diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 1ce2aa7c8f..e9aa3a3ec8 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -39,18 +39,18 @@ MODULE input_cp2k_dft do_admm_exch_scaling_none, do_admm_purify_cauchy, do_admm_purify_cauchy_subspace, & do_admm_purify_mcweeny, do_admm_purify_mo_diag, do_admm_purify_mo_no_diag, & do_admm_purify_none, do_admm_purify_none_dm, do_arnoldi, do_bch, do_cn, do_em, do_etrs, & - do_pade, do_taylor, ehrenfest, gaussian, kg_color_dsatur, kg_color_greedy, & + do_exact, do_pade, do_taylor, ehrenfest, gaussian, kg_color_dsatur, kg_color_greedy, & kg_tnadd_atomic, kg_tnadd_embed, kg_tnadd_embed_ri, kg_tnadd_none, no_admm_type, & no_excitations, numerical, oe_gllb, oe_lb, oe_none, oe_saop, oe_sic, plus_u_lowdin, & plus_u_mulliken, plus_u_mulliken_charges, real_time_propagation, rel_dkh, rel_none, & rel_pot_erfc, rel_pot_full, rel_sczora_mp, rel_trans_atom, rel_trans_full, & - rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, sccs_andreussi, & - sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, sccs_derivative_fft, & - sccs_fattebert_gygi, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, & - sic_mauri_us, sic_none, slater, tddfpt_davidson, tddfpt_excitations, tddfpt_lanczos, & - tddfpt_singlet, tddfpt_triplet, use_mom_ref_coac, use_mom_ref_com, use_mom_ref_user, & - use_mom_ref_zero, use_restart_wfn, use_rt_restart, use_scf_wfn, weight_type_mass, & - weight_type_unit + rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, rtp_bse_ham_g0w0, & + rtp_bse_ham_ks, rtp_method_bse, rtp_method_tddft, sccs_andreussi, sccs_derivative_cd3, & + sccs_derivative_cd5, sccs_derivative_cd7, sccs_derivative_fft, sccs_fattebert_gygi, & + sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, & + slater, tddfpt_davidson, tddfpt_excitations, tddfpt_lanczos, tddfpt_singlet, & + tddfpt_triplet, use_mom_ref_coac, use_mom_ref_com, use_mom_ref_user, use_mom_ref_zero, & + use_restart_wfn, use_rt_restart, use_scf_wfn, weight_type_mass, weight_type_unit USE input_cp2k_almo, ONLY: create_almo_scf_section USE input_cp2k_as, ONLY: create_active_space_section USE input_cp2k_ec, ONLY: create_ec_section @@ -1630,15 +1630,17 @@ SUBROUTINE create_rtp_section(section) " in the propagator. It is recommended to use BCH when employing density_propagation "// & "and ARNOLDI otherwise.", & usage="MAT_EXP TAYLOR", default_i_val=do_arnoldi, & - enum_c_vals=s2a("TAYLOR", "PADE", "ARNOLDI", "BCH"), & - enum_i_vals=(/do_taylor, do_pade, do_arnoldi, do_bch/), & + enum_c_vals=s2a("TAYLOR", "PADE", "ARNOLDI", "BCH", "EXACT"), & + enum_i_vals=(/do_taylor, do_pade, do_arnoldi, do_bch, do_exact/), & enum_desc=s2a("exponential is evaluated using scaling and squaring in combination"// & " with a taylor expansion of the exponential.", & "uses scaling and squaring together with the pade approximation", & "uses arnoldi subspace algorithm to compute exp(H)*MO directly, can't be used in "// & "combination with Crank Nicholson or density propagation", & "Uses a Baker-Campbell-Hausdorff expansion to propagate the density matrix,"// & - " only works for density propagation")) + " only works for density propagation", & + "Uses diagonalisation of the exponent matrices to determine the "// & + "matrix exponential exactly. Only implemented for GWBSE.")) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -1668,6 +1670,33 @@ SUBROUTINE create_rtp_section(section) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + ! Marek : Controlling flow to GWBSE + CALL keyword_create(keyword, __LOCATION__, name="RTP_METHOD", & + description="Which method is used for the time propagation of electronic structure. "// & + "By default, use the TDDFT method. Can also choose RT-BSE method, which propagates the lesser "// & + "Green's function instead of density matrix/molecular orbitals.", & + usage="RTP_METHOD RTBSE", & + default_i_val=rtp_method_tddft, & + enum_c_vals=s2a("TDDFT", "RTBSE"), & + enum_i_vals=(/rtp_method_tddft, rtp_method_bse/), & + enum_desc=s2a("Use TDDFT for density matrix/MO propagation.", & + "Use RT-BSE for Green's function propagation")) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + ! Marek : Development option - run GWBSE starting from the KS Hamiltonian + CALL keyword_create(keyword, __LOCATION__, name="RTBSE_HAMILTONIAN", & + description="Which Hamiltonian to use as the single-particle Hamiltonian"// & + " in the Green's propagator.", & + usage="RTBSE_HAMILTONIAN G0W0", & + default_i_val=rtp_bse_ham_g0w0, & + enum_c_vals=s2a("KS", "G0W0"), & + enum_i_vals=(/rtp_bse_ham_ks, rtp_bse_ham_g0w0/), & + enum_desc=s2a("Use Kohn-Sham Hamiltonian for Green's propagation.", & + "Use G0W0 Hamiltonian for Green's function propagation")) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="PROPAGATOR", & description="Which propagator should be used for the orbitals", & usage="PROPAGATOR ETRS", default_i_val=do_etrs, & @@ -1860,7 +1889,8 @@ SUBROUTINE create_rtp_section(section) CALL cp_print_key_section_create(print_key, __LOCATION__, "RESTART", & description="Controls the dumping of the MO restart file during rtp. "// & "By default keeps a short history of three restarts. "// & - "See also RESTART_HISTORY. In density propagation this controls the printing of P.", & + "See also RESTART_HISTORY. In density propagation this controls the printing of "// & + "density matrix.", & print_level=low_print_level, common_iter_levels=3, & each_iter_names=s2a("MD"), each_iter_values=(/20/), & add_last=add_last_numeric, filename="RESTART") @@ -1926,6 +1956,72 @@ SUBROUTINE create_rtp_section(section) CALL section_add_subsection(print_section, print_key) CALL section_release(print_key) + ! Marek : Add print option for ASCII density files - DEVELPMENT ONLY? + CALL cp_print_key_section_create(print_key, __LOCATION__, "DENSITY_MATRIX", & + description="Prints the density matrix at iterations in clear text to a file", & + print_level=high_print_level, common_iter_levels=0, & + each_iter_names=s2a("MD"), & + each_iter_values=(/1/), & + filename="rho") + CALL section_add_subsection(print_section, print_key) + CALL section_release(print_key) + ! Marek : Moments ASCII print + CALL cp_print_key_section_create(print_key, __LOCATION__, "MOMENTS", & + description="Prints the time-dependent electronic moments at "// & + "iterations in clear text to a file.", & + print_level=high_print_level, common_iter_levels=0, & + each_iter_names=s2a("MD"), & + each_iter_values=(/1/), & + filename="__STD_OUT__") + CALL section_add_subsection(print_section, print_key) + CALL section_release(print_key) + ! Marek : Fourier transform of MOMENTS ASCII print + CALL cp_print_key_section_create(print_key, __LOCATION__, "MOMENTS_FT", & + description="Prints the calculated Fourier transform of "// & + "time-dependent moments. For calculations with real time pulse (not delta kick) "// & + "can be supplied with starting time.", & + print_level=medium_print_level, common_iter_levels=0, & + each_iter_names=s2a("MD"), & + each_iter_values=(/1/), & + filename="MOMENTS_FT") + CALL keyword_create(keyword, __LOCATION__, "DAMPING", & + description="Sets the exponential damping of the time trace before carrying out Fourier transform. "// & + "For negative values (the default), calculates the damping at the run time so that the last point "// & + "in the time trace is reduced by factor e^(-4). Uses atomic units of inverse time.", & + type_of_var=real_t, & + default_r_val=-1.0_dp) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, "START_TIME", & + description="The starting time from which damping is applied and from which on the trace is "// & + "considered for the Fourier transform. Useful for real-time pulse - "// & + "one can specify the center of the pulse as the starting point.", & + type_of_var=real_t, & + default_r_val=0.0_dp) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL section_add_subsection(print_section, print_key) + CALL section_release(print_key) + ! Marek : Chosen element of (Fourier transformed) polarizability tensor (energy dependent) - text format + CALL cp_print_key_section_create(print_key, __LOCATION__, "POLARIZABILITY", & + description="Prints the chosen element of the energy dependent polarizability tensor "// & + "to a specified file. The tensor is calculated as ratio of "// & + "Fourier transform of the dipole "// & + "moment trace and Fourier transform of the applied field "// & + "(for delta kick, constant real field is applied.", & + print_level=medium_print_level, common_iter_levels=0, & + each_iter_names=s2a("MD"), & + each_iter_values=(/1/), & + filename="POLARIZABILITY") + CALL keyword_create(keyword, __LOCATION__, "ELEMENT", & + description="Specifies the element of polarizability which is to be printed out "// & + "(indexing starts at 1).", & + type_of_var=integer_t, default_i_vals=(/1, 1/), n_var=2, usage="ELEMENT 1 1") + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL section_add_subsection(print_section, print_key) + CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, "E_CONSTITUENTS", & description="Print the energy constituents (relevant to RTP) which make up "// & "the Total Energy", & diff --git a/src/post_scf_bandstructure_types.F b/src/post_scf_bandstructure_types.F index 46d12b0405..fc24486c6f 100644 --- a/src/post_scf_bandstructure_types.F +++ b/src/post_scf_bandstructure_types.F @@ -124,20 +124,21 @@ MODULE post_scf_bandstructure_types INTEGER, DIMENSION(3) :: periodic = -1 REAL(KIND=dp), DIMENSION(3, 3) :: hmat = -1.0_dp - ! imaginary time and imaginary frequency grids - INTEGER :: num_time_freq_points = -1, & - num_freq_points_fit = -1 - REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: imag_time_points, & - imag_freq_points, & - imag_freq_points_fit - REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: weights_cos_t_to_w, & - weights_cos_w_to_t, & - weights_sin_t_to_w - INTEGER :: nparam_pade = -1, & - num_points_per_magnitude = -1 - REAL(KIND=dp) :: freq_max_fit = -1.0_dp, & - regularization_minimax = -1.0_dp, & - stabilize_exp = -1.0_dp + ! imaginary time and frequency grids + INTEGER :: num_time_freq_points = -1, & + num_freq_points_fit = -1 + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: imag_time_points, & + imag_time_weights_freq_zero, & + imag_freq_points, & + imag_freq_points_fit + REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: weights_cos_t_to_w, & + weights_cos_w_to_t, & + weights_sin_t_to_w + INTEGER :: nparam_pade = -1, & + num_points_per_magnitude = -1 + REAL(KIND=dp) :: freq_max_fit = -1.0_dp, & + regularization_minimax = -1.0_dp, & + stabilize_exp = -1.0_dp ! filter threshold for matrix-tensor operations REAL(KIND=dp) :: eps_filter = -1.0_dp, & @@ -161,7 +162,9 @@ MODULE post_scf_bandstructure_types fm_chi_Gamma_freq, & fm_W_MIC_freq, & fm_W_MIC_freq_1_extra, & - fm_W_MIC_freq_1_no_extra + fm_W_MIC_freq_1_no_extra, & + fm_W_MIC_freq_zero, & + fm_h_G0W0_Gamma TYPE(cp_cfm_type) :: cfm_work_mo, & cfm_work_mo_2 @@ -209,6 +212,8 @@ MODULE post_scf_bandstructure_types atoms_j_t_group LOGICAL, DIMENSION(:, :), ALLOCATABLE :: skip_Sigma_occ, & skip_Sigma_vir + ! Marek : rtp_method + INTEGER :: rtp_method ! check-arrays and names for restarting LOGICAL, DIMENSION(:), ALLOCATABLE :: read_chi, & @@ -317,6 +322,7 @@ SUBROUTINE bs_env_release(bs_env) IF (ALLOCATED(bs_env%l_RI)) DEALLOCATE (bs_env%l_RI) IF (ALLOCATED(bs_env%xkp_special)) DEALLOCATE (bs_env%xkp_special) IF (ALLOCATED(bs_env%imag_time_points)) DEALLOCATE (bs_env%imag_time_points) + IF (ALLOCATED(bs_env%imag_time_weights_freq_zero)) DEALLOCATE (bs_env%imag_time_weights_freq_zero) IF (ALLOCATED(bs_env%imag_freq_points)) DEALLOCATE (bs_env%imag_freq_points) IF (ALLOCATED(bs_env%eigenval_scf_Gamma)) DEALLOCATE (bs_env%eigenval_scf_Gamma) IF (ALLOCATED(bs_env%eigenval_scf)) DEALLOCATE (bs_env%eigenval_scf) @@ -363,6 +369,8 @@ SUBROUTINE bs_env_release(bs_env) CALL cp_fm_release(bs_env%fm_RI_RI) CALL cp_fm_release(bs_env%fm_chi_Gamma_freq) CALL cp_fm_release(bs_env%fm_W_MIC_freq) + ! TODO : Test whether IF (rtp_method == rtp_method_gw) is needed + CALL cp_fm_release(bs_env%fm_W_MIC_freq_zero) CALL cp_fm_release(bs_env%fm_W_MIC_freq_1_extra) CALL cp_fm_release(bs_env%fm_W_MIC_freq_1_no_extra) CALL cp_cfm_release(bs_env%cfm_work_mo) diff --git a/src/post_scf_bandstructure_utils.F b/src/post_scf_bandstructure_utils.F index 362289ee3e..6bbb474c05 100644 --- a/src/post_scf_bandstructure_utils.F +++ b/src/post_scf_bandstructure_utils.F @@ -832,6 +832,7 @@ SUBROUTINE get_parameters_from_qs_env(qs_env, bs_env) TYPE(mp_para_env_type), POINTER :: para_env TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(scf_control_type), POINTER :: scf_control + TYPE(section_vals_type), POINTER :: input CALL timeset(routineN, handle) @@ -883,6 +884,10 @@ SUBROUTINE get_parameters_from_qs_env(qs_env, bs_env) u = bs_env%unit_nr + ! Marek : Get and save the rtp method + CALL get_qs_env(qs_env=qs_env, input=input) + CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTP_METHOD", i_val=bs_env%rtp_method) + IF (u > 0) THEN WRITE (u, FMT="(T2,2A,T73,I8)") "Number of occupied molecular orbitals (MOs) ", & "= Number of occupied bands", homo diff --git a/src/start/PACKAGE b/src/start/PACKAGE index cd22196069..d1907e7c10 100644 --- a/src/start/PACKAGE +++ b/src/start/PACKAGE @@ -16,5 +16,6 @@ "../swarm", "../motion", "../motion/mc", + "../emd", ], } diff --git a/src/start/cp2k_runs.F b/src/start/cp2k_runs.F index 8977fb2638..eabc653ca9 100644 --- a/src/start/cp2k_runs.F +++ b/src/start/cp2k_runs.F @@ -80,7 +80,7 @@ MODULE cp2k_runs do_sirius, do_swarm, do_tamc, do_test, do_tree_mc, do_tree_mc_ana, driver_run, ehrenfest, & electronic_spectra_run, energy_force_run, energy_run, geo_opt_run, linear_response_run, & mol_dyn_run, mon_car_run, negf_run, none_run, pint_run, real_time_propagation, & - tree_mc_run, vib_anal + rtp_method_bse, tree_mc_run, vib_anal USE input_cp2k, ONLY: create_cp2k_root_section USE input_cp2k_check, ONLY: check_cp2k_input USE input_cp2k_global, ONLY: create_global_section @@ -127,6 +127,7 @@ MODULE cp2k_runs USE qs_linres_module, ONLY: linres_calculation USE qs_tddfpt_module, ONLY: tddfpt_calculation USE reference_manager, ONLY: export_references_as_xml + USE rt_bse, ONLY: run_propagation_bse USE rt_propagation, ONLY: rt_prop_setup USE sirius_interface, ONLY: cp_sirius_finalize,& cp_sirius_init @@ -383,7 +384,14 @@ RECURSIVE SUBROUTINE cp2k_run(input_declaration, input_file_name, output_unit, m CPABORT("Real time propagation needs METHOD QS. ") CALL get_qs_env(force_env%qs_env, dft_control=dft_control) dft_control%rtp_control%fixed_ions = .TRUE. - CALL rt_prop_setup(force_env) + SELECT CASE (dft_control%rtp_control%rtp_method) + CASE (rtp_method_bse) + ! Run the TD-BSE method + CALL run_propagation_bse(force_env%qs_env, force_env) + CASE default + ! Run the TDDFT method + CALL rt_prop_setup(force_env) + END SELECT CASE (ehrenfest) IF (method_name_id /= do_qs) & CPABORT("Ehrenfest dynamics needs METHOD QS ") diff --git a/tests/QS/regtest-rtbse/TEST_FILES b/tests/QS/regtest-rtbse/TEST_FILES new file mode 100644 index 0000000000..215e425dc7 --- /dev/null +++ b/tests/QS/regtest-rtbse/TEST_FILES @@ -0,0 +1,6 @@ +h2_rt.inp 115 1e-3 0.14574977E-005 +h2_rt.inp 116 1e-5 0.00000000E+000 +h2_rt_cont.inp 117 1e-3 0.86507260E-006 +h2_delta.inp 115 1e-3 -0.21190408E-001 +h2_delta_cont.inp 117 1e-3 0.88156804E-002 +h2_ft.inp 118 1e-3 0.91127346E+001 diff --git a/tests/QS/regtest-rtbse/h2_delta.inp b/tests/QS/regtest-rtbse/h2_delta.inp new file mode 100644 index 0000000000..557a169e9c --- /dev/null +++ b/tests/QS/regtest-rtbse/h2_delta.inp @@ -0,0 +1,101 @@ +&GLOBAL + PROJECT_NAME H2_DELTA + RUN_TYPE RT_PROPAGATION +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 50 + TEMPERATURE [K] 0.0 + TIMESTEP [fs] 0.002 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + BASIS_SET_FILE_NAME BASIS_RI_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + ! Filter for density matrix + &LS_SCF + EPS_FILTER 1.0e-20 + &END LS_SCF + &MGRID + CUTOFF 100 + NGRIDS 1 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-7 + METHOD GAPW + &END QS + &REAL_TIME_PROPAGATION + APPLY_DELTA_PULSE + ASPC_ORDER 1 + DELTA_PULSE_DIRECTION 1 0 0 + DELTA_PULSE_SCALE 0.01 + DENSITY_PROPAGATION ON + EPS_ITER 1.0E-8 + EXP_ACCURACY 1.0E-10 + INITIAL_WFN SCF_WFN + MAT_EXP BCH + PERIODIC .FALSE. + RTP_METHOD RTBSE + &PRINT + &MOMENTS + &END MOMENTS + &END PRINT + &END REAL_TIME_PROPAGATION + &SCF + ADDED_MOS -1 + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &DIAGONALIZATION ON + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.4 + METHOD BROYDEN_MIXING + NBROYDEN 4 + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + ! Include G0W0 bandstructure correction + &PROPERTIES + &BANDSTRUCTURE + &GW + ! Simplest setting for quick checks + NUM_TIME_FREQ_POINTS 20 + &END GW + &END BANDSTRUCTURE + &END PROPERTIES + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + PERIODIC NONE + &END CELL + &COORD + H 0.0 0.0 0.0 + H 0.74 0.0 0.0 + &END COORD + &KIND H + BASIS_SET ORB SZV-MOLOPT-GTH + BASIS_SET RI_AUX RI_11_2_3_0_0_0_0_0_3.5e-03 + POTENTIAL GTH-PBE + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rtbse/h2_delta_cont.inp b/tests/QS/regtest-rtbse/h2_delta_cont.inp new file mode 100644 index 0000000000..0ee53ffbe6 --- /dev/null +++ b/tests/QS/regtest-rtbse/h2_delta_cont.inp @@ -0,0 +1,101 @@ +&GLOBAL + PROJECT_NAME H2_DELTA + RUN_TYPE RT_PROPAGATION +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 100 + TEMPERATURE [K] 0.0 + TIMESTEP [fs] 0.002 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + BASIS_SET_FILE_NAME BASIS_RI_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + ! Filter for density matrix + &LS_SCF + EPS_FILTER 1.0e-20 + &END LS_SCF + &MGRID + CUTOFF 100 + NGRIDS 1 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-7 + METHOD GAPW + &END QS + &REAL_TIME_PROPAGATION + APPLY_DELTA_PULSE + ASPC_ORDER 1 + DELTA_PULSE_DIRECTION 1 0 0 + DELTA_PULSE_SCALE 0.01 + DENSITY_PROPAGATION ON + EPS_ITER 1.0E-8 + EXP_ACCURACY 1.0E-10 + INITIAL_WFN RT_RESTART + MAT_EXP BCH + PERIODIC .FALSE. + RTP_METHOD RTBSE + &PRINT + &MOMENTS + &END MOMENTS + &END PRINT + &END REAL_TIME_PROPAGATION + &SCF + ADDED_MOS -1 + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &DIAGONALIZATION ON + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.4 + METHOD BROYDEN_MIXING + NBROYDEN 4 + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + ! Include G0W0 bandstructure correction + &PROPERTIES + &BANDSTRUCTURE + &GW + ! Simplest setting for quick checks + NUM_TIME_FREQ_POINTS 20 + &END GW + &END BANDSTRUCTURE + &END PROPERTIES + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + PERIODIC NONE + &END CELL + &COORD + H 0.0 0.0 0.0 + H 0.74 0.0 0.0 + &END COORD + &KIND H + BASIS_SET ORB SZV-MOLOPT-GTH + BASIS_SET RI_AUX RI_11_2_3_0_0_0_0_0_3.5e-03 + POTENTIAL GTH-PBE + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rtbse/h2_ft.inp b/tests/QS/regtest-rtbse/h2_ft.inp new file mode 100644 index 0000000000..a9839a8ca0 --- /dev/null +++ b/tests/QS/regtest-rtbse/h2_ft.inp @@ -0,0 +1,102 @@ +&GLOBAL + PROJECT_NAME H2_FT + RUN_TYPE RT_PROPAGATION +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 100 + TEMPERATURE [K] 0.0 + TIMESTEP [fs] 0.01 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + BASIS_SET_FILE_NAME BASIS_RI_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + ! Filter for density matrix + &LS_SCF + EPS_FILTER 1.0e-20 + &END LS_SCF + &MGRID + CUTOFF 100 + NGRIDS 1 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-7 + METHOD GAPW + &END QS + &REAL_TIME_PROPAGATION + APPLY_DELTA_PULSE + ASPC_ORDER 1 + DELTA_PULSE_DIRECTION 1 0 0 + DELTA_PULSE_SCALE 0.01 + DENSITY_PROPAGATION ON + EPS_ITER 1.0E-8 + EXP_ACCURACY 1.0E-10 + INITIAL_WFN SCF_WFN + MAT_EXP BCH + PERIODIC .FALSE. + RTP_METHOD RTBSE + &PRINT + &POLARIZABILITY + FILENAME __STD_OUT__ + &END POLARIZABILITY + &END PRINT + &END REAL_TIME_PROPAGATION + &SCF + ADDED_MOS -1 + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &DIAGONALIZATION ON + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.4 + METHOD BROYDEN_MIXING + NBROYDEN 4 + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + ! Include G0W0 bandstructure correction + &PROPERTIES + &BANDSTRUCTURE + &GW + ! Simplest setting for quick checks + NUM_TIME_FREQ_POINTS 20 + &END GW + &END BANDSTRUCTURE + &END PROPERTIES + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + PERIODIC NONE + &END CELL + &COORD + H 0.0 0.0 0.0 + H 0.74 0.0 0.0 + &END COORD + &KIND H + BASIS_SET ORB SZV-MOLOPT-GTH + BASIS_SET RI_AUX RI_11_2_3_0_0_0_0_0_3.5e-03 + POTENTIAL GTH-PBE + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rtbse/h2_rt.inp b/tests/QS/regtest-rtbse/h2_rt.inp new file mode 100644 index 0000000000..fc925425d5 --- /dev/null +++ b/tests/QS/regtest-rtbse/h2_rt.inp @@ -0,0 +1,109 @@ +&GLOBAL + PROJECT_NAME H2_RT + RUN_TYPE RT_PROPAGATION +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 50 + TEMPERATURE [K] 0.0 + TIMESTEP [fs] 0.002 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + BASIS_SET_FILE_NAME BASIS_RI_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + &EFIELD + ENVELOP GAUSSIAN + INTENSITY 1.0E+5 + PHASE 0.0 + POLARISATION 1 0 0 + WAVELENGTH 1000.0 + &GAUSSIAN_ENV + SIGMA [fs] 0.005 + T0 [fs] 0.05 + &END GAUSSIAN_ENV + &END EFIELD + ! Filter for density matrix + &LS_SCF + EPS_FILTER 1.0e-20 + &END LS_SCF + &MGRID + CUTOFF 100 + NGRIDS 1 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-7 + METHOD GAPW + &END QS + &REAL_TIME_PROPAGATION + ASPC_ORDER 1 + DENSITY_PROPAGATION ON + EPS_ITER 1.0E-8 + EXP_ACCURACY 1.0E-10 + INITIAL_WFN SCF_WFN + MAT_EXP BCH + PERIODIC .FALSE. + RTP_METHOD RTBSE + &PRINT + &MOMENTS + &END MOMENTS + &END PRINT + &END REAL_TIME_PROPAGATION + &SCF + ADDED_MOS -1 + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &DIAGONALIZATION ON + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.4 + METHOD BROYDEN_MIXING + NBROYDEN 4 + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + ! Include G0W0 bandstructure correction + &PROPERTIES + &BANDSTRUCTURE + &GW + ! Simplest setting for quick checks + NUM_TIME_FREQ_POINTS 20 + &END GW + &END BANDSTRUCTURE + &END PROPERTIES + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + PERIODIC NONE + &END CELL + &COORD + H 0.0 0.0 0.0 + H 0.74 0.0 0.0 + &END COORD + &KIND H + BASIS_SET ORB SZV-MOLOPT-GTH + BASIS_SET RI_AUX RI_11_2_3_0_0_0_0_0_3.5e-03 + POTENTIAL GTH-PBE + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rtbse/h2_rt_cont.inp b/tests/QS/regtest-rtbse/h2_rt_cont.inp new file mode 100644 index 0000000000..19ae5d38a8 --- /dev/null +++ b/tests/QS/regtest-rtbse/h2_rt_cont.inp @@ -0,0 +1,109 @@ +&GLOBAL + PROJECT_NAME H2_RT + RUN_TYPE RT_PROPAGATION +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 100 + TEMPERATURE [K] 0.0 + TIMESTEP [fs] 0.002 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD QS + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT + BASIS_SET_FILE_NAME BASIS_RI_MOLOPT + POTENTIAL_FILE_NAME GTH_POTENTIALS + &EFIELD + ENVELOP GAUSSIAN + INTENSITY 1.0E+5 + PHASE 0.0 + POLARISATION 1 0 0 + WAVELENGTH 1000.0 + &GAUSSIAN_ENV + SIGMA [fs] 0.005 + T0 [fs] 0.05 + &END GAUSSIAN_ENV + &END EFIELD + ! Filter for density matrix + &LS_SCF + EPS_FILTER 1.0e-20 + &END LS_SCF + &MGRID + CUTOFF 100 + NGRIDS 1 + REL_CUTOFF 10 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-7 + METHOD GAPW + &END QS + &REAL_TIME_PROPAGATION + ASPC_ORDER 1 + DENSITY_PROPAGATION ON + EPS_ITER 1.0E-8 + EXP_ACCURACY 1.0E-10 + INITIAL_WFN RT_RESTART + MAT_EXP BCH + PERIODIC .FALSE. + RTP_METHOD RTBSE + &PRINT + &MOMENTS + &END MOMENTS + &END PRINT + &END REAL_TIME_PROPAGATION + &SCF + ADDED_MOS -1 + EPS_SCF 1.0E-8 + SCF_GUESS ATOMIC + &DIAGONALIZATION ON + ALGORITHM STANDARD + &END DIAGONALIZATION + &MIXING + ALPHA 0.4 + METHOD BROYDEN_MIXING + NBROYDEN 4 + &END MIXING + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + ! Include G0W0 bandstructure correction + &PROPERTIES + &BANDSTRUCTURE + &GW + ! Simplest setting for quick checks + NUM_TIME_FREQ_POINTS 20 + &END GW + &END BANDSTRUCTURE + &END PROPERTIES + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + PERIODIC NONE + &END CELL + &COORD + H 0.0 0.0 0.0 + H 0.74 0.0 0.0 + &END COORD + &KIND H + BASIS_SET ORB SZV-MOLOPT-GTH + BASIS_SET RI_AUX RI_11_2_3_0_0_0_0_0_3.5e-03 + POTENTIAL GTH-PBE + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 597f85a04b..8dbd0fa6ec 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -365,3 +365,4 @@ Fist/regtest-deepmd deepmd QS/regtest-elpa elpa mpiranks==1||mpiranks%2==0 Fist/regtest-plumed2 plumed2 Fist/regtest-quip quip +QS/regtest-rtbse libint diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 4763398eca..932b2f6584 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -1,4 +1,4 @@ -114 +118 Total energy:!3 MD| Potential energy!5 Total energy \[eV\]:!4 @@ -113,6 +113,10 @@ HOMO-LUMO gap in evGW0 iteration 3 (eV)!8 HOMO-LUMO gap in evGW iteration 3 (eV)!8 BSE| 1 -TDA- !7 BSE| 1 -ABBA- !7 +MOMENTS_TRACE_RE| 0.10000000E+000!3 +MOMENTS_TRACE_IM| 0.10000000E+000!3 +MOMENTS_TRACE_RE| 0.20000000E+000!3 +POLARIZABILITY| 0.28663041E+002!4 # # these are the tests the can be selected for regtesting. # do regtest will grep for test_grep (first column) and look if the numeric value