import logging
from typing import Optional
import numpy as np
from scipy.linalg import block_diag
from .base_calculator import BaseCalculator
from ..response.excitations import Excitations
from ..response.utilities import Broadening, NoArtificialBroadening
from .. import CoupledSystem
from .utilities import get_dipole_dipole_tensor
from ..types import Array
logger = logging.getLogger(__name__)
[docs]
class CasidaCalculator(BaseCalculator):
"""Instances of this class enable the calculation of correlation energy
and spectrum of a coupled system, all polarizable objects in which have an
internal representation in the form of a Casida matrix.
The calculations are done for discrete excitations, and continuous spectra
calculated with this calculator use the given broadening.
Parameters
----------
coupled_system
Coupled system for which to carry out calculations.
broadening
Artificial broadening used for the continuous response;
defaults to no broadening.
name
Name of response.
Examples
--------
The following code snippet shows the calculation of the correlation energy
for a simple demo system:
>>> from strongcoca import CoupledSystem, PolarizableUnit
>>> from strongcoca.response import build_random_casida
>>> from strongcoca.calculators import CasidaCalculator
>>>
>>> # construct coupled system
>>> response = build_random_casida(n_states=2, name='Tiger', random_seed=42)
>>> pu1 = PolarizableUnit(response, [0, 0, 0], name='Cheetah')
>>> pu2 = PolarizableUnit(response, [3, 0, 0], name='Leopard')
>>> cs = CoupledSystem([pu1, pu2])
>>>
>>> # set up calculator
>>> calc = CasidaCalculator(cs)
>>>
>>> # ... and compute the correlation energy
>>> calc.get_correlation_energy()
-0.050798240...
"""
def __init__(self,
coupled_system: CoupledSystem,
broadening: Broadening = NoArtificialBroadening(),
name: str = 'CasidaCalculator') -> None:
super().__init__(coupled_system, broadening=broadening, name=name)
self._excitations: Optional[Excitations] = None
def _discard_data(self) -> None:
super()._discard_data()
self._excitations = None
@property
def excitations(self) -> Excitations:
"""Excitations as :class:`~strongcoca.response.excitations.Excitations`."""
self._verify_state()
if self._excitations is None:
self._excitations = self._build_casida_system()
return self._excitations
def _calculate_correlation_energy(self) -> float:
"""Returns the correlation energy of the coupled system in atomic units.
"""
logger.debug('Entering CasidaCalculator._calculate_correlation_energy')
energy = 0.5 * (np.sum(self.excitations._energies) - np.trace(self._D)) # type: float
return energy
def _build_casida_system(self) -> Excitations:
"""Compile the Casida matrices representing the coupled system.
"""
logger.debug('Entering CasidaCalculator._build_casida_system')
for k, pu in enumerate(self._coupled_system):
if not hasattr(pu.response, '_D_n') or \
not hasattr(pu.response, '_K_nn') or \
not hasattr(pu.response, '_mu_nv'):
raise ValueError(f'Polarizable unit #{k} has no internal Casida representation.')
# The 'type: ignore' statements on the following lines are necessary
# to apease mypy since it does not pick up that the loop above ensures
# that the components of the Casida representation are in fact available.
D = block_diag(*[np.diag(pu.response._D_n) # type: ignore
for pu in self._coupled_system])
K_blocks = []
for ii, pu1 in enumerate(self._coupled_system):
row = []
mu1 = np.array([pu1.orientation.apply(m)
for m in pu1.response._mu_nv]) # type: ignore
for jj, pu2 in enumerate(self._coupled_system):
if ii == jj:
K_sub = pu1.response._K_nn # type: ignore
else:
# todo: include cell and pbc
distance_vector = pu2._position - pu1._position
mu2 = np.array([pu2.orientation.apply(m)
for m in pu2.response._mu_nv]) # type: ignore
K_sub = np.dot(mu1, np.dot(get_dipole_dipole_tensor(distance_vector), mu2.T))
if ii > jj:
K_sub = K_sub.conj()
row.append(K_sub)
K_blocks.append(row)
K = np.block(K_blocks)
C = D ** 2 + 2 * np.dot(np.sqrt(D), np.dot(K, np.sqrt(D)))
U = block_diag(*[pu.response.U for pu in self._coupled_system]) # type: ignore
C = np.dot(U.conj().T, np.dot(C, U))
self._D = np.sqrt(np.diag(np.diag(C)))
eigval_I, eigvec_NI = np.linalg.eig(C)
srt_I = np.argsort(eigval_I)
eigval_I = eigval_I[srt_I]
eigvec_NI = eigvec_NI[:, srt_I]
# Calculate transition dipole moments
omega_I = np.sqrt(eigval_I)
F_NI = np.dot(U, eigvec_NI)
mu_vNN = []
for v in range(3):
mu_NN = block_diag(*[np.diag(pu.response._mu_nv[:, v]) # type: ignore
for pu in self._coupled_system])
mu_vNN.append(mu_NN)
mu_Iv = np.einsum('vNN,NN,NI,I->Iv', mu_vNN, np.sqrt(D), F_NI,
1 / np.sqrt(omega_I),
optimize=True)
return Excitations(omega_I, mu_Iv, broadening=self._broadening,
name=self._name, units='au')
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)
