from typing import Optional
import logging
import numpy as np
from .base import BaseResponse
from ..types import Array
from ..units import au_to_eV, eV_to_au, au_to_eA, eA_to_au
from ..utilities import ClassFormatter
from .excitations import Excitations
from .utilities import Broadening, NoArtificialBroadening
logger = logging.getLogger(__name__)
[docs]
class CasidaResponse(BaseResponse):
"""Objects of this class hold different representations of the response
functions for a Casida system.
Parameters
----------
D
Diagonal part of decomposition of Casida matrix
(initial-final energy difference of uncoupled single particle levels).
Units are eV by default, optionally atomic units (see :attr:`units`).
K
Coupling matrix of decomposition of Casida matrix
(electron-hole coupling).
Units are eV by default, optionally atomic units (see :attr:`units`).
mu
Transition dipoles.
Units are eÅ by default, optionally atomic units (see :attr:`units`).
occ
Occupation number differences, defaults to one, unitless.
broadening
Artificial broadening used for the continuous response;
defaults to no broadening.
units
`eVA` to specify D, K in eV and mu in eÅ or `au` to specify inputs
in atomic units.
This parameter determines whether conversion should be performed during
initialization and has no effect on instance methods and variables.
name
Name of response functioN.
Examples
--------
The following snippet illustrates how a Casida response object can be constructed.
Here, artificial Casida matrix components are used. In practice the Casida matrix
would usually be imported from, e.g., NWChem calculations:
>>> import numpy as np
>>> from strongcoca.response import CasidaResponse
>>>
>>> D = np.array([0.1, 0.3])
>>> K = np.zeros((2, 2))
>>> mu = np.array([[0.1, 0.0, 0.0], [0.0, 0.1, 0.0]])
>>> response = CasidaResponse(D, K, mu, name='Panda')
>>> print(response)
---------------------- CasidaResponse ----------------------
name : Panda
n_states : 2
D (eV) : [0.1 0.3]
K (eV) : [[0. 0.]
[0. 0.]]
mu (eÅ) : [[0.1 0. 0. ]
[0. 0.1 0. ]]
broadening : NoArtificialBroadening
atoms : None
"""
def __init__(self,
D: np.ndarray,
K: np.ndarray,
mu: np.ndarray,
broadening: Broadening = NoArtificialBroadening(),
occ: Optional[np.ndarray] = None,
units: str = 'eVA',
name: str = 'Casida') -> None:
logger.debug(f'Entering {self.__class__.__name__}.__init__')
super().__init__(broadening=broadening, pbc=False, name=name)
# TODO: explore way to generate these automatically based on type hints
if not isinstance(D, np.ndarray):
raise TypeError(f'D matrix must be provided as numpy array: {D}')
if not isinstance(K, np.ndarray):
raise TypeError(f'K matrix must be provided as numpy array: {K}')
if not isinstance(mu, np.ndarray):
raise TypeError(f'mu matrix must be provided as numpy array: {mu}')
if len(D.shape) != 1:
raise ValueError(f'D matrix has the wrong shape: {D.shape}')
n = len(D)
if K.shape != (n, n):
raise ValueError(f'K matrix must be {n}x{n}; actual shape: {K.shape}')
if mu.shape != (n, 3):
raise ValueError(f'mu matrix must be {n}x{3}; actual shape: {mu.shape}')
if occ is None:
occ = np.ones_like(D)
else:
# TODO: check that all code supports occupation numbers
raise NotImplementedError('Occupation numbers not implemented')
evals = np.linalg.eigvalsh(K)
if not np.all(evals > -1e-8):
raise ValueError(f'K matrix must be positive definite; eigenvalues: {evals}')
if units == 'eVA':
D = D * eV_to_au
K = K * eV_to_au
mu = mu * eA_to_au
elif units != 'au':
raise ValueError(f"units has to be 'eVA' or 'au', not '{units}'")
self._D_n = D
self._K_nn = K
self._mu_nv = mu
self._f_n = occ
# TODO: do not diagonalize and evaluate the following in init but
# use lazy evaluation
# See also issue #34
# compose and diagonalize Casida matrix
sq_fD_n = np.sqrt(self._f_n * self._D_n)
C_nn = np.diag(self._D_n)**2 + 2 * np.outer(sq_fD_n, sq_fD_n) * self._K_nn
eigval_I, F_nI = np.linalg.eig(C_nn)
# Calculate transition dipole moments
omega_I = np.sqrt(eigval_I)
mu_Iv = np.einsum('nv,n,nI,I->Iv', self._mu_nv, sq_fD_n,
F_nI, 1. / np.sqrt(omega_I),
optimize=True)
self._excitations = Excitations(omega_I, mu_Iv, broadening=self._broadening,
name=self._name, units='au')
self._eigvec_nI = F_nI
def __str__(self) -> str:
fmt = ClassFormatter(self, pad=15)
fmt.append_class_name()
fmt.append_attr('name')
fmt.append_attr('n_states')
fmt.append_attr('D', unit='eV')
fmt.append_attr('K', unit='eV')
fmt.append_attr('mu', unit='eÅ')
fmt.append_attr('broadening')
formula = self.atoms.get_chemical_formula() if self.atoms is not None else None
fmt.append_attr('atoms', formula)
return fmt.to_string()
@property
def n_states(self) -> int:
"""Number of states in Casida system."""
return len(self._D_n)
@property
def D(self) -> np.ndarray:
"""D matrix of Casida system in units of eV."""
return self._D_n * au_to_eV # type: ignore
@property
def K(self) -> np.ndarray:
"""K matrix of Casida system in units of eV."""
return self._K_nn * au_to_eV # type: ignore
@property
def mu(self) -> np.ndarray:
"""mu matrix of Casida system in units of eÅ."""
return self._mu_nv * au_to_eA # type: ignore
@property
def U(self) -> np.ndarray:
"""Unitless eigenvectors of Casida matrix."""
return self._eigvec_nI # type: ignore
@property
def excitations(self) -> Excitations:
"""Excitations as :class:`~strongcoca.response.excitations.Excitations`."""
return self._excitations
def _get_dynamic_polarizability(self, frequencies: Array) -> np.ndarray:
return self._excitations._get_dynamic_polarizability(frequencies)
def _get_dynamic_polarizability_imaginary_frequency(
self, frequencies: Array) -> np.ndarray:
return self._excitations._get_dynamic_polarizability_imaginary_frequency(
frequencies)
