Source code for strongcoca.response.time_domain_casida

import logging
from typing import Optional, Tuple
import numpy as np
from . import CasidaResponse, TimeDomainResponse
from ..units import fs_to_au


logger = logging.getLogger(__name__)


def _check_propagation_parameters(eigenvals: np.ndarray,
                                  dt: Optional[float] = None,
                                  n_steps: Optional[int] = None) -> Tuple[float, int]:
    """Checks propagation parameters based on eigenvalues.

    Parameters
    ----------
    dt
        Time step in atomic units.

        The time step is used to propagate the Casida system and
        generate the dipole tensor in the time domain; by default the time
        step is set based on the highest frequency of the eigenvalue spectrum.
    n_steps
        Number of time steps to propagate the Casida system and generate the
        dipole tensor in the time domain; by default the value is set
        based on the lowest frequency of the eigenvalue spectrum.

    Returns
    -------
    dt
        Time step in atomic units (the same as given or recommended).
    n_steps
        Number of time steps (the same as given or recommended).
    """

    # determine suitable integration parameters and
    # check sensibility of user-provided values (if applicable)
    # - time step
    dt_recommended = 1 / sorted(eigenvals)[-1]
    logger.debug(f'dt_recommended: {dt_recommended} a.u.')
    if dt is None:
        dt = dt_recommended
    else:
        if dt > 2 * dt_recommended or dt < 0.5 * dt_recommended:
            logger.warning(f'dt ({dt} a.u.) deviates strongly from the'
                           f' recommended value ({dt_recommended} a.u.).')
    logger.debug(f'dt_actual: {dt} a.u.')
    # - maximum time and number of time steps
    t_max_recommended = 100 / sorted(eigenvals)[0]
    logger.debug(f't_max_recommended: {t_max_recommended} a.u.')
    if n_steps is None:
        n_steps = int(t_max_recommended / dt)
    else:
        t_max = dt * (n_steps - 1)
        if t_max < 0.5 * t_max_recommended or t_max > 2 * t_max_recommended:
            logger.warning(f't_max ({t_max} a.u.) deviates strongly from the'
                           f' recommended value ({t_max_recommended} a.u.).')
    logger.debug(f'n_steps: {n_steps}')

    return dt, n_steps




def build_time_domain_response_from_casida(
        casida_response: CasidaResponse,
        dt: Optional[float] = None,
        n_steps: Optional[int] = None,
        units: str = 'fs',
        name: str = 'PropagatedCasida') -> TimeDomainResponse:
    """Builds a :class:`~strongcoca.response.TimeDomainResponse` instance
    from an existing Casida representation via time propagation.

    Parameters
    ----------
    casida_response
        Response in Casida representation.
    dt
        Time step; units are fs by default, optionally atomic units (see :attr:`units`).
    n_steps
        Number of time steps for which to propagate the equations of motion.
    units
        `fs` to specify dt in fs or `au` to specify `dt` in atomic units.
    """

    if units == 'fs':
        if dt is not None:
            dt = dt * fs_to_au
    elif units != 'au':
        raise ValueError(f"units has to be 'fs' or 'au', not '{units}'")

    # TODO: The argument below should probably be just `energies`. Test this. See #41.
    dt, n_steps = _check_propagation_parameters(
            casida_response.excitations._energies**2, dt, n_steps)

    # Matrices
    D_n = casida_response._D_n
    K_nn = casida_response._K_nn
    mu_nv = casida_response._mu_nv

    # Propagate
    time_t = np.linspace(0, dt * (n_steps - 1), n_steps)
    dm_tvv = np.zeros((n_steps, 3, 3))

    # Integrate equations of motion
    Re_rho_nv = 0 * mu_nv
    Im_rho_nv = -mu_nv
    for t in range(0, n_steps):
        # Evaluate dipole moment
        dm_tvv[t] = 2 * np.dot(Re_rho_nv.T, mu_nv)

        # Symplectic propagator for real-time TDDFT
        Im_rho_nv += dt * D_n[:, np.newaxis] * Re_rho_nv + 2 * dt * np.dot(K_nn, Re_rho_nv)
        Re_rho_nv -= dt * D_n[:, np.newaxis] * Im_rho_nv

    # Broadening
    broadening = casida_response.broadening
    return TimeDomainResponse(time_t, dm_tvv, broadening=broadening,
                              name=name, units='au')