diff --git a/adcc/ElectronicTransition.py b/adcc/ElectronicTransition.py
new file mode 100644
index 00000000..c78f1b2d
--- /dev/null
+++ b/adcc/ElectronicTransition.py
@@ -0,0 +1,210 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2020 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see .
+##
+## ---------------------------------------------------------------------
+import warnings
+import numpy as np
+
+from .misc import cached_property
+from .timings import timed_member_call
+from .visualisation import ExcitationSpectrum
+from .OneParticleOperator import product_trace
+
+from scipy import constants
+from matplotlib import pyplot as plt
+from .Excitation import mark_excitation_property
+
+
+class ElectronicTransition:
+ """
+ Documentation
+ """
+ @cached_property
+ @mark_excitation_property()
+ @timed_member_call(timer="_property_timer")
+ def transition_dipole_moment(self):
+ """List of transition dipole moments of all computed states"""
+ if self.property_method.level == 0:
+ warnings.warn("ADC(0) transition dipole moments are known to be "
+ "faulty in some cases.")
+ dipole_integrals = self.operators.electric_dipole
+ return np.array([
+ [product_trace(comp, tdm) for comp in dipole_integrals]
+ for tdm in self.transition_dm
+ ])
+
+ @cached_property
+ @mark_excitation_property()
+ @timed_member_call(timer="_property_timer")
+ def transition_dipole_moment_velocity(self):
+ """List of transition dipole moments in the
+ velocity gauge of all computed states"""
+ if self.property_method.level == 0:
+ warnings.warn("ADC(0) transition velocity dipole moments "
+ "are known to be faulty in some cases.")
+ dipole_integrals = self.operators.nabla
+ return np.array([
+ [product_trace(comp, tdm) for comp in dipole_integrals]
+ for tdm in self.transition_dm
+ ])
+
+ @cached_property
+ @mark_excitation_property()
+ @timed_member_call(timer="_property_timer")
+ def transition_magnetic_dipole_moment(self):
+ """List of transition magnetic dipole moments of all computed states"""
+ if self.property_method.level == 0:
+ warnings.warn("ADC(0) transition magnetic dipole moments "
+ "are known to be faulty in some cases.")
+ mag_dipole_integrals = self.operators.magnetic_dipole
+ return np.array([
+ [product_trace(comp, tdm) for comp in mag_dipole_integrals]
+ for tdm in self.transition_dm
+ ])
+
+ @cached_property
+ @mark_excitation_property()
+ def oscillator_strength(self):
+ """List of oscillator strengths of all computed states"""
+ return 2. / 3. * np.array([
+ np.linalg.norm(tdm)**2 * np.abs(ev)
+ for tdm, ev in zip(self.transition_dipole_moment,
+ self.excitation_energy)
+ ])
+
+ @cached_property
+ @mark_excitation_property()
+ def oscillator_strength_velocity(self):
+ """List of oscillator strengths in
+ velocity gauge of all computed states"""
+ return 2. / 3. * np.array([
+ np.linalg.norm(tdm)**2 / np.abs(ev)
+ for tdm, ev in zip(self.transition_dipole_moment_velocity,
+ self.excitation_energy)
+ ])
+
+ @cached_property
+ @mark_excitation_property()
+ def rotatory_strength(self):
+ """List of rotatory strengths of all computed states"""
+ return np.array([
+ np.dot(tdm, magmom) / ee
+ for tdm, magmom, ee in zip(self.transition_dipole_moment_velocity,
+ self.transition_magnetic_dipole_moment,
+ self.excitation_energy)
+ ])
+
+ @property
+ @mark_excitation_property()
+ def cross_section(self):
+ """List of one-photon absorption cross sections of all computed states"""
+ # TODO Source?
+ fine_structure = constants.fine_structure
+ fine_structure_au = 1 / fine_structure
+ prefac = 2.0 * np.pi ** 2 / fine_structure_au
+ return prefac * self.oscillator_strength
+
+ def plot_spectrum(self, broadening="lorentzian", xaxis="eV",
+ yaxis="cross_section", width=0.01, **kwargs):
+ """One-shot plotting function for the spectrum generated by all states
+ known to this class.
+
+ Makes use of the :class:`adcc.visualisation.ExcitationSpectrum` class
+ in order to generate and format the spectrum to be plotted, using
+ many sensible defaults.
+
+ Parameters
+ ----------
+ broadening : str or None or callable, optional
+ The broadening type to used for the computed excitations.
+ A value of None disables broadening any other value is passed
+ straight to
+ :func:`adcc.visualisation.ExcitationSpectrum.broaden_lines`.
+ xaxis : str
+ Energy unit to be used on the x-Axis. Options:
+ ["eV", "au", "nm", "cm-1"]
+ yaxis : str
+ Quantity to plot on the y-Axis. Options are "cross_section",
+ "osc_strength", "dipole" (plots norm of transition dipole),
+ "rotational_strength" (ECD spectrum with rotational strength)
+ width : float, optional
+ Gaussian broadening standard deviation or Lorentzian broadening
+ gamma parameter. The value should be given in atomic units
+ and will be converted to the unit of the energy axis.
+ """
+ if xaxis == "eV":
+ eV = constants.value("Hartree energy in eV")
+ energies = self.excitation_energy * eV
+ width = width * eV
+ xlabel = "Energy (eV)"
+ elif xaxis in ["au", "Hartree", "a.u."]:
+ energies = self.excitation_energy
+ xlabel = "Energy (au)"
+ elif xaxis == "nm":
+ hc = constants.h * constants.c
+ Eh = constants.value("Hartree energy")
+ energies = hc / (self.excitation_energy * Eh) * 1e9
+ xlabel = "Wavelength (nm)"
+ if broadening is not None and not callable(broadening):
+ raise ValueError("xaxis=nm and broadening enabled is "
+ "not supported.")
+ elif xaxis in ["cm-1", "cm^-1", "cm^{-1}"]:
+ towvn = constants.value("hartree-inverse meter relationship") / 100
+ energies = self.excitation_energy * towvn
+ width = width * towvn
+ xlabel = "Wavenumbers (cm^{-1})"
+ else:
+ raise ValueError("Unknown xaxis specifier: {}".format(xaxis))
+
+ if yaxis in ["osc", "osc_strength", "oscillator_strength", "f"]:
+ absorption = self.oscillator_strength
+ ylabel = "Oscillator strengths (au)"
+ elif yaxis in ["dipole", "dipole_norm", "μ"]:
+ absorption = np.linalg.norm(self.transition_dipole_moment, axis=1)
+ ylabel = "Modulus of transition dipole (au)"
+ elif yaxis in ["cross_section", "σ"]:
+ absorption = self.cross_section
+ ylabel = "Cross section (au)"
+ elif yaxis in ["rot", "rotational_strength", "rotatory_strength"]:
+ absorption = self.rotatory_strength
+ ylabel = "Rotatory strength (au)"
+ else:
+ raise ValueError("Unknown yaxis specifier: {}".format(yaxis))
+
+ sp = ExcitationSpectrum(energies, absorption)
+ sp.xlabel = xlabel
+ sp.ylabel = ylabel
+ if not broadening:
+ plots = sp.plot(style="discrete", **kwargs)
+ else:
+ kwdisc = kwargs.copy()
+ kwdisc.pop("label", "")
+ plots = sp.plot(style="discrete", **kwdisc)
+
+ kwargs.pop("color", "")
+ sp_broad = sp.broaden_lines(width, shape=broadening)
+ plots.extend(sp_broad.plot(color=plots[0].get_color(),
+ style="continuous", **kwargs))
+
+ if xaxis in ["nm"]:
+ # Invert x axis
+ plt.xlim(plt.xlim()[::-1])
+ return plots
diff --git a/adcc/ExcitedStates.py b/adcc/ExcitedStates.py
index 7bd666f2..a7148e70 100644
--- a/adcc/ExcitedStates.py
+++ b/adcc/ExcitedStates.py
@@ -30,17 +30,17 @@
from .AdcMethod import AdcMethod
from .FormatIndex import (FormatIndexAdcc, FormatIndexBase,
FormatIndexHfProvider, FormatIndexHomoLumo)
-from .visualisation import ExcitationSpectrum
from .OneParticleOperator import product_trace
from .FormatDominantElements import FormatDominantElements
from adcc import dot
from scipy import constants
-from matplotlib import pyplot as plt
from .solver.SolverStateBase import EigenSolverStateBase
from .Excitation import Excitation, mark_excitation_property
+from .ElectronicTransition import ElectronicTransition
+
class FormatExcitationVector:
def __init__(self, matrix, tolerance=0.01, index_format=None):
@@ -132,7 +132,7 @@ def format(self, vector):
return "\n".join(ret)
-class ExcitedStates:
+class ExcitedStates(ElectronicTransition):
def __init__(self, data, method=None, property_method=None,
excitation_energy_corrections={}):
"""Construct an ExcitedStates class from some data obtained
@@ -290,81 +290,6 @@ def transition_dm(self):
evec, self.matrix.intermediates)
for evec in self.excitation_vector]
- @cached_property
- @mark_excitation_property()
- @timed_member_call(timer="_property_timer")
- def transition_dipole_moment(self):
- """List of transition dipole moments of all computed states"""
- if self.property_method.level == 0:
- warnings.warn("ADC(0) transition dipole moments are known to be "
- "faulty in some cases.")
- dipole_integrals = self.operators.electric_dipole
- return np.array([
- [product_trace(comp, tdm) for comp in dipole_integrals]
- for tdm in self.transition_dm
- ])
-
- @cached_property
- @mark_excitation_property()
- @timed_member_call(timer="_property_timer")
- def transition_dipole_moment_velocity(self):
- """List of transition dipole moments in the
- velocity gauge of all computed states"""
- if self.property_method.level == 0:
- warnings.warn("ADC(0) transition velocity dipole moments "
- "are known to be faulty in some cases.")
- dipole_integrals = self.operators.nabla
- return np.array([
- [product_trace(comp, tdm) for comp in dipole_integrals]
- for tdm in self.transition_dm
- ])
-
- @cached_property
- @mark_excitation_property()
- @timed_member_call(timer="_property_timer")
- def transition_magnetic_dipole_moment(self):
- """List of transition magnetic dipole moments of all computed states"""
- if self.property_method.level == 0:
- warnings.warn("ADC(0) transition magnetic dipole moments "
- "are known to be faulty in some cases.")
- mag_dipole_integrals = self.operators.magnetic_dipole
- return np.array([
- [product_trace(comp, tdm) for comp in mag_dipole_integrals]
- for tdm in self.transition_dm
- ])
-
- @cached_property
- @mark_excitation_property()
- def oscillator_strength(self):
- """List of oscillator strengths of all computed states"""
- return 2. / 3. * np.array([
- np.linalg.norm(tdm)**2 * np.abs(ev)
- for tdm, ev in zip(self.transition_dipole_moment,
- self.excitation_energy)
- ])
-
- @cached_property
- @mark_excitation_property()
- def oscillator_strength_velocity(self):
- """List of oscillator strengths in
- velocity gauge of all computed states"""
- return 2. / 3. * np.array([
- np.linalg.norm(tdm)**2 / np.abs(ev)
- for tdm, ev in zip(self.transition_dipole_moment_velocity,
- self.excitation_energy)
- ])
-
- @cached_property
- @mark_excitation_property()
- def rotatory_strength(self):
- """List of rotatory strengths of all computed states"""
- return np.array([
- np.dot(tdm, magmom) / ee
- for tdm, magmom, ee in zip(self.transition_dipole_moment_velocity,
- self.transition_magnetic_dipole_moment,
- self.excitation_energy)
- ])
-
@cached_property
@mark_excitation_property(transform_to_ao=True)
@timed_member_call(timer="_property_timer")
@@ -398,103 +323,6 @@ def state_dipole_moment(self):
for ddm in self.state_diffdm
])
- @property
- @mark_excitation_property()
- def cross_section(self):
- """List of one-photon absorption cross sections of all computed states"""
- # TODO Source?
- fine_structure = constants.fine_structure
- fine_structure_au = 1 / fine_structure
- prefac = 2.0 * np.pi ** 2 / fine_structure_au
- return prefac * self.oscillator_strength
-
- def plot_spectrum(self, broadening="lorentzian", xaxis="eV",
- yaxis="cross_section", width=0.01, **kwargs):
- """One-shot plotting function for the spectrum generated by all states
- known to this class.
-
- Makes use of the :class:`adcc.visualisation.ExcitationSpectrum` class
- in order to generate and format the spectrum to be plotted, using
- many sensible defaults.
-
- Parameters
- ----------
- broadening : str or None or callable, optional
- The broadening type to used for the computed excitations.
- A value of None disables broadening any other value is passed
- straight to
- :func:`adcc.visualisation.ExcitationSpectrum.broaden_lines`.
- xaxis : str
- Energy unit to be used on the x-Axis. Options:
- ["eV", "au", "nm", "cm-1"]
- yaxis : str
- Quantity to plot on the y-Axis. Options are "cross_section",
- "osc_strength", "dipole" (plots norm of transition dipole),
- "rotational_strength" (ECD spectrum with rotational strength)
- width : float, optional
- Gaussian broadening standard deviation or Lorentzian broadening
- gamma parameter. The value should be given in atomic units
- and will be converted to the unit of the energy axis.
- """
- if xaxis == "eV":
- eV = constants.value("Hartree energy in eV")
- energies = self.excitation_energy * eV
- width = width * eV
- xlabel = "Energy (eV)"
- elif xaxis in ["au", "Hartree", "a.u."]:
- energies = self.excitation_energy
- xlabel = "Energy (au)"
- elif xaxis == "nm":
- hc = constants.h * constants.c
- Eh = constants.value("Hartree energy")
- energies = hc / (self.excitation_energy * Eh) * 1e9
- xlabel = "Wavelength (nm)"
- if broadening is not None and not callable(broadening):
- raise ValueError("xaxis=nm and broadening enabled is "
- "not supported.")
- elif xaxis in ["cm-1", "cm^-1", "cm^{-1}"]:
- towvn = constants.value("hartree-inverse meter relationship") / 100
- energies = self.excitation_energy * towvn
- width = width * towvn
- xlabel = "Wavenumbers (cm^{-1})"
- else:
- raise ValueError("Unknown xaxis specifier: {}".format(xaxis))
-
- if yaxis in ["osc", "osc_strength", "oscillator_strength", "f"]:
- absorption = self.oscillator_strength
- ylabel = "Oscillator strengths (au)"
- elif yaxis in ["dipole", "dipole_norm", "μ"]:
- absorption = np.linalg.norm(self.transition_dipole_moment, axis=1)
- ylabel = "Modulus of transition dipole (au)"
- elif yaxis in ["cross_section", "σ"]:
- absorption = self.cross_section
- ylabel = "Cross section (au)"
- elif yaxis in ["rot", "rotational_strength", "rotatory_strength"]:
- absorption = self.rotatory_strength
- ylabel = "Rotatory strength (au)"
- else:
- raise ValueError("Unknown yaxis specifier: {}".format(yaxis))
-
- sp = ExcitationSpectrum(energies, absorption)
- sp.xlabel = xlabel
- sp.ylabel = ylabel
- if not broadening:
- plots = sp.plot(style="discrete", **kwargs)
- else:
- kwdisc = kwargs.copy()
- kwdisc.pop("label", "")
- plots = sp.plot(style="discrete", **kwdisc)
-
- kwargs.pop("color", "")
- sp_broad = sp.broaden_lines(width, shape=broadening)
- plots.extend(sp_broad.plot(color=plots[0].get_color(),
- style="continuous", **kwargs))
-
- if xaxis in ["nm"]:
- # Invert x axis
- plt.xlim(plt.xlim()[::-1])
- return plots
-
def describe(self, oscillator_strengths=True, rotatory_strengths=False,
state_dipole_moments=False, transition_dipole_moments=False,
block_norms=True):
diff --git a/adcc/State2States.py b/adcc/State2States.py
new file mode 100644
index 00000000..668154d6
--- /dev/null
+++ b/adcc/State2States.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2020 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see .
+##
+## ---------------------------------------------------------------------
+import numpy as np
+
+from . import adc_pp
+from .misc import cached_property
+from .timings import timed_member_call
+
+from .Excitation import mark_excitation_property
+from .ElectronicTransition import ElectronicTransition
+
+
+class State2States(ElectronicTransition):
+ """
+ Documentation
+ """
+ def __init__(self, parent_state, initial=0):
+ self.parent_state = parent_state
+ self.reference_state = self.parent_state.reference_state
+ self.ground_state = self.parent_state.ground_state
+ self.matrix = self.parent_state.matrix
+ self.property_method = self.parent_state.property_method
+ self.operators = self.parent_state.operators
+ self.initial = initial
+ # TODO: describe?!
+
+ @property
+ def excitation_energy(self):
+ return np.array([
+ e.excitation_energy
+ - self.parent_state.excitation_energy[self.initial]
+ for e in self.parent_state.excitations if e.index > self.initial
+ ])
+
+ @cached_property
+ @mark_excitation_property(transform_to_ao=True)
+ @timed_member_call(timer="_property_timer")
+ def transition_dm(self):
+ """
+ List of transition density matrices from
+ initial state to final state/s
+ """
+ return [
+ adc_pp.state2state_transition_dm(
+ self.property_method, self.ground_state,
+ self.parent_state.excitation_vector[self.initial],
+ e.excitation_vector, self.matrix.intermediates
+ )
+ for e in self.parent_state.excitations
+ if e.index > self.initial
+ ]
diff --git a/adcc/__init__.py b/adcc/__init__.py
index be7d2dec..d0215d01 100644
--- a/adcc/__init__.py
+++ b/adcc/__init__.py
@@ -36,6 +36,7 @@
from .memory_pool import memory_pool
from .AdcBlockView import AdcBlockView
from .ExcitedStates import ExcitedStates
+from .State2States import State2States
from .DataHfProvider import DataHfProvider, DictHfProvider
from .ReferenceState import ReferenceState
from .caching_policy import DefaultCachingPolicy, GatherStatisticsPolicy
@@ -55,7 +56,7 @@
"lincomb", "nosym_like", "ones_like", "transpose",
"linear_combination", "zeros_like", "direct_sum",
"memory_pool", "set_n_threads", "get_n_threads", "AmplitudeVector",
- "HartreeFockProvider", "ExcitedStates",
+ "HartreeFockProvider", "ExcitedStates", "State2States",
"Tensor", "DictHfProvider", "DataHfProvider", "OneParticleOperator",
"guesses_singlet", "guesses_triplet", "guesses_any",
"guess_symmetries", "guesses_spin_flip", "guess_zero",
diff --git a/adcc/adc_pp/state2state_transition_dm.py b/adcc/adc_pp/state2state_transition_dm.py
index a4a3c4c1..b0c8eeac 100644
--- a/adcc/adc_pp/state2state_transition_dm.py
+++ b/adcc/adc_pp/state2state_transition_dm.py
@@ -20,12 +20,118 @@
## along with adcc. If not, see .
##
## ---------------------------------------------------------------------
+from adcc import block as b
from adcc.AdcMethod import AdcMethod
+from adcc.functions import einsum
from adcc.AmplitudeVector import AmplitudeVector
+from adcc.OneParticleOperator import OneParticleOperator
import libadcc
-DISPATCH = {} # None implemented
+
+def s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates):
+ x1 = amplitude_l["s"]
+ y1 = amplitude_r["s"]
+
+ dm = OneParticleOperator(mp, is_symmetric=False)
+ dm[b.oo] = -einsum('ja,ia->ij', x1, y1)
+ dm[b.vv] = einsum('ia,ib->ab', x1, y1)
+ return dm
+
+
+def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
+ dm = s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates)
+
+ x1 = amplitude_l["s"]
+ y1 = amplitude_r["s"]
+ x2 = amplitude_l["d"]
+ y2 = amplitude_r["d"]
+
+ t2 = mp.t2(b.oovv)
+ p0_ov = mp.mp2_diffdm[b.ov]
+ p0_oo = mp.mp2_diffdm[b.oo]
+ p0_vv = mp.mp2_diffdm[b.vv]
+ p1_oo = dm[b.oo].evaluate() # ADC(1) tdm
+ p1_vv = dm[b.vv].evaluate() # ADC(1) tdm
+
+ rx1 = einsum('ijab,jb->ia', t2, x1)
+ ry1 = einsum('ijab,jb->ia', t2, y1)
+
+ dm[b.oo] = (
+ p1_oo - 2.0 * einsum('ikab,jkab->ij', y2, x2) + (
+ + 0.5 * (
+ einsum('ik,kj->ij', p1_oo, p0_oo)
+ + einsum('ik,kj->ij', p0_oo, p1_oo)
+ )
+ - einsum(
+ 'ikcd,jkcd->ij', t2,
+ einsum('lk,jlcd->jkcd', 0.5 * p1_oo, t2)
+ - einsum('jkcb,db->jkcd', t2, p1_vv)
+ )
+ - 0.5 * (
+ einsum('ia,ja->ij', y1, einsum('jkac,kc->ja', t2, rx1))
+ + einsum('ikac,kc,ja->ij', t2, ry1, x1)
+ )
+ - einsum('ia,ja->ij', rx1, ry1)
+ )
+ )
+ dm[b.vv] = (
+ p1_vv + 2.0 * einsum('ijac,ijbc->ab', x2, y2)
+ + (
+ - 0.5 * (
+ einsum("ac,cb->ab", p1_vv, p0_vv)
+ + einsum("ac,cb->ab", p0_vv, p1_vv)
+ )
+ - einsum(
+ "klbc,klac->ab", t2,
+ 0.5 * einsum("klad,cd->klac", t2, p1_vv)
+ - einsum("jk,jlac->klac", p1_oo, t2)
+ )
+ + 0.5 * (
+ einsum("ikac,kc,ib->ab", t2, rx1, y1)
+ + einsum(
+ "ia,ib->ab", x1,
+ einsum("ikbc,kc->ib", t2, ry1)
+ )
+ )
+ + einsum("ia,ib->ab", ry1, rx1)
+ )
+ )
+
+ p1_ov = -2.0 * einsum("jb,ijab->ia", x1, y2)
+ p1_vo = -2.0 * einsum("ijab,jb->ai", x2, y1)
+
+ dm[b.ov] = (
+ p1_ov
+ - einsum("ijab,bj->ia", t2, p1_vo)
+ - einsum("ib,ba->ia", p0_ov, p1_vv)
+ + einsum("ij,ja->ia", p1_oo, p0_ov)
+ - einsum(
+ "ib,ab->ia", y1,
+ einsum("klca,klcb->ab", t2, x2)
+ )
+ - einsum("ikcd,jkcd,ja->ia", t2, x2, y1)
+ )
+ dm[b.vo] = (
+ p1_vo
+ - einsum("ijab,jb->ai", t2, p1_ov)
+ - einsum("ib,ab->ai", p0_ov, p1_vv)
+ + einsum("ji,ja->ai", p1_oo, p0_ov)
+ - einsum(
+ "ib,ab->ai", x1,
+ einsum("klca,klcb->ab", t2, y2)
+ )
+ - einsum("ikcd,jkcd,ja->ai", t2, y2, x1)
+ )
+ return dm
+
+
+# Ref: https://doi.org/10.1080/00268976.2013.859313
+DISPATCH = {"adc0": s2s_tdm_adc0,
+ "adc1": s2s_tdm_adc0, # same as ADC(0)
+ "adc2": s2s_tdm_adc2,
+ "adc2x": s2s_tdm_adc2, # same as ADC(2)
+ }
def state2state_transition_dm(method, ground_state, amplitude_from,
diff --git a/adcc/test_properties.py b/adcc/test_properties.py
index 8728478f..f982fd7d 100644
--- a/adcc/test_properties.py
+++ b/adcc/test_properties.py
@@ -23,13 +23,15 @@
import unittest
import numpy as np
-from .misc import assert_allclose_signfix
-from .test_state_densities import Runners
-
from numpy.testing import assert_allclose
+
+from adcc.State2States import State2States
from adcc.testdata.cache import cache
-from pytest import approx
+from .misc import assert_allclose_signfix
+from .test_state_densities import Runners
+
+from pytest import approx, skip
class TestTransitionDipoleMoments(unittest.TestCase, Runners):
@@ -81,3 +83,25 @@ def base_test(self, system, method, kind):
ref = refdata[method][kind]
n_ref = len(state.excitation_vector)
assert_allclose(res_dms, ref["state_dipole_moments"][:n_ref], atol=1e-4)
+
+
+class TestState2StateTransitionDipoleMoments(unittest.TestCase, Runners):
+ def base_test(self, system, method, kind):
+ method = method.replace("_", "-")
+ if "cvs" in method:
+ skip("State-to-state transition dms not yet implemented for CVS.")
+
+ refdata = cache.reference_data[system]
+ state = cache.adc_states[system][method][kind]
+
+ state_to_state = refdata[method][kind]["state_to_state"]
+ refevals = refdata[method][kind]["eigenvalues"]
+ for i, exci in enumerate(state.excitations):
+ assert exci.excitation_energy == refevals[i]
+ fromi_ref = state_to_state[f"from_{i}"]["transition_dipole_moments"]
+
+ state2state = State2States(state, initial=i)
+ for ii, j in enumerate(range(i + 1, state.size)):
+ assert state.excitation_energy[j] == refevals[j]
+ assert_allclose_signfix(state2state.transition_dipole_moment[ii],
+ fromi_ref[ii], atol=1e-4)
diff --git a/adcc/test_state_densities.py b/adcc/test_state_densities.py
index e48341ba..9d14d71a 100644
--- a/adcc/test_state_densities.py
+++ b/adcc/test_state_densities.py
@@ -21,11 +21,13 @@
##
## ---------------------------------------------------------------------
import unittest
+import numpy as np
-from .misc import expand_test_templates
+from adcc.State2States import State2States
from adcc.testdata.cache import cache
-from pytest import approx
+from .misc import expand_test_templates
+from pytest import approx, skip
# The methods to test
basemethods = ["adc0", "adc1", "adc2", "adc2x", "adc3"]
@@ -119,10 +121,8 @@ def base_test(self, system, method, kind):
assert state.excitation_energy[i] == refevals[i]
dm_ao_a, dm_ao_b = state.state_diffdm[i].to_ao_basis()
- dm_ao_a = dm_ao_a.to_ndarray()
- dm_ao_b = dm_ao_b.to_ndarray()
- assert dm_ao_a == approx(refdens_a[i])
- assert dm_ao_b == approx(refdens_b[i])
+ assert dm_ao_a.to_ndarray() == approx(refdens_a[i])
+ assert dm_ao_b.to_ndarray() == approx(refdens_b[i])
class TestStateGroundToExcitedTdm(unittest.TestCase, Runners):
@@ -140,9 +140,35 @@ def base_test(self, system, method, kind):
# comparing reference and computed
assert state.excitation_energy[i] == refevals[i]
- tdms = state.transition_dm[i].to_ao_basis()
- dm_ao_a, dm_ao_b = tdms
- dm_ao_a = dm_ao_a.to_ndarray()
- dm_ao_b = dm_ao_b.to_ndarray()
- assert dm_ao_a == approx(refdens_a[i])
- assert dm_ao_b == approx(refdens_b[i])
+ dm_ao_a, dm_ao_b = state.transition_dm[i].to_ao_basis()
+ assert dm_ao_a.to_ndarray() == approx(refdens_a[i])
+ assert dm_ao_b.to_ndarray() == approx(refdens_b[i])
+
+
+class TestStateExcitedToExcitedTdm(unittest.TestCase, Runners):
+ def base_test(self, system, method, kind):
+ method = method.replace("_", "-")
+ if "cvs" in method:
+ skip("State-to-state transition dms not yet implemented for CVS.")
+ refdata = cache.reference_data[system]
+ state = cache.adc_states[system][method][kind]
+ state_to_state = refdata[method][kind]["state_to_state"]
+ refevals = refdata[method][kind]["eigenvalues"]
+
+ for i, exci in enumerate(state.excitations):
+ # Check that we are talking about the same state when
+ # comparing reference and computed
+ assert exci.excitation_energy == refevals[i]
+ fromi_ref_a = state_to_state[f"from_{i}"]["state_to_excited_tdm_bb_a"]
+ fromi_ref_b = state_to_state[f"from_{i}"]["state_to_excited_tdm_bb_b"]
+
+ state2state = State2States(state, initial=i)
+ for ii, j in enumerate(range(i + 1, state.size)):
+ assert state.excitation_energy[j] == refevals[j]
+ ee_ref = refevals[j] - refevals[i]
+ assert state2state.excitation_energy[ii] == ee_ref
+ dm_ao_a, dm_ao_b = state2state.transition_dm[ii].to_ao_basis()
+ np.testing.assert_allclose(fromi_ref_a[ii],
+ dm_ao_a.to_ndarray().T, atol=1e-4)
+ np.testing.assert_allclose(fromi_ref_b[ii],
+ dm_ao_b.to_ndarray().T, atol=1e-4)
diff --git a/examples/water/tpa.py b/examples/water/tpa.py
index 9ac78ea8..bf2a975b 100644
--- a/examples/water/tpa.py
+++ b/examples/water/tpa.py
@@ -11,8 +11,8 @@
from adcc.AmplitudeVector import AmplitudeVector
from adcc.solver import IndexSymmetrisation
from adcc.solver.conjugate_gradient import conjugate_gradient, default_print
-from adcc.modified_transition_moments import compute_modified_transition_moments
-from adcc.state_densities import compute_state2state_optdm
+from adcc.adc_pp.modified_transition_moments import modified_transition_moments
+from adcc.adc_pp.state2state_transition_dm import state2state_transition_dm
from adcc.OneParticleOperator import product_trace
@@ -47,10 +47,7 @@ def __matmul__(self, other):
# setup modified transition moments
dips = state.reference_state.operators.electric_dipole
-rhss = [
- compute_modified_transition_moments(state.matrix, dip, "adc2")
- for dip in dips
-]
+rhss = modified_transition_moments("adc2", state.ground_state, dips)
S = np.zeros((len(state.excitation_energy), 3, 3))
for f, ee in enumerate(state.excitation_energy):
@@ -72,11 +69,11 @@ def __matmul__(self, other):
response.append(res)
for mu in range(3):
for nu in range(mu, 3):
- tdm_mu_f = compute_state2state_optdm(
+ tdm_mu_f = state2state_transition_dm(
"adc2", matrix.ground_state, response[mu].solution,
state.excitation_vector[f]
)
- tdm_nu_f = compute_state2state_optdm(
+ tdm_nu_f = state2state_transition_dm(
"adc2", matrix.ground_state, response[nu].solution,
state.excitation_vector[f]
)