diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py
index 016475b4..6776f81b 100644
--- a/adcc/AdcMatrix.py
+++ b/adcc/AdcMatrix.py
@@ -75,7 +75,6 @@ class AdcMatrixlike:
class AdcMatrix(AdcMatrixlike):
# Default perturbation-theory orders for the matrix blocks (== standard ADC-PP).
default_block_orders = {
- # ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None),
"adc0": dict(ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None), # noqa: E501
"adc1": dict(ph_ph=1, ph_pphh=None, pphh_ph=None, pphh_pphh=None), # noqa: E501
"adc2": dict(ph_ph=2, ph_pphh=1, pphh_ph=1, pphh_pphh=0), # noqa: E501
@@ -87,7 +86,6 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
diagonal_precomputed=None):
"""
Initialise an ADC matrix.
-
Parameters
----------
method : str or AdcMethod
@@ -130,6 +128,14 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
self.ndim = 2
self.extra_terms = []
+ if self.reference_state.is_qed:
+ if method.base_method.name in ["adc2x", "adc3"] or self.is_core_valence_separated: # noqa: E501
+ raise NotImplementedError("Neither adc2x and adc3 nor cvs methods "
+ "are implemented for QED-ADC")
+ elif method.base_method.name == "adc2" and not self.reference_state.qed_hf: # noqa: E501
+ raise NotImplementedError("QED-ADC(2) is only available for a "
+ "QED-HF reference")
+
self.intermediates = intermediates
if self.intermediates is None:
self.intermediates = Intermediates(self.ground_state)
@@ -137,6 +143,8 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
# Determine orders of PT in the blocks
if block_orders is None:
block_orders = self.default_block_orders[method.base_method.name]
+ if self.reference_state.is_qed and not self.reference_state.approx:
+ block_orders["ph_gs"] = block_orders["ph_ph"]
else:
tmp_orders = self.default_block_orders[method.base_method.name].copy()
tmp_orders.update(block_orders)
@@ -144,7 +152,7 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
# Sanity checks on block_orders
for block in block_orders.keys():
- if block not in ("ph_ph", "ph_pphh", "pphh_ph", "pphh_pphh"):
+ if block not in ("ph_gs", "ph_ph", "ph_pphh", "pphh_ph", "pphh_pphh"):
raise ValueError(f"Invalid block order key: {block}")
if block_orders["ph_pphh"] != block_orders["pphh_ph"]:
raise ValueError("ph_pphh and pphh_ph should always have "
@@ -154,19 +162,53 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
raise ValueError("pphh_pphh cannot be None if ph_pphh isn't.")
self.block_orders = block_orders
+ self.qed_dispatch_dict = {
+ "elec": "",
+ "elec_couple": "_couple", "phot_couple": "_phot_couple",
+ "phot": "_phot", "phot2": "_phot2",
+ "elec_couple_inner": "_couple_inner",
+ "elec_couple_edge": "_couple_edge",
+ "phot_couple_inner": "_phot_couple_inner",
+ "phot_couple_edge": "_phot_couple_edge"
+ }
+
+ def get_pp_blocks(disp_str):
+ """
+ Extraction of blocks from the adc_pp matrix
+ ----------
+ disp_str : string
+ key specifying block to return
+ """
+ return { # TODO Rename to self.block in 0.16.0
+ block: ppmatrix.block(self.ground_state, block.split("_"),
+ order=str(order) +\
+ self.qed_dispatch_dict[disp_str],
+ intermediates=self.intermediates,
+ variant=variant)
+ for block, order in block_orders.items() if order is not None
+ }
+
# Build the blocks and diagonals
+
with self.timer.record("build"):
variant = None
- if self.is_core_valence_separated:
+ if method.is_core_valence_separated:
variant = "cvs"
- blocks = {
- block: ppmatrix.block(self.ground_state, block.split("_"),
- order=order, intermediates=self.intermediates,
- variant=variant)
- for block, order in self.block_orders.items() if order is not None
- }
- # TODO Rename to self.block in 0.16.0
- self.blocks_ph = {bl: blocks[bl].apply for bl in blocks}
+ # Build full QED-matrix
+ if self.reference_state.is_qed and not self.reference_state.approx:
+ blocks = {
+ bl + "_" + key: get_pp_blocks(key)[bl]
+ for bl, order in block_orders.items()
+ if order is not None
+ for key in self.qed_dispatch_dict
+ }
+ else: # Build "standard" ADC-matrix
+ blocks = {
+ bl: get_pp_blocks("elec")[bl]
+ for bl, order in block_orders.items()
+ if order is not None
+ }
+
if diagonal_precomputed:
self.__diagonal = diagonal_precomputed
else:
@@ -175,9 +217,11 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
self.__diagonal.evaluate()
self.__init_space_data(self.__diagonal)
+ # TODO Rename to self.block in 0.16.0
+ self.blocks_ph = {bl: blocks[bl].apply for bl in blocks}
+
def __iadd__(self, other):
"""In-place addition of an :py:class:`AdcExtraTerm`
-
Parameters
----------
other : AdcExtraTerm
@@ -205,12 +249,10 @@ def patched_apply(ampl, original=orig_app, other=other_app):
def __add__(self, other):
"""Addition of an :py:class:`AdcExtraTerm`, creating
a copy of self and adding the term to the new matrix
-
Parameters
----------
other : AdcExtraTerm
the extra term to be added
-
Returns
-------
AdcMatrix
@@ -233,10 +275,17 @@ def __init_space_data(self, diagonal):
self.axis_spaces = {}
self.axis_lengths = {}
for block in diagonal.blocks_ph:
- self.axis_spaces[block] = getattr(diagonal, block).subspaces
- self.axis_lengths[block] = np.prod([
- self.mospaces.n_orbs(sp) for sp in self.axis_spaces[block]
- ])
+ if "gs" in block:
+ # Either include g1 in whole libadcc backend, or use this
+ # approach for now, which is only required for functionalities,
+ # which should not be used with the full qed matrix yet anyway
+ self.axis_spaces[block] = ['g1']
+ self.axis_lengths[block] = 1
+ else:
+ self.axis_spaces[block] = getattr(diagonal, block).subspaces
+ self.axis_lengths[block] = np.prod([
+ self.mospaces.n_orbs(sp) for sp in self.axis_spaces[block]
+ ])
self.shape = (sum(self.axis_lengths.values()),
sum(self.axis_lengths.values()))
@@ -339,7 +388,7 @@ def __matmul__(self, other):
if isinstance(other, AmplitudeVector):
return self.matvec(other)
if isinstance(other, list):
- if all(isinstance(elem, AmplitudeVector) for elem in other):
+ if all(isinstance(elem, AmplitudeVector) for elem in other): # noqa: E501
return [self.matvec(ov) for ov in other]
return NotImplemented
@@ -370,11 +419,9 @@ def construct_symmetrisation_for_blocks(self):
applied to relevant blocks of an AmplitudeVector in order
to symmetrise it to the right symmetry in order to be used
with the various matrix-vector-products of this function.
-
Most importantly the returned functions antisymmetrise
the occupied and virtual parts of the doubles parts
if this is sensible for the method behind this adcmatrix.
-
Returns a dictionary block identifier -> function
"""
ret = {}
@@ -394,7 +441,6 @@ def dense_basis(self, axis_blocks=None, ordering="adcc"):
"""
Return the list of indices and their values
of the dense basis representation
-
ordering: adcc, spin, spatial
"""
ret = []
@@ -484,19 +530,15 @@ def to_ndarray(self, out=None):
Return the ADC matrix object as a dense numpy array. Converts the sparse
internal representation of the ADC matrix to a dense matrix and return
as a numpy array.
-
Notes
-----
-
This method is only intended to be used for debugging and
visualisation purposes as it involves computing a large amount of
matrix-vector products and the returned array consumes a considerable
amount of memory.
-
The resulting matrix has no spin symmetry imposed, which means that
its eigenspectrum may contain non-physical excitations (e.g. with linear
combinations of α->β and α->α components in the excitation vector).
-
This function has not been sufficiently tested to be considered stable.
"""
# TODO Update to ph / pphh
@@ -589,7 +631,6 @@ def __init__(self, matrix, shift=0.0):
"""
Initialise a shifted ADC matrix. Applying this class to a vector ``v``
represents an efficient version of ``matrix @ v + shift * v``.
-
Parameters
----------
matrix : AdcMatrix
@@ -638,18 +679,15 @@ def __init__(self, matrix, excitation_blocks, core_orbitals=None,
Initialise a projected ADC matrix, i.e. represents the expression
``P @ M @ P`` where ``P`` is a projector onto a subset of
``excitation_blocks``.
-
The ``excitation_blocks`` are defined by partitioning the ``o1`` occupied
and ``v1`` virtual space of the ``matrix.mospaces`` into a core-occupied
``c``, valence-occupied ``o``, inner-virtual ``v`` and outer-virtual ``w``.
This matrix will only keep selected blocks in the amplitudes non-zero, which
are selected in the ``excitation_blocks`` list
(e.g. ``["cv", "ccvv", "ocvv"]``).
-
For details on the option how to select the spaces, see the documentation
in :py:`adcc.ReferenceState.__init__` (``outer_virtuals`` follows the same
rules as ``frozen_virtuals``).
-
Parameters
----------
matrix : AdcMatrix
diff --git a/adcc/ElectronicTransition.py b/adcc/ElectronicTransition.py
index c2f7e4a0..a147d724 100644
--- a/adcc/ElectronicTransition.py
+++ b/adcc/ElectronicTransition.py
@@ -28,6 +28,10 @@
from .visualisation import ExcitationSpectrum
from .OneParticleOperator import product_trace
from .AdcMethod import AdcMethod
+from adcc.functions import einsum
+from adcc import block as b
+from adcc.adc_pp.state2state_transition_dm import state2state_transition_dm
+from adcc.adc_pp.transition_dm import transition_dm
from scipy import constants
from .Excitation import mark_excitation_property
@@ -165,6 +169,122 @@ def transition_dipole_moment(self):
for tdm in self.transition_dm
])
+ @cached_property
+ @mark_excitation_property()
+ @timed_member_call(timer="_property_timer")
+ def transition_dipole_moments_qed(self):
+ """
+ List of transition dipole moments of all computed states
+ to build the QED-matrix in the basis of the diagonal
+ purely electric subblock
+ """
+ if self.reference_state.approx:
+
+ dipole_integrals = self.operators.electric_dipole
+
+ def tdm(i, prop_level):
+ self.ground_state.tdm_contribution = prop_level
+ return transition_dm(self.method, self.ground_state,
+ self.excitation_vector[i])
+
+ prop_level = "adc" + str(self.property_method.level - 1)
+ return np.array([
+ [product_trace(comp, tdm(i, prop_level))
+ for comp in dipole_integrals]
+ for i in np.arange(len(self.excitation_energy))
+ ])
+ else:
+ return ("transition_dipole_moments_qed are only calculated,"
+ "if reference_state contains 'approx' attribute")
+
+ @cached_property
+ @mark_excitation_property()
+ @timed_member_call(timer="_property_timer")
+ def s2s_dipole_moments_qed(self):
+ """
+ List of s2s transition dipole moments of all computed states
+ to build the QED-matrix in the basis of the diagonal
+ purely electric subblock
+ """
+ if self.reference_state.approx:
+ dipole_integrals = self.operators.electric_dipole
+ print("note, that only the z coordinate of the "
+ "dipole integrals is calculated")
+ n_states = len(self.excitation_energy)
+
+ def s2s(i, f, s2s_contribution):
+ self.ground_state.s2s_contribution = s2s_contribution
+ vec = self.excitation_vector
+ return state2state_transition_dm(self.method, self.ground_state,
+ vec[i], vec[f])
+
+ def final_block(name):
+ return np.array([[[product_trace(comp, s2s(i, j, name))
+ for comp in dipole_integrals]
+ for j in np.arange(n_states)]
+ for i in np.arange(n_states)])
+
+ block_dict = {"qed_adc1_off_diag": final_block("adc1")}
+
+ if self.method.name == "adc2":
+ keys = ("qed_adc2_diag", "qed_adc2_edge_couple",
+ "qed_adc2_edge_phot_couple", "qed_adc2_ph_pphh",
+ "qed_adc2_pphh_ph")
+ for key in keys:
+ block_dict[key] = final_block(key)
+ return block_dict
+ else:
+ return ("s2s_dipole_moments_qed are only calculated,"
+ "if reference_state contains 'approx' attribute")
+
+ @cached_property
+ @mark_excitation_property()
+ @timed_member_call(timer="_property_timer")
+ def qed_second_order_ph_ph_couplings(self):
+ """
+ List of blocks containing the expectation value of the perturbation
+ of the Hamiltonian for all computed states required
+ to build the QED-matrix in the basis of the diagonal
+ purely electric subblock
+ """
+ if self.reference_state.approx:
+ qed_t1 = self.ground_state.qed_t1(b.ov)
+
+ def couple(qed_t1, ul, ur):
+ return {
+ b.ooov: einsum("kc,ia,ja->kjic", qed_t1, ul, ur)
+ + einsum("ka,ia,jb->jkib", qed_t1, ul, ur),
+ b.ovvv: einsum("kc,ia,ib->kacb", qed_t1, ul, ur)
+ + einsum("ic,ia,jb->jabc", qed_t1, ul, ur)
+ }
+
+ def phot_couple(qed_t1, ul, ur):
+ return {
+ b.ooov: einsum("kc,ia,ja->kijc", qed_t1, ul, ur)
+ + einsum("kb,ia,jb->ikja", qed_t1, ul, ur),
+ b.ovvv: einsum("kc,ia,ib->kbca", qed_t1, ul, ur)
+ + einsum("jc,ia,jb->ibac", qed_t1, ul, ur)
+ }
+
+ def prod_sum(hf, two_p_op):
+ return (einsum("ijka,ijka->", hf.ooov, two_p_op[b.ooov])
+ + einsum("iabc,iabc->", hf.ovvv, two_p_op[b.ovvv]))
+
+ def final_block(func):
+ return np.array([
+ [prod_sum(self.reference_state, func(qed_t1, i.ph, j.ph))
+ for i in self.excitation_vector]
+ for j in self.excitation_vector])
+
+ block_dict = {}
+ block_dict["couple"] = final_block(couple)
+ block_dict["phot_couple"] = final_block(phot_couple)
+
+ return block_dict
+ else:
+ return ("qed_second_order_ph_ph_couplings are only calculated,"
+ "if reference_state contains 'approx' attribute")
+
@cached_property
@mark_excitation_property()
@timed_member_call(timer="_property_timer")
diff --git a/adcc/LazyMp.py b/adcc/LazyMp.py
index 70e4eb43..e441cc10 100644
--- a/adcc/LazyMp.py
+++ b/adcc/LazyMp.py
@@ -47,6 +47,12 @@ def __init__(self, hf):
self.mospaces = hf.mospaces
self.timer = Timer()
self.has_core_occupied_space = hf.has_core_occupied_space
+ # for qed mp
+ self.get_qed_total_dip = OneParticleOperator(self.mospaces, is_symmetric=True) # noqa: E501
+ self.get_qed_total_dip.oo = hf.get_qed_total_dip(b.oo)
+ self.get_qed_total_dip.ov = hf.get_qed_total_dip(b.ov)
+ self.get_qed_total_dip.vv = hf.get_qed_total_dip(b.vv)
+ self.omega = hf.get_qed_omega()
def __getattr__(self, attr):
# Shortcut some quantities, which are needed most often
@@ -122,7 +128,7 @@ def t2eri(self, space, contraction):
@timed_member_call(timer="timer")
def mp2_diffdm(self):
"""
- Return the MP2 differensce density in the MO basis.
+ Return the MP2 difference density in the MO basis.
"""
hf = self.reference_state
ret = OneParticleOperator(self.mospaces, is_symmetric=True)
@@ -182,6 +188,66 @@ def density(self, level=2):
raise NotImplementedError("Only densities for level 1 and 2"
" are implemented.")
+ @cached_member_function
+ def qed_t1_df(self, space):
+ """
+ qed_t1 amplitude times (df+omega)
+ """
+ total_dip = OneParticleOperator(self.mospaces, is_symmetric=True)
+ if space == b.oo:
+ total_dip.oo = self.get_qed_total_dip.oo
+ return total_dip.oo
+ elif space == b.ov:
+ total_dip.ov = self.get_qed_total_dip.ov
+ return total_dip.ov
+ elif space == b.vv:
+ total_dip.vv = self.get_qed_total_dip.vv
+ return total_dip.vv
+
+ @cached_member_function
+ def qed_t1(self, space):
+ """
+ Return new electronic singly excited amplitude in the first
+ order correction to the wavefunction for qed
+ """
+ if space != b.ov:
+ raise NotImplementedError("qed_t1 term not implemented "
+ f"for space {space}.")
+ omega = ReferenceState.get_qed_omega(self.reference_state)
+ return self.qed_t1_df(b.ov) / (self.df(b.ov) + omega)
+
+ @cached_member_function
+ def qed_t0_df(self, space):
+ """
+ qed_t0 amplitude times df
+ """
+ total_dip = OneParticleOperator(self.mospaces, is_symmetric=True)
+ total_dip.oo = self.get_qed_total_dip.oo
+ total_dip.ov = self.get_qed_total_dip.ov
+ total_dip.vv = self.get_qed_total_dip.vv
+ if space == b.ov:
+ occ_sum = einsum("ka,ki->ia", total_dip.ov, total_dip.oo)
+ virt_sum = einsum("ac,ic->ia", total_dip.vv, total_dip.ov)
+ elif space == b.oo:
+ occ_sum = einsum("ki,kj->ij", total_dip.oo, total_dip.oo)
+ virt_sum = einsum("ic,jc->ij", total_dip.ov, total_dip.ov)
+ elif space == b.vv:
+ occ_sum = einsum("ka,kb->ab", total_dip.ov, total_dip.ov)
+ virt_sum = einsum("ac,bc->ab", total_dip.vv, total_dip.vv)
+ return occ_sum - virt_sum
+
+ @cached_member_function
+ def qed_t0(self, space):
+ """
+ Return second new electronic singly excited amplitude in the first
+ order correction to the wavefunction for qed from the standard
+ HF reference
+ """
+ if space != b.ov:
+ raise NotImplementedError("qed_t0 term not implemented "
+ f"for space {space}.")
+ return self.qed_t0_df(b.ov) / self.df(b.ov)
+
def dipole_moment(self, level=2):
"""
Return the MP dipole moment at the specified level of
@@ -201,12 +267,41 @@ def energy_correction(self, level=2):
if level > 3:
raise NotImplementedError(f"MP({level}) energy correction "
"not implemented.")
+ # For qed_mp1 from non-qed-hf also first order corrections come into play,
+ # which is for now done in mp2 part here
if level < 2:
return 0.0
hf = self.reference_state
is_cvs = self.has_core_occupied_space
+ qed_mp2_correction = 0
if level == 2 and not is_cvs:
terms = [(1.0, hf.oovv, self.t2oo)]
+ if hf.is_qed:
+ total_dip = OneParticleOperator(self.mospaces, is_symmetric=True)
+ omega = ReferenceState.get_qed_omega(hf)
+ total_dip.ov = self.get_qed_total_dip.ov
+ qed_terms = [(omega / 2, total_dip.ov, self.qed_t1(b.ov))]
+ qed_mp2_correction_1 = sum(
+ -pref * lambda_dip.dot(qed_t)
+ for pref, lambda_dip, qed_t in qed_terms
+ )
+ if hf.qed_hf:
+ qed_mp2_correction = qed_mp2_correction_1
+ else:
+ qed_terms_0 = [(1.0, self.qed_t0(b.ov), self.qed_t0_df(b.ov))]
+ qed_mp2_correction_0 = sum(
+ -0.25 * pref * ampl_t0.dot(ampl_t0_df)
+ for pref, ampl_t0, ampl_t0_df in qed_terms_0
+ )
+ # mp1 terms:
+ qed_mp1_additional_terms = [(0.5, total_dip.ov)]
+ qed_mp1_correction = sum(
+ pref * lambda_dip.dot(lambda_dip)
+ for pref, lambda_dip in qed_mp1_additional_terms
+ )
+ qed_mp2_correction = qed_mp2_correction_1 +\
+ qed_mp2_correction_0 +\
+ qed_mp1_correction
elif level == 2 and is_cvs:
terms = [(1.0, hf.oovv, self.t2oo),
(2.0, hf.ocvv, self.t2oc),
@@ -218,7 +313,7 @@ def energy_correction(self, level=2):
return sum(
-0.25 * pref * eri.dot(t2)
for pref, eri, t2 in terms
- )
+ ) + qed_mp2_correction
def energy(self, level=2):
"""
diff --git a/adcc/ReferenceState.py b/adcc/ReferenceState.py
index a3274c61..befddf71 100644
--- a/adcc/ReferenceState.py
+++ b/adcc/ReferenceState.py
@@ -28,6 +28,7 @@
from .backends import import_scf_results
from .OperatorIntegrals import OperatorIntegrals
from .OneParticleOperator import OneParticleOperator, product_trace
+from .misc import cached_member_function
import libadcc
@@ -37,7 +38,8 @@
class ReferenceState(libadcc.ReferenceState):
def __init__(self, hfdata, core_orbitals=None, frozen_core=None,
frozen_virtual=None, symmetry_check_on_import=False,
- import_all_below_n_orbs=10):
+ import_all_below_n_orbs=10, is_qed=False, coupl=None,
+ freq=None, qed_hf=True, qed_approx=False):
"""Construct a ReferenceState holding information about the employed
SCF reference.
@@ -138,6 +140,13 @@ def __init__(self, hfdata, core_orbitals=None, frozen_core=None,
which would place the 2nd and 3rd alpha and the 1st and second
beta orbital into the core space.
"""
+ self.is_qed = is_qed
+ self.coupling = coupl
+ self.frequency = np.real(freq)
+ self.freq_with_loss = freq
+ self.qed_hf = qed_hf
+ self.approx = qed_approx
+
if not isinstance(hfdata, libadcc.HartreeFockSolution_i):
hfdata = import_scf_results(hfdata)
@@ -165,9 +174,98 @@ def __getattr__(self, attr):
if attr.startswith("f"):
return self.fock(b.__getattr__(attr[1:]))
+ elif attr.startswith("get_qed_total_dip"):
+ return self.get_qed_total_dip(b.__getattr__(attr))
+ elif attr.startswith("get_qed_omega"):
+ return self.get_qed_omega
else:
return self.eri(b.__getattr__(attr))
+ @cached_member_function
+ def get_qed_total_dip(self, block):
+ """
+ Return qed coupling strength times dipole operator
+ """
+ # TODO: Here we always multiply with sqrt(2 * omega), since this
+ # eases up the Hamiltonian and is required if you provide a hilbert
+ # package QED-HF input. This can differ between QED-HF implementations,
+ # e.g. the psi4numpy QED-RHF helper does not do that. Therefore, this
+ # factor needs to be adjusted depending on the input, but since the
+ # hilbert package is currently the best in terms of performance, at
+ # least to my knowledge, the factor should be included here.
+ if self.is_qed:
+ dips = self.operators.electric_dipole
+ couplings = self.coupling
+ freqs = self.frequency
+ total_dip = OneParticleOperator(self.mospaces, is_symmetric=True)
+ for coupling, dip in zip(couplings, dips):
+ total_dip += coupling * dip
+ total_dip *= np.sqrt(2 * np.linalg.norm(freqs))
+ total_dip.evaluate()
+ return total_dip[block]
+
+ @cached_member_function
+ def get_qed_omega(self):
+ """
+ Return the cavity frequency
+ """
+ if self.is_qed:
+ freqs = self.frequency
+ return np.linalg.norm(freqs)
+
+ @cached_member_function
+ def qed_D_object(self, block):
+ """
+ Return the object, which is added to the ERIs in a PT QED calculation
+ """
+ if self.is_qed:
+ from . import block as b
+ from .functions import einsum
+ total_dip = OneParticleOperator(self.mospaces, is_symmetric=True)
+ total_dip.oo = ReferenceState.get_qed_total_dip(self, b.oo)
+ total_dip.ov = ReferenceState.get_qed_total_dip(self, b.ov)
+ total_dip.vv = ReferenceState.get_qed_total_dip(self, b.vv)
+ # We have to define all the blocks from the
+ # D_{pqrs} = d_{pr} d_{qs} - d_{ps} d_{qr} object, which has the
+ # same symmetry properties as the ERI object
+ # Actually in b.ovov: second term: ib,ja would be ib,aj ,
+ # but d_{ia} = d_{ai}
+ ds = {
+ b.oooo: einsum('ik,jl->ijkl', total_dip.oo, total_dip.oo)
+ - einsum('il,jk->ijkl', total_dip.oo, total_dip.oo),
+ b.ooov: einsum('ik,ja->ijka', total_dip.oo, total_dip.ov)
+ - einsum('ia,jk->ijka', total_dip.ov, total_dip.oo),
+ b.oovv: einsum('ia,jb->ijab', total_dip.ov, total_dip.ov)
+ - einsum('ib,ja->ijab', total_dip.ov, total_dip.ov),
+ b.ovvv: einsum('ib,ac->iabc', total_dip.ov, total_dip.vv)
+ - einsum('ic,ab->iabc', total_dip.ov, total_dip.vv),
+ b.ovov: einsum('ij,ab->iajb', total_dip.oo, total_dip.vv)
+ - einsum('ib,ja->iajb', total_dip.ov, total_dip.ov),
+ b.vvvv: einsum('ac,bd->abcd', total_dip.vv, total_dip.vv)
+ - einsum('ad,bc->abcd', total_dip.vv, total_dip.vv),
+ }
+ return ds[block]
+
+ def eri(self, block):
+ if self.is_qed:
+ from . import block as b
+ from .functions import einsum
+ # Since there is no TwoParticleOperator object,
+ # we initialize it like this
+ ds_init = OneParticleOperator(self.mospaces, is_symmetric=True)
+ ds = {
+ b.oooo: einsum('ik,jl->ijkl', ds_init.oo, ds_init.oo),
+ b.ooov: einsum('ik,ja->ijka', ds_init.oo, ds_init.ov),
+ b.oovv: einsum('ia,jb->ijab', ds_init.ov, ds_init.ov),
+ b.ovvv: einsum('ib,ac->iabc', ds_init.ov, ds_init.vv),
+ b.ovov: einsum('ij,ab->iajb', ds_init.oo, ds_init.vv),
+ b.vvvv: einsum('ac,bd->abcd', ds_init.vv, ds_init.vv),
+ }
+ ds[block] = ReferenceState.qed_D_object(self, block)
+ return super().eri(block) + ds[block]
+ else:
+ return super().eri(block)
+
@property
def mospaces(self):
return self._mospaces
diff --git a/adcc/adc_pp/matrix.py b/adcc/adc_pp/matrix.py
index 80c489ca..d86500ef 100644
--- a/adcc/adc_pp/matrix.py
+++ b/adcc/adc_pp/matrix.py
@@ -27,6 +27,8 @@
from adcc.functions import direct_sum, einsum, zeros_like
from adcc.Intermediates import Intermediates, register_as_intermediate
from adcc.AmplitudeVector import AmplitudeVector
+from adcc.ReferenceState import ReferenceState
+from libadcc import set_lt_scalar
__all__ = ["block"]
@@ -85,33 +87,133 @@ def block(ground_state, spaces, order, variant=None, intermediates=None):
return globals()[fn](reference_state, ground_state, intermediates)
+# For QED-ADC (up to double photon dispersion) we build the matrix as follows:
+# elec phot_couple phot_couple_edge
+# elec_couple phot phot_couple_inner
+# elec_couple_edge elec_couple_inner phot2
+# where each block is a "standard" ADC matrix itself, including the groundstate
+# and the groundstate couplings. However, the gs_ph and gs_gs blocks are merged
+# into the ph_ph and gs_gs blocks, respectively, since we calculate matrix vector
+# products anyway. Note, that the purely electric groundstate never appears, since
+# it is always zero.
+
+
+#
+# 0th order gs blocks
+#
+
+def block_ph_gs_0(hf, mp, intermediates):
+ return AdcBlock(lambda ampl: 0, 0)
+
+
+def block_ph_gs_0_phot(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(gs1=omega * ampl.gs1)
+ return AdcBlock(apply, AmplitudeVector(gs1=set_lt_scalar(omega)))
+
+
+def block_ph_gs_0_phot2(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(gs2=2 * omega * ampl.gs2)
+ return AdcBlock(apply, AmplitudeVector(gs2=set_lt_scalar(2 * omega)))
+
+
+block_ph_gs_0_couple = block_ph_gs_0_phot_couple =\
+ block_ph_gs_0_couple_edge = block_ph_gs_0_phot_couple_edge =\
+ block_ph_gs_0_couple_inner = block_ph_gs_0_phot_couple_inner =\
+ block_ph_gs_0
+
+block_cvs_ph_gs_0 = block_cvs_ph_gs_1 = block_cvs_ph_gs_2 = block_ph_gs_0
+
+
#
# 0th order main
#
+
def block_ph_ph_0(hf, mp, intermediates):
fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
- diagonal = AmplitudeVector(ph=direct_sum("a-i->ia", hf.fvv.diagonal(),
- fCC.diagonal()))
-
- def apply(ampl):
- return AmplitudeVector(ph=(
- + einsum("ib,ab->ia", ampl.ph, hf.fvv)
- - einsum("IJ,Ja->Ia", fCC, ampl.ph)
- ))
+ if hf.is_qed:
+ diagonal = AmplitudeVector(ph=direct_sum("a-i->ia", hf.fvv.diagonal(),
+ fCC.diagonal()))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ + einsum("ib,ab->ia", ampl.ph, hf.fvv)
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph)
+ ))
+ else:
+ diagonal = AmplitudeVector(ph=direct_sum("a-i->ia", hf.fvv.diagonal(),
+ fCC.diagonal()))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ + einsum("ib,ab->ia", ampl.ph, hf.fvv)
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph)
+ ))
return AdcBlock(apply, diagonal)
block_cvs_ph_ph_0 = block_ph_ph_0
-def diagonal_pphh_pphh_0(hf):
+def block_ph_ph_0_couple(hf, mp, intermediates):
+ return AdcBlock(lambda ampl: 0, 0)
+
+
+block_ph_ph_0_phot_couple = block_ph_ph_0_phot_couple_edge =\
+ block_ph_ph_0_phot_couple_inner = block_ph_ph_0_couple_edge =\
+ block_ph_ph_0_couple_inner = block_ph_ph_0_couple
+
+
+def block_ph_ph_0_phot(hf, mp, intermediates):
+ fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
+ if hf.is_qed:
+ diagonal = AmplitudeVector(ph1=direct_sum("a-i->ia", hf.fvv.diagonal(),
+ fCC.diagonal()))
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ + einsum("ib,ab->ia", ampl.ph1, hf.fvv)
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph1)
+ ))
+ else:
+ raise NotImplementedError("block_ph_ph_0_phot is requested, "
+ "but ReferenceState has no coupling attribute")
+ return AdcBlock(apply, diagonal)
+
+
+def block_ph_ph_0_phot2(hf, mp, intermediates):
+ fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
+ diagonal = AmplitudeVector(ph2=direct_sum("a-i->ia", hf.fvv.diagonal(),
+ fCC.diagonal()))
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=(
+ + einsum("ib,ab->ia", ampl.ph2, hf.fvv)
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph2)
+ ))
+ return AdcBlock(apply, diagonal)
+
+
+def diagonal_pphh_pphh_0(hf, n_omega=None):
# Note: adcman similarly does not symmetrise the occupied indices
# (for both CVS and general ADC)
fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
- res = direct_sum("-i-J+a+b->iJab",
- hf.foo.diagonal(), fCC.diagonal(),
- hf.fvv.diagonal(), hf.fvv.diagonal())
- return AmplitudeVector(pphh=res.symmetrise(2, 3))
+ if n_omega is not None:
+ omega = float(ReferenceState.get_qed_omega(hf))
+ qed = n_omega * omega
+ res = direct_sum("-i-J+a+b->iJab",
+ hf.foo.diagonal() + qed, fCC.diagonal() + qed,
+ hf.fvv.diagonal() + qed, hf.fvv.diagonal() + qed)
+ else:
+ res = direct_sum("-i-J+a+b->iJab",
+ hf.foo.diagonal(), fCC.diagonal(),
+ hf.fvv.diagonal(), hf.fvv.diagonal())
+ return res.symmetrise(2, 3)
def block_pphh_pphh_0(hf, mp, intermediates):
@@ -120,7 +222,28 @@ def apply(ampl):
+ 2 * einsum("ijac,bc->ijab", ampl.pphh, hf.fvv).antisymmetrise(2, 3)
- 2 * einsum("ik,kjab->ijab", hf.foo, ampl.pphh).antisymmetrise(0, 1)
))
- return AdcBlock(apply, diagonal_pphh_pphh_0(hf))
+ return AdcBlock(apply, AmplitudeVector(pphh=diagonal_pphh_pphh_0(hf)))
+
+
+def block_pphh_pphh_0_couple(hf, mp, intermediates):
+ return AdcBlock(lambda ampl: 0, 0)
+
+
+block_pphh_pphh_0_phot_couple = block_pphh_pphh_0_phot_couple_edge =\
+ block_pphh_pphh_0_phot_couple_inner = block_pphh_pphh_0_couple_edge =\
+ block_pphh_pphh_0_couple_inner = block_pphh_pphh_0_couple
+
+
+def block_pphh_pphh_0_phot(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(pphh1=(
+ + 2 * einsum("ijac,bc->ijab", ampl.pphh1, hf.fvv).antisymmetrise(2, 3)
+ - 2 * einsum("ik,kjab->ijab", hf.foo, ampl.pphh1).antisymmetrise(0, 1)
+ + omega * ampl.pphh1
+ ))
+ return AdcBlock(apply, AmplitudeVector(pphh1=diagonal_pphh_pphh_0(hf, 1)))
def block_cvs_pphh_pphh_0(hf, mp, intermediates):
@@ -130,12 +253,25 @@ def apply(ampl):
- einsum("ik,kJab->iJab", hf.foo, ampl.pphh)
- einsum("JK,iKab->iJab", hf.fcc, ampl.pphh)
))
- return AdcBlock(apply, diagonal_pphh_pphh_0(hf))
+ return AdcBlock(apply, AmplitudeVector(pphh=diagonal_pphh_pphh_0(hf)))
+
+
+def block_pphh_pphh_0_phot2(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(pphh2=(
+ + 2 * einsum("ijac,bc->ijab", ampl.pphh2, hf.fvv).antisymmetrise(2, 3)
+ - 2 * einsum("ik,kjab->ijab", hf.foo, ampl.pphh2).antisymmetrise(0, 1)
+ + 2 * omega * ampl.pphh2
+ ))
+ return AdcBlock(apply, AmplitudeVector(pphh2=diagonal_pphh_pphh_0(hf, 2)))
#
# 0th order coupling
#
+
def block_ph_pphh_0(hf, mp, intermediates):
return AdcBlock(lambda ampl: 0, 0)
@@ -147,30 +283,249 @@ def block_pphh_ph_0(hf, mp, intermediates):
block_cvs_ph_pphh_0 = block_ph_pphh_0
block_cvs_pphh_ph_0 = block_pphh_ph_0
+block_pphh_ph_0_couple = block_pphh_ph_0_phot_couple = block_pphh_ph_0_phot =\
+ block_pphh_ph_0_couple_edge = block_pphh_ph_0_couple_inner =\
+ block_pphh_ph_0_phot_couple_edge = block_pphh_ph_0_phot_couple_inner =\
+ block_pphh_ph_0_phot2 = block_pphh_ph_0
+
+block_ph_pphh_0_couple = block_ph_pphh_0_phot_couple = block_ph_pphh_0_phot =\
+ block_ph_pphh_0_couple_edge = block_ph_pphh_0_couple_inner =\
+ block_ph_pphh_0_phot_couple_edge = block_ph_pphh_0_phot_couple_inner =\
+ block_ph_pphh_0_phot2 = block_ph_pphh_0
+
+
+#
+# 1st order gs blocks
+#
+
+block_ph_gs_1_phot = block_ph_gs_0_phot
+block_ph_gs_1_phot2 = block_ph_gs_0_phot2
+
+
+block_ph_gs_1_phot_couple = block_ph_gs_1_phot_couple_edge =\
+ block_ph_gs_1_couple_edge = block_ph_gs_1_phot_couple_inner =\
+ block_ph_gs_1 = block_ph_gs_0
+
+
+def block_ph_gs_1_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(gs1=(
+ (-1) * sqrt(0.5 * omega) * mp.qed_t1_df(b.ov).dot(ampl.ph)
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_gs_1_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(gs2=(
+ (-1) * sqrt(omega) * mp.qed_t1_df(b.ov).dot(ampl.ph1)
+ ))
+ return AdcBlock(apply, 0)
+
#
# 1st order main
#
+
def block_ph_ph_1(hf, mp, intermediates):
fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
CvCv = hf.cvcv if hf.has_core_occupied_space else hf.ovov
- diagonal = AmplitudeVector(ph=(
- + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # order 0
- - einsum("IaIa->Ia", CvCv) # order 1
- ))
+ if hf.is_qed and not hf.qed_hf:
+ diagonal = AmplitudeVector(ph=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # order 0
+ - einsum("IaIa->Ia", CvCv) # order 1
+ + (1 / 2) * direct_sum("i-a->ia", einsum("ii->i", mp.qed_t0_df(b.oo)),
+ einsum("aa->a", mp.qed_t0_df(b.vv)))
+ ))
- def apply(ampl):
- return AmplitudeVector(ph=( # PT order
- + einsum("ib,ab->ia", ampl.ph, hf.fvv) # 0
- - einsum("IJ,Ja->Ia", fCC, ampl.ph) # 0
- - einsum("JaIb,Jb->Ia", CvCv, ampl.ph) # 1
+ def apply(ampl):
+ return AmplitudeVector(ph=( # PT order
+ + einsum("ib,ab->ia", ampl.ph, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph) # 1
+ + (1 / 2) * einsum("ij,ja->ia", mp.qed_t0_df(b.oo), ampl.ph)
+ - (1 / 2) * einsum("ib,ab->ia", ampl.ph, mp.qed_t0_df(b.vv))
+ ))
+ elif hf.is_qed and hf.qed_hf:
+ diagonal = AmplitudeVector(ph=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # 0
+ - einsum("IaIa->Ia", CvCv) # 1
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=( # PT order
+ + einsum("ib,ab->ia", ampl.ph, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph) # 1
+ ))
+
+ else:
+ diagonal = AmplitudeVector(ph=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # order 0
+ - einsum("IaIa->Ia", CvCv) # order 1
))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=( # PT order
+ + einsum("ib,ab->ia", ampl.ph, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph) # 1
+ ))
return AdcBlock(apply, diagonal)
block_cvs_ph_ph_1 = block_ph_ph_1
+def block_ph_ph_1_phot(hf, mp, intermediates):
+ fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
+ CvCv = hf.cvcv if hf.has_core_occupied_space else hf.ovov
+ if hf.is_qed and not hf.qed_hf:
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ diagonal = AmplitudeVector(ph1=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # order 0
+ - einsum("IaIa->Ia", CvCv) # order 1
+ + (1 / 2) * direct_sum("i-a->ia", einsum("ii->i", mp.qed_t0_df(b.oo)),
+ einsum("aa->a", mp.qed_t0_df(b.vv)))
+ + intermediates.delta_ia_omega
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=( # PT order
+ + einsum("ib,ab->ia", ampl.ph1, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph1) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph1) # 1
+ + (1 / 2) * einsum("ij,ja->ia", mp.qed_t0_df(b.oo), ampl.ph1)
+ - (1 / 2) * einsum("ib,ab->ia", ampl.ph1, mp.qed_t0_df(b.vv))
+ + omega * ampl.ph1
+ ))
+ elif hf.is_qed and hf.qed_hf:
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ diagonal = AmplitudeVector(ph1=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # 0
+ - einsum("IaIa->Ia", CvCv) # 1
+ + intermediates.delta_ia_omega
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=( # PT order
+ + einsum("ib,ab->ia", ampl.ph1, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph1) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph1) # 1
+ + omega * ampl.ph1
+ ))
+ else:
+ raise NotImplementedError("block_ph_ph_1_phot is requested, "
+ "but ReferenceState has no coupling attribute")
+ return AdcBlock(apply, diagonal)
+
+
+def block_ph_ph_1_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ if hf.is_qed:
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ sqrt(omega / 2) * (
+ - einsum("ib,ab->ia", ampl.ph, mp.qed_t1_df(b.vv))
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph))
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_1_phot_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ if hf.is_qed:
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ sqrt(omega / 2) * (
+ - einsum("ib,ab->ia", ampl.ph1, mp.qed_t1_df(b.vv))
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph1)
+ - mp.qed_t1_df(b.ov) * ampl.gs1.to_ndarray())
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_1_couple_edge(hf, mp, intermediates):
+ return AdcBlock(lambda ampl: 0, 0)
+
+
+block_ph_ph_1_phot_couple_edge = block_ph_ph_1_couple_edge
+
+
+def block_ph_ph_1_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ if hf.is_qed:
+ def apply(ampl):
+ return AmplitudeVector(ph2=(
+ sqrt(omega) * (- einsum("ib,ab->ia", ampl.ph1, mp.qed_t1_df(b.vv))
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph1))
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_1_phot_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ if hf.is_qed:
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ sqrt(omega) * (- einsum("ib,ab->ia", ampl.ph2, mp.qed_t1_df(b.vv))
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph2)
+ - mp.qed_t1_df(b.ov) * ampl.gs2.to_ndarray())
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_1_phot2(hf, mp, intermediates):
+ fCC = hf.fcc if hf.has_core_occupied_space else hf.foo
+ CvCv = hf.cvcv if hf.has_core_occupied_space else hf.ovov
+ if hf.is_qed and not hf.qed_hf:
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ diagonal = AmplitudeVector(ph2=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # order 0
+ - einsum("IaIa->Ia", CvCv) # order 1
+ + (1 / 2) * direct_sum("i-a->ia", einsum("ii->i", mp.qed_t0_df(b.oo)),
+ einsum("aa->a", mp.qed_t0_df(b.vv)))
+ + intermediates.delta_ia_omega * 2
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=( # PT order
+ + einsum("ib,ab->ia", ampl.ph2, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph2) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph2) # 1
+ + (1 / 2) * einsum("ij,ja->ia", mp.qed_t0_df(b.oo), ampl.ph2)
+ - (1 / 2) * einsum("ib,ab->ia", ampl.ph2, mp.qed_t0_df(b.vv))
+ + 2 * omega * ampl.ph2
+ ))
+ elif hf.is_qed and hf.qed_hf:
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ diagonal = AmplitudeVector(ph2=(
+ + direct_sum("a-i->ia", hf.fvv.diagonal(), fCC.diagonal()) # 0
+ - einsum("IaIa->Ia", CvCv) # 1
+ + intermediates.delta_ia_omega * 2
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=( # PT order
+ + einsum("ib,ab->ia", ampl.ph2, hf.fvv) # 0
+ - einsum("IJ,Ja->Ia", fCC, ampl.ph2) # 0
+ - einsum("JaIb,Jb->Ia", CvCv, ampl.ph2) # 1
+ + 2 * omega * ampl.ph2
+ ))
+ else:
+ raise NotImplementedError("block_ph_ph_1_phot2 is requested, "
+ "but ReferenceState has no coupling attribute")
+ return AdcBlock(apply, diagonal)
+
+
def diagonal_pphh_pphh_1(hf):
# Fock matrix and ovov diagonal term (sometimes called "intermediate diagonal")
dinterm_ov = (direct_sum("a-i->ia", hf.fvv.diagonal(), hf.foo.diagonal())
@@ -227,6 +582,7 @@ def apply(ampl):
#
# 1st order coupling
#
+
def block_ph_pphh_1(hf, mp, intermediates):
def apply(ampl):
return AmplitudeVector(ph=(
@@ -236,6 +592,24 @@ def apply(ampl):
return AdcBlock(apply, 0)
+def block_ph_pphh_1_phot(hf, mp, intermediates):
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ + einsum("jkib,jkab->ia", hf.ooov, ampl.pphh1)
+ + einsum("ijbc,jabc->ia", ampl.pphh1, hf.ovvv)
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_pphh_1_phot2(hf, mp, intermediates):
+ def apply(ampl):
+ return AmplitudeVector(ph2=(
+ + einsum("jkib,jkab->ia", hf.ooov, ampl.pphh2)
+ + einsum("ijbc,jabc->ia", ampl.pphh2, hf.ovvv)
+ ))
+ return AdcBlock(apply, 0)
+
+
def block_cvs_ph_pphh_1(hf, mp, intermediates):
def apply(ampl):
return AmplitudeVector(ph=(
@@ -254,6 +628,24 @@ def apply(ampl):
return AdcBlock(apply, 0)
+def block_pphh_ph_1_phot(hf, mp, intermediates):
+ def apply(ampl):
+ return AmplitudeVector(pphh1=(
+ + einsum("ic,jcab->ijab", ampl.ph1, hf.ovvv).antisymmetrise(0, 1)
+ - einsum("ijka,kb->ijab", hf.ooov, ampl.ph1).antisymmetrise(2, 3)
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_pphh_ph_1_phot2(hf, mp, intermediates):
+ def apply(ampl):
+ return AmplitudeVector(pphh2=(
+ + einsum("ic,jcab->ijab", ampl.ph2, hf.ovvv).antisymmetrise(0, 1)
+ - einsum("ijka,kb->ijab", hf.ooov, ampl.ph2).antisymmetrise(2, 3)
+ ))
+ return AdcBlock(apply, 0)
+
+
def block_cvs_pphh_ph_1(hf, mp, intermediates):
def apply(ampl):
return AmplitudeVector(pphh=(
@@ -264,31 +656,151 @@ def apply(ampl):
return AdcBlock(apply, 0)
+def block_ph_pphh_1_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ 2 * sqrt(omega / 2) * einsum("kc,ikac->ia", mp.qed_t1_df(b.ov), ampl.pphh) # noqa: E501
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_pphh_1_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=(
+ 2 * sqrt(omega) * einsum("kc,ikac->ia", mp.qed_t1_df(b.ov), ampl.pphh1)
+ ))
+ return AdcBlock(apply, 0)
+
+
+block_pphh_ph_1_couple = block_pphh_ph_1_couple_inner =\
+ block_pphh_ph_1_phot_couple_edge = block_pphh_ph_1_couple_edge =\
+ block_pphh_ph_0_couple
+
+block_ph_pphh_1_phot_couple = block_ph_pphh_1_phot_couple_inner =\
+ block_ph_pphh_1_phot_couple_edge = block_ph_pphh_1_couple_edge =\
+ block_ph_pphh_0_phot_couple
+
+
+def block_pphh_ph_1_phot_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(pphh=(
+ 2 * sqrt(omega / 2) * einsum(
+ "jb,ia->ijab", mp.qed_t1_df(b.ov),
+ ampl.ph1).antisymmetrise(0, 1).antisymmetrise(2, 3)
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_pphh_ph_1_phot_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(pphh1=(
+ 2 * sqrt(omega) * einsum(
+ "jb,ia->ijab", mp.qed_t1_df(b.ov),
+ ampl.ph2).antisymmetrise(0, 1).antisymmetrise(2, 3)
+ ))
+ return AdcBlock(apply, 0)
+
+
+#
+# 2nd order gs blocks
+#
+
+block_ph_gs_2_phot_couple = block_ph_gs_2_phot_couple_edge =\
+ block_ph_gs_2_couple_edge = block_ph_gs_2_phot_couple_inner =\
+ block_ph_gs_2 = block_ph_gs_0
+
+
+def block_ph_gs_2_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(gs1=(sqrt(omega / 2) * einsum(
+ "jkbc,kc->jb", mp.t2oo, mp.qed_t1_df(b.ov)).dot(ampl.ph)
+ - sqrt(0.5 * omega) * mp.qed_t1_df(b.ov).dot(ampl.ph))) # 1. order
+ return AdcBlock(apply, 0)
+
+
+block_ph_gs_2_phot = block_ph_gs_0_phot
+block_ph_gs_2_phot2 = block_ph_gs_0_phot2
+
+
+def block_ph_gs_2_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(gs2=(sqrt(omega) * einsum(
+ "jkbc,kc->jb", mp.t2oo, mp.qed_t1_df(b.ov)).dot(ampl.ph1)
+ - sqrt(omega) * mp.qed_t1_df(b.ov).dot(ampl.ph1))) # 1. order
+ return AdcBlock(apply, 0)
+
+
#
# 2nd order main
#
+
def block_ph_ph_2(hf, mp, intermediates):
i1 = intermediates.adc2_i1
i2 = intermediates.adc2_i2
- diagonal = AmplitudeVector(ph=(
- + direct_sum("a-i->ia", i1.diagonal(), i2.diagonal())
- - einsum("IaIa->Ia", hf.ovov)
- - einsum("ikac,ikac->ia", mp.t2oo, hf.oovv)
- ))
- # Not used anywhere else, so kept as an anonymous intermediate
- term_t2_eri = (
- + einsum("ijab,jkbc->ikac", mp.t2oo, hf.oovv)
- + einsum("ijab,jkbc->ikac", hf.oovv, mp.t2oo)
- ).evaluate()
-
- def apply(ampl):
- return AmplitudeVector(ph=(
- + einsum("ib,ab->ia", ampl.ph, i1)
- - einsum("ij,ja->ia", i2, ampl.ph)
- - einsum("jaib,jb->ia", hf.ovov, ampl.ph) # 1
- - 0.5 * einsum("ikac,kc->ia", term_t2_eri, ampl.ph) # 2
+ term_t2_eri = intermediates.term_t2_eri
+
+ if hf.is_qed and not hf.approx:
+ omega = float(ReferenceState.get_qed_omega(hf))
+ if hf.qed_hf:
+ qed_i1 = intermediates.adc2_qed_i1
+ qed_i2 = intermediates.adc2_qed_i2
+ diagonal = AmplitudeVector(ph=(
+ + direct_sum("a-i->ia", i1.diagonal(), i2.diagonal())
+ - einsum("IaIa->Ia", hf.ovov)
+ - einsum("ikac,ikac->ia", mp.t2oo, hf.oovv)
+ + (-omega / 2) * (
+ - direct_sum("a+i->ia", qed_i1.diagonal(), qed_i2.diagonal())
+ + (1 / 2) * 2 * einsum(
+ "ia,ia->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)
+ )
+ )
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ + einsum("ib,ab->ia", ampl.ph, i1)
+ - einsum("ij,ja->ia", i2, ampl.ph)
+ - einsum("jaib,jb->ia", hf.ovov, ampl.ph) # 1
+ - 0.5 * einsum("ikac,kc->ia", term_t2_eri, ampl.ph) # 2
+ + (-omega / 2) * (
+ - einsum("ib,ab->ia", ampl.ph, qed_i1)
+ - einsum("ij,ja->ia", qed_i2, ampl.ph)
+ + (1 / 2) * (
+ mp.qed_t1(b.ov) * mp.qed_t1_df(b.ov).dot(ampl.ph)
+ + mp.qed_t1_df(b.ov) * mp.qed_t1(b.ov).dot(ampl.ph)
+ )
+ )
+ ))
+ else:
+ raise NotImplementedError("QED-ADC(2) from non-QED-HF reference"
+ "is not implemented")
+ else:
+ diagonal = AmplitudeVector(ph=(
+ + direct_sum("a-i->ia", i1.diagonal(), i2.diagonal())
+ - einsum("IaIa->Ia", hf.ovov)
+ - einsum("ikac,ikac->ia", mp.t2oo, hf.oovv)
))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ + einsum("ib,ab->ia", ampl.ph, i1)
+ - einsum("ij,ja->ia", i2, ampl.ph)
+ - einsum("jaib,jb->ia", hf.ovov, ampl.ph) # 1
+ - 0.5 * einsum("ikac,kc->ia", term_t2_eri, ampl.ph) # 2
+ ))
return AdcBlock(apply, diagonal)
@@ -308,9 +820,192 @@ def apply(ampl):
return AdcBlock(apply, diagonal)
+def block_ph_ph_2_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ qed_i1 = intermediates.adc2_qed_couple_i1
+ qed_i2 = intermediates.adc2_qed_couple_i2
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ + einsum("ib,ab->ia", ampl.ph, qed_i1)
+ + einsum("ij,ja->ia", qed_i2, ampl.ph)
+ + sqrt(omega / 2) * (
+ + einsum("ka,jkib,jb->ia", mp.qed_t1(b.ov), hf.ooov, ampl.ph)
+ + einsum("ic,jabc,jb->ia", mp.qed_t1(b.ov), hf.ovvv, ampl.ph))
+ + sqrt(omega / 2) * (
+ - einsum("ib,ab->ia", ampl.ph, mp.qed_t1_df(b.vv)) # 1. order
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph)) # 1. order
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_2_phot_couple(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ gs_part = intermediates.adc2_qed_ph_ph_2_phot_couple_gs_part
+ qed_i1 = intermediates.adc2_qed_phot_couple_i1
+ qed_i2 = intermediates.adc2_qed_phot_couple_i2
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(
+ + einsum("ib,ab->ia", ampl.ph1, qed_i1)
+ + einsum("ij,ja->ia", qed_i2, ampl.ph1)
+ + sqrt(omega / 2) * (
+ + einsum("kb,ikja,jb->ia", mp.qed_t1(b.ov), hf.ooov, ampl.ph1)
+ + einsum("jc,ibac,jb->ia", mp.qed_t1(b.ov), hf.ovvv, ampl.ph1))
+ + gs_part * ampl.gs1.to_ndarray()
+ + sqrt(omega / 2) * (
+ - einsum("ib,ab->ia", ampl.ph1, mp.qed_t1_df(b.vv)) # 1. order
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph1) # 1. order
+ - mp.qed_t1_df(b.ov) * ampl.gs1.to_ndarray()) # 1. order
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_2_phot(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ i1 = intermediates.adc2_i1
+ i2 = intermediates.adc2_i2
+
+ term_t2_eri = intermediates.term_t2_eri
+
+ qed_i1 = intermediates.adc2_qed_i1
+ qed_i2 = intermediates.adc2_qed_i2
+
+ diagonal = AmplitudeVector(ph1=(
+ + direct_sum("a-i->ia", i1.diagonal(), i2.diagonal())
+ - einsum("IaIa->Ia", hf.ovov)
+ - einsum("ikac,ikac->ia", mp.t2oo, hf.oovv)
+ + (-omega / 2) * 2 * (
+ - direct_sum("a+i->ia", qed_i1.diagonal(), qed_i2.diagonal())
+ + einsum("ia,ia->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)))
+ + intermediates.delta_ia_omega # 0. order
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ + einsum("ib,ab->ia", ampl.ph1, i1)
+ - einsum("ij,ja->ia", i2, ampl.ph1)
+ - einsum("jaib,jb->ia", hf.ovov, ampl.ph1) # 1
+ - 0.5 * einsum("ikac,kc->ia", term_t2_eri, ampl.ph1) # 2
+ + (-omega / 2) * 2 * (
+ - einsum("ib,ab->ia", ampl.ph1, qed_i1)
+ - einsum("ij,ja->ia", qed_i2, ampl.ph1)
+ + 0.5 * (mp.qed_t1(b.ov) * mp.qed_t1_df(b.ov).dot(ampl.ph1)
+ + mp.qed_t1_df(b.ov) * mp.qed_t1(b.ov).dot(ampl.ph1)))
+ + omega * ampl.ph1 # 0. order
+ ))
+
+ if not hasattr(hf, "coupling"):
+ raise NotImplementedError("block_ph_ph_2_phot is requested, "
+ "but ReferenceState has no coupling attribute")
+ return AdcBlock(apply, diagonal)
+
+
+def block_ph_ph_2_phot2(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ i1 = intermediates.adc2_i1
+ i2 = intermediates.adc2_i2
+
+ term_t2_eri = intermediates.term_t2_eri
+
+ qed_i1 = intermediates.adc2_qed_i1
+ qed_i2 = intermediates.adc2_qed_i2
+
+ diagonal = AmplitudeVector(ph2=(
+ + direct_sum("a-i->ia", i1.diagonal(), i2.diagonal())
+ - einsum("IaIa->Ia", hf.ovov)
+ - einsum("ikac,ikac->ia", mp.t2oo, hf.oovv)
+ + (-omega / 2) * 3 * (
+ - direct_sum("a+i->ia", qed_i1.diagonal(), qed_i2.diagonal())
+ + einsum("ia,ia->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)))
+ + intermediates.delta_ia_omega * 2 # 0. order
+ ))
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=(
+ + einsum("ib,ab->ia", ampl.ph2, i1)
+ - einsum("ij,ja->ia", i2, ampl.ph2)
+ - einsum("jaib,jb->ia", hf.ovov, ampl.ph2) # 1
+ - 0.5 * einsum("ikac,kc->ia", term_t2_eri, ampl.ph2) # 2
+ + (-omega / 2) * 3 * (
+ - einsum("ib,ab->ia", ampl.ph2, qed_i1)
+ - einsum("ij,ja->ia", qed_i2, ampl.ph2)
+ + 0.5 * (mp.qed_t1(b.ov) * mp.qed_t1_df(b.ov).dot(ampl.ph2)
+ + mp.qed_t1_df(b.ov) * mp.qed_t1(b.ov).dot(ampl.ph2)))
+ + 2 * omega * ampl.ph2 # 0. order
+ ))
+ return AdcBlock(apply, diagonal)
+
+
+def block_ph_ph_2_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ qed_i1 = intermediates.adc2_qed_couple_i1
+ qed_i2 = intermediates.adc2_qed_couple_i2
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=(
+ + sqrt(2) * einsum("ib,ab->ia", ampl.ph1, qed_i1)
+ + sqrt(2) * einsum("ij,ja->ia", qed_i2, ampl.ph1)
+ + sqrt(omega) * (
+ + einsum("ka,jkib,jb->ia", mp.qed_t1(b.ov), hf.ooov, ampl.ph1)
+ + einsum("ic,jabc,jb->ia", mp.qed_t1(b.ov), hf.ovvv, ampl.ph1))
+ + sqrt(omega) * (
+ - einsum("ib,ab->ia", ampl.ph1, mp.qed_t1_df(b.vv)) # 1. order
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph1)) # 1. order
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_2_phot_couple_inner(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ gs_part = intermediates.adc2_qed_ph_ph_2_phot_couple_inner_gs_part
+ qed_i1 = intermediates.adc2_qed_phot_couple_i1
+ qed_i2 = intermediates.adc2_qed_phot_couple_i2
+
+ def apply(ampl):
+ return AmplitudeVector(ph1=(
+ + sqrt(2) * einsum("ib,ab->ia", ampl.ph2, qed_i1)
+ + sqrt(2) * einsum("ij,ja->ia", qed_i2, ampl.ph2)
+ + sqrt(omega) * (
+ + einsum("kb,ikja,jb->ia", mp.qed_t1(b.ov), hf.ooov, ampl.ph2)
+ + einsum("jc,ibac,jb->ia", mp.qed_t1(b.ov), hf.ovvv, ampl.ph2))
+ + gs_part * ampl.gs2.to_ndarray()
+ + sqrt(omega) * (
+ - einsum("ib,ab->ia", ampl.ph2, mp.qed_t1_df(b.vv)) # 1. order
+ + einsum("ij,ja->ia", mp.qed_t1_df(b.oo), ampl.ph2) # 1. order
+ - mp.qed_t1_df(b.ov) * ampl.gs2.to_ndarray()) # 1. order
+ ))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_2_couple_edge(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(ph2=(- (omega / 2) * sqrt(2) * (
+ einsum("kc,kc->", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)) * ampl.ph
+ - einsum("ka,kb,ib->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov), ampl.ph)
+ - einsum("ic,jc,ja->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov), ampl.ph)
+ )))
+ return AdcBlock(apply, 0)
+
+
+def block_ph_ph_2_phot_couple_edge(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ def apply(ampl):
+ return AmplitudeVector(ph=(- (omega / 2) * sqrt(2) * (
+ einsum("kc,kc->", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)) * ampl.ph2
+ - einsum("kb,ka,ib->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov), ampl.ph2)
+ - einsum("jc,ic,ja->ia", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov), ampl.ph2)
+ )))
+ return AdcBlock(apply, 0)
+
+
#
# 2nd order coupling
#
+
def block_ph_pphh_2(hf, mp, intermediates):
pia_ooov = intermediates.adc3_pia
pib_ovvv = intermediates.adc3_pib
@@ -346,11 +1041,11 @@ def apply(ampl):
return AmplitudeVector(pphh=(
(
+ einsum("ic,jcab->ijab", ampl.ph, pib_ovvv)
- + einsum("lkic,kc,jlab->ijab", hf.ooov, ampl.ph, mp.t2oo) # 2st
+ + einsum("lkic,kc,jlab->ijab", hf.ooov, ampl.ph, mp.t2oo) # 2nd
).antisymmetrise(0, 1)
+ (
- einsum("ijka,kb->ijab", pia_ooov, ampl.ph)
- - einsum("ijac,kbcd,kd->ijab", mp.t2oo, hf.ovvv, ampl.ph) # 2st
+ - einsum("ijac,kbcd,kd->ijab", mp.t2oo, hf.ovvv, ampl.ph) # 2nd
).antisymmetrise(2, 3)
))
return AdcBlock(apply, 0)
@@ -403,6 +1098,82 @@ def adc2_i2(hf, mp, intermediates):
return hf.foo - 0.5 * einsum("ikab,jkab->ij", mp.t2oo, hf.oovv).symmetrise()
+# qed intermediates for adc2, without the factor of (omega/2),
+# which is added in the actual matrix builder
+@register_as_intermediate
+def adc2_qed_i1(hf, mp, intermediates): # maybe do this with symmetrise
+ return (1 / 2) * (einsum("kb,ka->ab", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov))
+ + einsum("ka,kb->ab", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)))
+
+
+@register_as_intermediate
+def adc2_qed_i2(hf, mp, intermediates): # maybe do this with symmetrise
+ return (1 / 2) * (einsum("jc,ic->ij", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov))
+ + einsum("ic,jc->ij", mp.qed_t1(b.ov), mp.qed_t1_df(b.ov)))
+
+
+@register_as_intermediate
+def qed_adc2_ph_gs_intermediate(hf, mp, intermediates):
+ return (0.5 * einsum("jkib,jkab->ia", hf.ooov, mp.t2oo)
+ + 0.5 * einsum("ijbc,jabc->ia", mp.t2oo, hf.ovvv))
+
+
+@register_as_intermediate
+def term_t2_eri(hf, mp, intermediates):
+ return (einsum("ijab,jkbc->ikac", mp.t2oo, hf.oovv)
+ + einsum("ijab,jkbc->ikac", hf.oovv, mp.t2oo))
+
+
+@register_as_intermediate
+def adc2_qed_ph_ph_2_phot_couple_gs_part(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ return sqrt(omega / 2) * (einsum("ikac,kc->ia", mp.t2oo, mp.qed_t1_df(b.ov)))
+
+
+@register_as_intermediate
+def adc2_qed_ph_ph_2_phot_couple_inner_gs_part(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+
+ return sqrt(omega) * (einsum("ikac,kc->ia", mp.t2oo, mp.qed_t1_df(b.ov)))
+
+
+@register_as_intermediate
+def adc2_qed_couple_i1(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ return (sqrt(omega / 2) * (einsum("kc,kacb->ab", mp.qed_t1(b.ov), hf.ovvv)))
+
+
+@register_as_intermediate
+def adc2_qed_couple_i2(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ return (sqrt(omega / 2) * (einsum("kc,kjic->ij", mp.qed_t1(b.ov), hf.ooov)))
+
+
+@register_as_intermediate
+def adc2_qed_phot_couple_i1(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ return (sqrt(omega / 2) * (einsum("kc,kbca->ab", mp.qed_t1(b.ov), hf.ovvv)))
+
+
+@register_as_intermediate
+def adc2_qed_phot_couple_i2(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ return (sqrt(omega / 2) * (einsum("kc,kijc->ij", mp.qed_t1(b.ov), hf.ooov)))
+
+
+@register_as_intermediate
+def delta_ia_omega(hf, mp, intermediates):
+ omega = float(ReferenceState.get_qed_omega(hf))
+ # Build two Kronecker deltas
+ d_oo = zeros_like(hf.foo)
+ d_vv = zeros_like(hf.fvv)
+ d_oo.set_mask("ii", 1.0)
+ d_vv.set_mask("aa", 1.0)
+
+ return einsum("ii,aa->ia", d_oo, d_vv) * omega
+
+
def adc3_i1(hf, mp, intermediates):
# Used for both CVS and general
td2 = mp.td2(b.oovv)
diff --git a/adcc/adc_pp/state2state_transition_dm.py b/adcc/adc_pp/state2state_transition_dm.py
index f7ea83e7..41793d7b 100644
--- a/adcc/adc_pp/state2state_transition_dm.py
+++ b/adcc/adc_pp/state2state_transition_dm.py
@@ -42,6 +42,88 @@ def s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates):
return dm
+def s2s_tdm_qed_adc2_diag_part(mp, amplitude_l, amplitude_r, intermediates):
+ dm = s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates)
+ ul1 = amplitude_l.ph
+ ur1 = amplitude_r.ph
+ p0_oo = dm.oo.evaluate()
+ p0_vv = dm.vv.evaluate()
+
+ dm_new = OneParticleOperator(mp, is_symmetric=False)
+
+ dm_new.ov = (
+ - einsum("kb,ab->ka", mp.qed_t1(b.ov), p0_vv)
+ + einsum("ji,ic->jc", p0_oo, mp.qed_t1(b.ov))
+ + ul1.dot(mp.qed_t1(b.ov)) * ur1
+ ) / 2
+
+ dm_new.vo = (
+ - einsum("kb,ba->ak", mp.qed_t1(b.ov), p0_vv)
+ + einsum("ij,ic->cj", p0_oo, mp.qed_t1(b.ov))
+ + ur1.dot(mp.qed_t1(b.ov)) * ul1.transpose()
+ ) / 2
+
+ return dm_new
+
+
+def s2s_tdm_qed_adc2_edge_part_couple(mp, amplitude_l, amplitude_r, intermediates):
+ dm = s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates)
+ ul1 = amplitude_l.ph
+ ur1 = amplitude_r.ph
+ p0_oo = dm.oo.evaluate()
+ p0_vv = dm.vv.evaluate()
+
+ dm_new = OneParticleOperator(mp, is_symmetric=False)
+
+ dm_new.ov = (
+ mp.qed_t1(b.ov) * ul1.dot(ur1)
+ - einsum("kb,ab->ka", mp.qed_t1(b.ov), p0_vv)
+ + einsum("ji,ic->jc", p0_oo, mp.qed_t1(b.ov))
+ )
+
+ return dm_new
+
+
+def s2s_tdm_qed_adc2_edge_part_phot_couple(mp, amplitude_l, amplitude_r, intermediates): # noqa: E501
+ dm = s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates)
+ ul1 = amplitude_l.ph
+ ur1 = amplitude_r.ph
+ p0_oo = dm.oo.evaluate()
+ p0_vv = dm.vv.evaluate()
+
+ dm_new = OneParticleOperator(mp, is_symmetric=False)
+
+ dm_new.vo = (
+ einsum("ia->ai", mp.qed_t1(b.ov)) * ul1.dot(ur1)
+ - einsum("kb,ba->ak", mp.qed_t1(b.ov), p0_vv)
+ + einsum("ij,ic->cj", p0_oo, mp.qed_t1(b.ov))
+ )
+
+ return dm_new
+
+
+def s2s_tdm_qed_adc2_ph_pphh_coupl_part(mp, amplitude_l, amplitude_r, intermediates): # noqa: E501
+ ul1 = amplitude_l.ph
+ ur2 = amplitude_r.pphh
+
+ dm = OneParticleOperator(mp, is_symmetric=False)
+
+ dm.ov = -2 * einsum("jb,ijab->ia", ul1, ur2)
+
+ return dm
+
+
+def s2s_tdm_qed_adc2_pphh_ph_phot_coupl_part(mp, amplitude_l, amplitude_r, intermediates): # noqa: E501
+ ul2 = amplitude_l.pphh
+ ur1 = amplitude_r.ph
+
+ dm = OneParticleOperator(mp, is_symmetric=False)
+
+ dm.vo = -2 * einsum("ijab,jb->ai", ul2, ur1)
+
+ return dm
+
+
def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude_l, amplitude_r)
dm = s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates)
@@ -104,6 +186,13 @@ def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
# Ref: https://doi.org/10.1080/00268976.2013.859313
DISPATCH = {"adc0": s2s_tdm_adc0,
"adc1": s2s_tdm_adc0, # same as ADC(0)
+ # The following qed terms are not the actual s2s densities,
+ # but those required for the approx function
+ "qed_adc2_diag": s2s_tdm_qed_adc2_diag_part,
+ "qed_adc2_edge_couple": s2s_tdm_qed_adc2_edge_part_couple,
+ "qed_adc2_edge_phot_couple": s2s_tdm_qed_adc2_edge_part_phot_couple,
+ "qed_adc2_ph_pphh": s2s_tdm_qed_adc2_ph_pphh_coupl_part,
+ "qed_adc2_pphh_ph": s2s_tdm_qed_adc2_pphh_ph_phot_coupl_part,
"adc2": s2s_tdm_adc2,
"adc2x": s2s_tdm_adc2, # same as ADC(2)
}
@@ -143,8 +232,12 @@ def state2state_transition_dm(method, ground_state, amplitude_from,
raise NotImplementedError("state2state_transition_dm is not implemented "
f"for {method.name}.")
else:
- # final state is on the bra side/left (complex conjugate)
- # see ref https://doi.org/10.1080/00268976.2013.859313, appendix A2
- ret = DISPATCH[method.name](ground_state, amplitude_to, amplitude_from,
- intermediates)
+ if hasattr(ground_state, "s2s_contribution"):
+ ret = DISPATCH[ground_state.s2s_contribution](
+ ground_state, amplitude_to, amplitude_from, intermediates)
+ else:
+ # final state is on the bra side/left (complex conjugate)
+ # see ref https://doi.org/10.1080/00268976.2013.859313, appendix A2
+ ret = DISPATCH[method.name](ground_state, amplitude_to, amplitude_from,
+ intermediates)
return ret.evaluate()
diff --git a/adcc/adc_pp/transition_dm.py b/adcc/adc_pp/transition_dm.py
index 111a1f89..8864641c 100644
--- a/adcc/adc_pp/transition_dm.py
+++ b/adcc/adc_pp/transition_dm.py
@@ -117,7 +117,6 @@ def transition_dm(method, ground_state, amplitude, intermediates=None):
"""
Compute the one-particle transition density matrix from ground to excited
state in the MO basis.
-
Parameters
----------
method : str, AdcMethod
@@ -142,5 +141,10 @@ def transition_dm(method, ground_state, amplitude, intermediates=None):
raise NotImplementedError("transition_dm is not implemented "
f"for {method.name}.")
else:
- ret = DISPATCH[method.name](ground_state, amplitude, intermediates)
+ if hasattr(ground_state, "tdm_contribution"):
+ ret = DISPATCH[ground_state.tdm_contribution](
+ ground_state, amplitude, intermediates
+ )
+ else:
+ ret = DISPATCH[method.name](ground_state, amplitude, intermediates)
return ret.evaluate()
diff --git a/adcc/backends/__init__.py b/adcc/backends/__init__.py
index 8abc3424..65efea51 100644
--- a/adcc/backends/__init__.py
+++ b/adcc/backends/__init__.py
@@ -83,7 +83,9 @@ def import_scf_results(res):
import psi4
from . import psi4 as backend_psi4
- if isinstance(res, psi4.core.HF):
+ # Here we also have to import the .core.Wavefunction object,
+ # since this is the one returned by the hilbert package
+ if isinstance(res, (psi4.core.HF, psi4.core.Wavefunction)):
return backend_psi4.import_scf(res)
from libadcc import HartreeFockSolution_i
diff --git a/adcc/backends/psi4.py b/adcc/backends/psi4.py
index c588fa83..b0cf11d1 100644
--- a/adcc/backends/psi4.py
+++ b/adcc/backends/psi4.py
@@ -167,7 +167,17 @@ def get_conv_tol(self):
return threshold
def get_restricted(self):
- return isinstance(self.wfn, (psi4.core.RHF, psi4.core.ROHF))
+ if isinstance(self.wfn, (psi4.core.RHF, psi4.core.ROHF)):
+ return True
+ elif isinstance(self.wfn, psi4.core.UHF):
+ return False
+ else:
+ # The hilbert package returns a more basic object, which does
+ # not provide a restricted indicator, so we determine it here
+ orben_a = np.asarray(self.wfn.epsilon_a())
+ orben_b = np.asarray(self.wfn.epsilon_b())
+ # TODO Maybe give a small tolerance here
+ return all(orben_a == orben_b)
def get_energy_scf(self):
return self.wfn.energy()
@@ -201,17 +211,26 @@ def fill_orbcoeff_fb(self, out):
np.hstack((mo_coeff[0], mo_coeff[1]))
)
- def fill_occupation_f(self, out):
- out[:] = np.hstack((
- np.asarray(self.wfn.occupation_a()),
- np.asarray(self.wfn.occupation_b())
- ))
-
def fill_orben_f(self, out):
orben_a = np.asarray(self.wfn.epsilon_a())
orben_b = np.asarray(self.wfn.epsilon_b())
out[:] = np.hstack((orben_a, orben_b))
+ def fill_occupation_f(self, out):
+ if not hasattr(self.wfn, "occupation_a"): # hilbert package
+ num_of_orbs = self.wfn.nmo()
+ nalpha_elec = self.wfn.nalpha()
+ nbeta_elec = self.wfn.nbeta()
+ occ_array_a = occ_array_b = np.zeros(num_of_orbs)
+ np.put(occ_array_a, np.arange(nalpha_elec), 1)
+ np.put(occ_array_b, np.arange(nbeta_elec), 1)
+ out[:] = np.hstack((occ_array_a, occ_array_b))
+ else:
+ out[:] = np.hstack((
+ np.asarray(self.wfn.occupation_a()),
+ np.asarray(self.wfn.occupation_b())
+ ))
+
def fill_fock_ff(self, slices, out):
diagonal = np.empty(self.n_orbs)
self.fill_orben_f(diagonal)
@@ -231,7 +250,7 @@ def flush_cache(self):
def import_scf(wfn):
- if not isinstance(wfn, psi4.core.HF):
+ if not isinstance(wfn, (psi4.core.HF, psi4.core.Wavefunction)):
raise InvalidReference(
"Only psi4.core.HF and its subtypes are supported references in "
"backends.psi4.import_scf. This indicates that you passed an "
@@ -239,7 +258,7 @@ def import_scf(wfn):
"unrestricted HF calculation."
)
- if not isinstance(wfn, (psi4.core.RHF, psi4.core.UHF)):
+ if not isinstance(wfn, (psi4.core.RHF, psi4.core.UHF, psi4.core.Wavefunction)):
raise InvalidReference("Right now only RHF and UHF references are "
"supported for Psi4.")
@@ -248,7 +267,9 @@ def import_scf(wfn):
# the actual set of options ... theoretically they could differ
scf_type = psi4.core.get_global_option('SCF_TYPE')
# CD = Choleski, DF = density-fitting
- unsupported_scf_types = ["CD", "DISK_DF", "MEM_DF"]
+ unsupported_scf_types = []
+ if not isinstance(wfn, psi4.core.Wavefunction): # hilbert package only uses DF
+ unsupported_scf_types += ["DISK_DF", "MEM_DF", "CD"]
if scf_type in unsupported_scf_types:
raise InvalidReference("Unsupported Psi4 SCF_TYPE, should not be one "
f"of {unsupported_scf_types}")
diff --git a/adcc/guess/guesses_from_diagonal.py b/adcc/guess/guesses_from_diagonal.py
index 2d94f972..5ea478b6 100644
--- a/adcc/guess/guesses_from_diagonal.py
+++ b/adcc/guess/guesses_from_diagonal.py
@@ -33,7 +33,8 @@
def guesses_from_diagonal(matrix, n_guesses, block="ph", spin_change=0,
spin_block_symmetrisation="none",
- degeneracy_tolerance=1e-14, max_diagonal_value=1000):
+ degeneracy_tolerance=1e-14, max_diagonal_value=1000,
+ qed_subblock=None):
"""
Obtain guesses by inspecting a block of the diagonal of the passed ADC
matrix. The symmetry of the returned vectors is already set-up properly.
@@ -90,7 +91,18 @@ def guesses_from_diagonal(matrix, n_guesses, block="ph", spin_change=0,
else:
raise ValueError(f"Don't know how to generate guesses for block {block}")
- return guessfunction(matrix, n_guesses, spin_change,
+ diag = matrix.diagonal()
+
+ if qed_subblock == "phot":
+ diag.ph = diag.ph1
+ if block == "pphh":
+ diag.pphh = diag.pphh1
+ elif qed_subblock == "phot2":
+ diag.ph = diag.ph2
+ if block == "pphh":
+ diag.pphh = diag.pphh2
+
+ return guessfunction(matrix, n_guesses, diag, spin_change,
spin_block_symmetrisation, degeneracy_tolerance,
max_diagonal_value)
@@ -197,7 +209,7 @@ def telem_nospin(telem):
return res
-def guesses_from_diagonal_singles(matrix, n_guesses, spin_change=0,
+def guesses_from_diagonal_singles(matrix, n_guesses, diag, spin_change=0,
spin_block_symmetrisation="none",
degeneracy_tolerance=1e-14,
max_diagonal_value=1000):
@@ -221,7 +233,7 @@ def pred_singles(telem):
and abs(telem.value) <= max_diagonal_value)
elements = find_smallest_matching_elements(
- pred_singles, matrix.diagonal().ph, motrans, n_guesses,
+ pred_singles, diag.ph, motrans, n_guesses,
degeneracy_tolerance=degeneracy_tolerance
)
if len(elements) == 0:
@@ -271,7 +283,7 @@ def telem_nospin(telem):
return [evaluate(v / np.sqrt(v @ v)) for v in ret[:ivec]]
-def guesses_from_diagonal_doubles(matrix, n_guesses, spin_change=0,
+def guesses_from_diagonal_doubles(matrix, n_guesses, diag, spin_change=0,
spin_block_symmetrisation="none",
degeneracy_tolerance=1e-14,
max_diagonal_value=1000):
@@ -299,7 +311,7 @@ def guesses_from_diagonal_doubles(matrix, n_guesses, spin_change=0,
ret = ret[:n_found]
# Filter out elements above the noted diagonal value
- diagonal_elements = [ret_d.pphh.dot(matrix.diagonal().pphh * ret_d.pphh)
+ diagonal_elements = [ret_d.pphh.dot(diag.pphh * ret_d.pphh)
for ret_d in ret]
return [ret[i] for (i, elem) in enumerate(diagonal_elements)
if elem <= max_diagonal_value]
diff --git a/adcc/qed_matrix_from_diag_adc.py b/adcc/qed_matrix_from_diag_adc.py
new file mode 100644
index 00000000..d3830049
--- /dev/null
+++ b/adcc/qed_matrix_from_diag_adc.py
@@ -0,0 +1,178 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2018 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
+#import scipy.linalg as sp
+
+
+class qed_matrix_from_diag_adc:
+ def __init__(self, exstates, refstate):
+ self.s2s = exstates.s2s_dipole_moments_qed
+ self.tdm = exstates.transition_dipole_moments_qed
+ self.h1 = exstates.qed_second_order_ph_ph_couplings
+ self.coupl = np.array(refstate.coupling)
+ self.freq = np.linalg.norm(refstate.frequency)
+ self.full_freq = np.linalg.norm(np.real(refstate.freq_with_loss)) +\
+ np.linalg.norm(np.imag(refstate.freq_with_loss)) * 1j
+ self.n_adc = len(exstates.excitation_energy)
+ self.exc_en = exstates.excitation_energy
+
+ def loop_helper(self, s2s_tensor):
+ ret = np.array([[self.coupl.dot(s2s)
+ for s2s in s2s_block]
+ for s2s_block in s2s_tensor])
+ ret *= - np.sqrt(self.freq / 2) * np.sqrt(2 * self.freq)
+ return ret
+
+ def first_order_coupling(self):
+
+ # build the blocks of the matrix
+
+ tdm_block = np.empty(self.n_adc)
+
+ for i, tdm in enumerate(self.tdm):
+ tdm_block[i] = self.coupl.dot(tdm)
+ tdm_block *= np.sqrt(2 * self.freq)
+
+ s2s_block = self.loop_helper(self.s2s["qed_adc1_off_diag"])
+
+ tdm_block *= - np.sqrt(self.freq / 2)
+
+ if np.iscomplex(self.full_freq):
+ self.freq = self.full_freq
+
+ elec_block = np.diag(self.exc_en)
+ phot_block = np.diag(self.exc_en + self.freq)
+
+ # build the matrix
+
+ matrix_upper = np.vstack((elec_block, tdm_block.reshape((1, self.n_adc)),
+ s2s_block))
+ matrix_middle = np.concatenate((tdm_block, np.array([self.freq]),
+ np.zeros(self.n_adc)))
+ matrix_lower = np.vstack((s2s_block, np.zeros((1, self.n_adc)),
+ phot_block))
+
+ matrix = np.hstack((matrix_upper,
+ matrix_middle.reshape((len(matrix_middle), 1)),
+ matrix_lower))
+
+ #if np.iscomplex(self.full_freq):
+ # return sp.eig(matrix)
+ #else:
+ # return sp.eigh(matrix)
+ return matrix
+
+ def second_order_coupling(self):
+
+ # tdm part
+
+ qed_adc2_tdm_vec = np.empty(self.n_adc)
+
+ for i, tdm in enumerate(self.tdm):
+ qed_adc2_tdm_vec[i] = self.coupl.dot(tdm)
+ qed_adc2_tdm_vec *= - np.sqrt(self.freq / 2) * np.sqrt(2 * self.freq)
+
+ # s2s_dipole parts of the ph_ph blocks
+
+ qed_adc1_off_diag_block = self.loop_helper(self.s2s["qed_adc1_off_diag"])
+
+ qed_adc2_diag_block = self.loop_helper(self.s2s["qed_adc2_diag"])
+ # missing factor from state.s2s_dipole_moments_qed_adc2_diag
+ # TODO: commit to one way of defining these factors
+ # within the approx method
+ qed_adc2_diag_block *= np.sqrt(self.freq / 2)
+
+ qed_adc2_edge_block_couple = self.loop_helper(
+ self.s2s["qed_adc2_edge_couple"])
+ # missing factor from state.s2s_dipole_moments_qed_adc2_edge
+ qed_adc2_edge_block_couple *= np.sqrt(self.freq)
+
+ qed_adc2_edge_block_phot_couple = self.loop_helper(
+ self.s2s["qed_adc2_edge_phot_couple"])
+ # missing factor from state.s2s_dipole_moments_qed_adc2_edge
+ qed_adc2_edge_block_phot_couple *= np.sqrt(self.freq)
+
+ # s2s_dipole parts of the pphh_ph and ph_pphh blocks
+
+ qed_adc2_ph_pphh_couple_block = self.loop_helper(
+ self.s2s["qed_adc2_ph_pphh"])
+
+ qed_adc2_pphh_ph_phot_couple_block = self.loop_helper(
+ self.s2s["qed_adc2_pphh_ph"])
+
+ # we still need the H_1 expectation value "as property"
+
+ qed_adc2_couple_block = np.sqrt(self.freq / 2) *\
+ self.h1["couple"]
+ qed_adc2_phot_couple_block = np.sqrt(self.freq / 2) *\
+ self.h1["phot_couple"]
+
+ # build the blocks of the matrix
+
+ if np.iscomplex(self.full_freq):
+ self.freq = self.full_freq
+
+ single_excitation_states = np.ones(self.n_adc)
+
+ elec_block = np.diag(self.exc_en) + qed_adc2_diag_block
+
+ phot_block = np.diag(self.exc_en) + qed_adc2_diag_block * 2 +\
+ np.diag(single_excitation_states) * self.freq
+
+ phot2_block = np.diag(self.exc_en) + qed_adc2_diag_block * 3 +\
+ np.diag(single_excitation_states) * 2 * self.freq
+
+ couple_block = qed_adc1_off_diag_block + qed_adc2_ph_pphh_couple_block +\
+ qed_adc2_couple_block
+
+ phot_couple_block = qed_adc1_off_diag_block +\
+ qed_adc2_pphh_ph_phot_couple_block +\
+ qed_adc2_phot_couple_block
+
+ # build the matrix
+
+ matrix_1 = np.vstack((elec_block, qed_adc2_tdm_vec.reshape((1, self.n_adc)), # noqa: E501
+ phot_couple_block, np.zeros((1, self.n_adc)),
+ qed_adc2_edge_block_phot_couple))
+ matrix_2 = np.concatenate((qed_adc2_tdm_vec, np.array([self.freq]),
+ np.zeros(self.n_adc), np.array([0]),
+ np.zeros(self.n_adc)))
+ matrix_3 = np.vstack((couple_block, np.zeros((1, self.n_adc)), phot_block,
+ np.sqrt(2) * qed_adc2_tdm_vec.reshape((1, self.n_adc)), # noqa: E501
+ np.sqrt(2) * phot_couple_block))
+ matrix_4 = np.concatenate((np.zeros(self.n_adc), np.array([0]),
+ np.sqrt(2) * qed_adc2_tdm_vec,
+ 2 * np.array([self.freq]), np.zeros(self.n_adc)))
+ matrix_5 = np.vstack((qed_adc2_edge_block_couple, np.zeros((1, self.n_adc)),
+ np.sqrt(2) * couple_block, np.zeros((1, self.n_adc)),
+ phot2_block))
+
+ matrix = np.hstack((matrix_1, matrix_2.reshape((len(matrix_2), 1)),
+ matrix_3, matrix_4.reshape((len(matrix_4), 1)),
+ matrix_5))
+
+ #if np.iscomplex(self.full_freq):
+ # return sp.eig(matrix)
+ #else:
+ # return sp.eigh(matrix)
+ return matrix
diff --git a/adcc/solver/explicit_symmetrisation.py b/adcc/solver/explicit_symmetrisation.py
index c11325da..81229482 100644
--- a/adcc/solver/explicit_symmetrisation.py
+++ b/adcc/solver/explicit_symmetrisation.py
@@ -49,9 +49,7 @@ def __init__(self, matrix):
def symmetrise(self, new_vectors):
"""
Symmetrise a set of new vectors to be added to the subspace.
-
new_vectors Vectors to symmetrise (updated in-place)
-
Returns:
The updated new_vectors
"""
diff --git a/adcc/solver/preconditioner.py b/adcc/solver/preconditioner.py
index 07cc1c0c..7f55dac6 100644
--- a/adcc/solver/preconditioner.py
+++ b/adcc/solver/preconditioner.py
@@ -75,7 +75,6 @@ def apply(self, invecs):
"with number of shifts stored inside "
"precoditioner. Update using the "
"'update_shifts' method.")
-
return [v / (self.diagonal - self.shifts[i])
for i, v in enumerate(invecs)]
else:
diff --git a/adcc/test_ExcitedStates.py b/adcc/test_ExcitedStates.py
index b4d1e0b3..9c13ef8e 100644
--- a/adcc/test_ExcitedStates.py
+++ b/adcc/test_ExcitedStates.py
@@ -41,8 +41,11 @@ def base_test(self, system, method, kind):
for key in dir(exci):
if key.startswith("_"):
continue
+ # TODO: Remove qed from blacklist, if the function
+ # qed_second_order_ph_ph_couplings, is removed from
+ # @mark_excitation_property()
blacklist = ["__", "index", "_ao", "excitation_vector",
- "method", "parent_state"]
+ "method", "parent_state", "qed"]
if any(b in key for b in blacklist):
continue
try:
diff --git a/adcc/test_qed.py b/adcc/test_qed.py
new file mode 100644
index 00000000..9da52fe2
--- /dev/null
+++ b/adcc/test_qed.py
@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2018 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 unittest
+
+from adcc.misc import expand_test_templates
+import adcc
+
+from numpy.testing import assert_allclose
+
+from adcc.testdata.cache import cache
+from adcc.testdata.cache import qed_data
+
+import itertools
+
+# In principle one could also test the approx method against the full
+# method, by expanding them to the full matrix dimension. The smallest
+# example would be HF sto-3g, since it contains only one virtual orbital,
+# which keeps the matrix dimension low. This is important since we
+# build most of the matrix in the approx method from properties, so this
+# is very slow for a lot of states, compared to the standard method.
+# However, even this test case would take quite some time, so we leave
+# it out for now.
+
+testcases = ["methox_sto3g", "h2o_sto3g"]
+methods = ["adc2"]
+
+
+@expand_test_templates(list(itertools.product(testcases, methods)))
+class qed_test(unittest.TestCase):
+ def set_refstate(self, case):
+ # Here we have to build the new refstate, instead of using the cached
+ # refstate, because otherwise the qed objects are not build
+ self.refstate = adcc.ReferenceState(cache.hfdata[case])
+ self.refstate.coupling = [0.0, 0.0, 0.05]
+ self.refstate.frequency = [0.0, 0.0, 0.5]
+ self.refstate.freq_with_loss = self.refstate.frequency
+ self.refstate.qed_hf = True
+ self.refstate.is_qed = True
+
+ def template_approx(self, case, method):
+ self.set_refstate(case)
+ self.refstate.approx = True
+
+ approx = adcc.adc2(self.refstate, n_singlets=5, conv_tol=1e-7)
+
+ ref_name = f"{case}_{method}_approx"
+ approx_ref = qed_data[ref_name]["excitation_energy"]
+
+ assert_allclose(approx.qed_excitation_energy,
+ approx_ref, atol=1e-6)
+
+ def template_full(self, case, method):
+ self.set_refstate(case)
+
+ full = adcc.adc2(self.refstate, n_singlets=3, conv_tol=1e-7)
+
+ ref_name = f"{case}_{method}_full"
+ full_ref = qed_data[ref_name]["excitation_energy"]
+
+ assert_allclose(full.excitation_energy,
+ full_ref, atol=1e-6)
diff --git a/adcc/testdata/cache.py b/adcc/testdata/cache.py
index 1f0e37f1..3880dbb1 100644
--- a/adcc/testdata/cache.py
+++ b/adcc/testdata/cache.py
@@ -260,3 +260,4 @@ def read_yaml_data(fname):
tmole_data = read_yaml_data("tmole_dump.yml")
psi4_data = read_yaml_data("psi4_dump.yml")
pyscf_data = read_yaml_data("pyscf_dump.yml")
+qed_data = read_yaml_data("qed_dump.yml")
diff --git a/adcc/testdata/qed_dump.yml b/adcc/testdata/qed_dump.yml
new file mode 100644
index 00000000..e19b0812
--- /dev/null
+++ b/adcc/testdata/qed_dump.yml
@@ -0,0 +1,28 @@
+h2o_sto3g_adc2_full:
+ basis: sto3g
+ method: adc2
+ molecule: h2o
+ approx: False
+ excitation_energy: [0.469789, 0.496679, 0.571507]
+h2o_sto3g_adc2_approx:
+ basis: sto3g
+ method: adc2
+ molecule: h2o
+ approx: True
+ excitation_energy: [0.469793, 0.496957, 0.571511, 0.593613, 0.712661, 0.844501,
+ 0.970428, 0.993996, 1.0714 , 1.094744, 1.212892, 1.346747,
+ 1.47356 , 1.574525, 1.597913, 1.716008, 1.841888]
+methox_sto3g_adc2_full:
+ basis: sto3g
+ method: adc2
+ molecule: cn
+ approx: False
+ excitation_energy: [0.341541, 0.415104, 0.486678]
+methox_sto3g_adc2_approx:
+ basis: sto3g
+ method: adc2
+ molecule: cn
+ approx: True
+ excitation_energy: [0.341774, 0.41538 , 0.491732, 0.501357, 0.511708, 0.546447,
+ 0.842004, 0.915563, 0.990758, 1.002418, 1.011967, 1.046802,
+ 1.342381, 1.416067, 1.493567, 1.512128, 1.54733]
\ No newline at end of file
diff --git a/adcc/workflow.py b/adcc/workflow.py
index a290032e..b672a90d 100644
--- a/adcc/workflow.py
+++ b/adcc/workflow.py
@@ -22,8 +22,10 @@
## ---------------------------------------------------------------------
import sys
import warnings
+import numpy as np
+import scipy.linalg as sp
-from libadcc import ReferenceState
+from libadcc import ReferenceState, set_lt_scalar
from . import solver
from .guess import (guesses_any, guesses_singlet, guesses_spin_flip,
@@ -38,6 +40,8 @@
from .solver.davidson import jacobi_davidson
from .solver.explicit_symmetrisation import (IndexSpinSymmetrisation,
IndexSymmetrisation)
+from .AmplitudeVector import AmplitudeVector
+from .qed_matrix_from_diag_adc import qed_matrix_from_diag_adc
__all__ = ["run_adc"]
@@ -47,7 +51,8 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None,
n_guesses_doubles=None, output=sys.stdout, core_orbitals=None,
frozen_core=None, frozen_virtual=None, method=None,
n_singlets=None, n_triplets=None, n_spin_flip=None,
- environment=None, **solverargs):
+ environment=None, coupl=None, freq=None, qed_hf=True,
+ qed_approx=False, **solverargs):
"""Run an ADC calculation.
Main entry point to run an ADC calculation. The reference to build the ADC
@@ -133,6 +138,23 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None,
The keywords to specify how coupling to an environment model,
e.g. PE, is treated. For details see :ref:`environment`.
+ coupl : list or tuple, optional
+ Specifies the coupling for a qed calculation as [x, y, z]
+
+ freq : list or tuple, optional
+ Specifies the photon energy corresponding to the coupling specified
+ in "coupl", which is required for a qed calculation, as [x, y, z]
+
+ qed_hf : bool, optional
+ Specifies, if the mean-field solution to the Pauli-Fierz Hamiltonian
+ is provided for the qed calculation. False expects the standard HF
+ solution.
+
+ qed_approx : bool, optional
+ Indicates whether the approximate solution for the qed method should
+ be calculated. (After paper is published, put link here)
+
+
Other parameters
----------------
max_subspace : int, optional
@@ -178,11 +200,13 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None,
"""
matrix = construct_adcmatrix(
data_or_matrix, core_orbitals=core_orbitals, frozen_core=frozen_core,
- frozen_virtual=frozen_virtual, method=method)
+ frozen_virtual=frozen_virtual, method=method, coupl=coupl,
+ freq=freq, qed_hf=qed_hf, qed_approx=qed_approx)
n_states, kind = validate_state_parameters(
matrix.reference_state, n_states=n_states, n_singlets=n_singlets,
- n_triplets=n_triplets, n_spin_flip=n_spin_flip, kind=kind)
+ n_triplets=n_triplets, n_spin_flip=n_spin_flip, kind=kind, coupl=coupl,
+ freq=freq, qed_hf=qed_hf, qed_approx=qed_approx)
# Determine spin change during excitation. If guesses is not None,
# i.e. user-provided, we cannot guarantee for obtaining a particular
@@ -212,6 +236,31 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None,
# add environment corrections to excited states
exstates += env_energy_corrections
+
+ # Build QED (approximated) matrix from "standard" ADC matrix
+ # and expectation values, if requested
+ if matrix.reference_state.approx:
+ qed_matrix = qed_matrix_from_diag_adc(exstates, matrix.reference_state)
+ if method == "adc2":
+ #qed_eigvals, qed_eigvecs = qed_matrix.second_order_coupling()
+ qed_approx_matrix = qed_matrix.second_order_coupling()
+ else:
+ #qed_eigvals, qed_eigvecs = qed_matrix.first_order_coupling()
+ qed_approx_matrix = qed_matrix.first_order_coupling()
+ # TODO: This is a bad solution, but filtering out the final excitation
+ # vectors and properly feeding them back into the corresponding libtensor is
+ # usually just unnecessary, especially since consistent qed properties
+ # are not implemented yet. This way the user is at least provided with
+ # the ADC result without polaritonic coupling between the states and can
+ # then request the energies and vectors from the exstates object.
+ exstates.qed_approx_matrix = qed_approx_matrix
+ if any(np.iscomplex(freq)):
+ qed_eigvals, qed_eigvecs = sp.eig(qed_approx_matrix)
+ else:
+ qed_eigvals, qed_eigvecs = sp.eigh(qed_approx_matrix)
+ exstates.qed_excitation_energy = qed_eigvals
+ exstates.qed_excitation_vector = qed_eigvecs
+
return exstates
@@ -219,7 +268,8 @@ def run_adc(data_or_matrix, n_states=None, kind="any", conv_tol=None,
# Individual steps
#
def construct_adcmatrix(data_or_matrix, core_orbitals=None, frozen_core=None,
- frozen_virtual=None, method=None):
+ frozen_virtual=None, method=None, coupl=None,
+ freq=None, qed_hf=True, qed_approx=False):
"""
Use the provided data or AdcMatrix object to check consistency of the
other passed parameters and construct the AdcMatrix object representing
@@ -242,10 +292,17 @@ def construct_adcmatrix(data_or_matrix, core_orbitals=None, frozen_core=None,
"to be specified via the parameter "
"core_orbitals.")
try:
+ if coupl is not None:
+ is_qed = True
+ else:
+ is_qed = False
refstate = adcc_ReferenceState(data_or_matrix,
core_orbitals=core_orbitals,
frozen_core=frozen_core,
- frozen_virtual=frozen_virtual)
+ frozen_virtual=frozen_virtual,
+ is_qed=is_qed, coupl=coupl,
+ qed_hf=qed_hf, freq=freq,
+ qed_approx=qed_approx)
except ValueError as e:
raise InputError(str(e)) # In case of an issue with the spaces
data_or_matrix = refstate
@@ -283,7 +340,9 @@ def construct_adcmatrix(data_or_matrix, core_orbitals=None, frozen_core=None,
def validate_state_parameters(reference_state, n_states=None, n_singlets=None,
- n_triplets=None, n_spin_flip=None, kind="any"):
+ n_triplets=None, n_spin_flip=None, kind="any",
+ coupl=None, freq=None, qed_hf=True,
+ qed_approx=False):
"""
Check the passed state parameters for consistency with itself and with
the passed reference and normalise them. In the end return the number of
@@ -342,6 +401,21 @@ def validate_state_parameters(reference_state, n_states=None, n_singlets=None,
raise InputError("kind==spin_flip is only valid for "
"ADC calculations in combination with an unrestricted "
"ground state.")
+
+ # qed sanity checks
+ if coupl is not None or freq is not None:
+ if not (coupl is not None and freq is not None):
+ raise InputError("qed calculation requires coupl and freq")
+ if len(coupl) != 3: # or len(freq) != 3:
+ raise InputError("freq and coupl must contain 3 elements, "
+ "i.e. x, y, z")
+ if any(np.iscomplex(freq)) and not qed_approx:
+ raise InputError("imaginary contribution to freq will only be "
+ "processed, if qed_approx=True")
+ if not qed_hf and qed_approx:
+ raise InputError("Approximate QED ADC method is only available for "
+ "QED-HF reference")
+
return n_states, kind
@@ -392,8 +466,11 @@ def diagonalise_adcmatrix(matrix, n_states, kind, eigensolver="davidson",
if guesses is None:
if n_guesses is None:
n_guesses = estimate_n_guesses(matrix, n_states, n_guesses_per_state)
- guesses = obtain_guesses_by_inspection(matrix, n_guesses, kind,
- n_guesses_doubles)
+ if matrix.reference_state.is_qed and not matrix.reference_state.approx:
+ guess_function = obtain_guesses_by_inspection_qed
+ else:
+ guess_function = obtain_guesses_by_inspection
+ guesses = guess_function(matrix, n_guesses, kind, n_guesses_doubles)
else:
if len(guesses) < n_states:
raise InputError("Less guesses provided via guesses (== {}) "
@@ -448,10 +525,10 @@ def estimate_n_guesses(matrix, n_states, singles_only=True,
return max(n_states, n_guesses)
-def obtain_guesses_by_inspection(matrix, n_guesses, kind, n_guesses_doubles=None):
+def obtain_guesses_by_inspection(matrix, n_guesses, kind, n_guesses_doubles=None, qed_subblock=None): # noqa: E501
"""
Obtain guesses by inspecting the diagonal matrix elements.
- If n_guesses_doubles is not None, this is number is always adhered to.
+ If n_guesses_doubles is not None, this number is always adhered to.
Otherwise the number of doubles guesses is adjusted to fill up whatever
the singles guesses cannot provide to reach n_guesses.
Internal function called from run_adc.
@@ -471,7 +548,8 @@ def obtain_guesses_by_inspection(matrix, n_guesses, kind, n_guesses_doubles=None
n_guess_singles = n_guesses
if n_guesses_doubles is not None:
n_guess_singles = n_guesses - n_guesses_doubles
- singles_guesses = guess_function(matrix, n_guess_singles, block="ph")
+ singles_guesses = guess_function(matrix, n_guess_singles, block="ph",
+ qed_subblock=qed_subblock)
doubles_guesses = []
if "pphh" in matrix.axis_blocks:
@@ -481,7 +559,8 @@ def obtain_guesses_by_inspection(matrix, n_guesses, kind, n_guesses_doubles=None
n_guesses_doubles = n_guesses - len(singles_guesses)
if n_guesses_doubles > 0:
doubles_guesses = guess_function(matrix, n_guesses_doubles,
- block="pphh")
+ block="pphh",
+ qed_subblock=qed_subblock)
total_guesses = singles_guesses + doubles_guesses
if len(total_guesses) < n_guesses:
@@ -490,6 +569,65 @@ def obtain_guesses_by_inspection(matrix, n_guesses, kind, n_guesses_doubles=None
return total_guesses
+def obtain_guesses_by_inspection_qed(matrix, n_guesses, kind, n_guesses_doubles=None): # noqa: E501
+ """
+ Obtain guesses for full qed method, by pushing the subblocks
+ into obtain_guesses_by_inspection.
+ Internal function called from run_adc.
+ """
+ guesses_elec = obtain_guesses_by_inspection(
+ matrix, n_guesses, kind, n_guesses_doubles, qed_subblock=None)
+ guesses_phot = obtain_guesses_by_inspection(
+ matrix, n_guesses, kind, n_guesses_doubles, qed_subblock="phot")
+ guesses_phot2 = obtain_guesses_by_inspection(
+ matrix, n_guesses, kind, n_guesses_doubles, qed_subblock="phot2")
+
+ # Usually only few states are requested and most of them are close
+ # to pure electronic states, so we initialize the guess vectors
+ # as almost purely electric guesses.
+ for i in np.arange(n_guesses):
+ # TODO: maybe make these values accessible by a keyword, since
+ # they can tune the performance. From my experience these work
+ # very well though
+ guesses_phot[i] *= 0.02
+ guesses_phot2[i] *= 0.001
+ if n_guesses != len(guesses_phot):
+ raise InputError("amount of guesses for electronic and photonic must be "
+ "equal, but are {} electronic and {} photonic "
+ "guesses".format(len(guesses_elec), len(guesses_phot)))
+
+ zero = set_lt_scalar(0.0)
+
+ if hasattr(guesses_elec[0], "pphh"):
+ final_guesses = [AmplitudeVector(**{
+ "ph": guesses_elec[guess_index].ph,
+ "pphh": guesses_elec[guess_index].pphh,
+ "gs1": zero.copy(), "ph1": guesses_phot[guess_index].ph,
+ "pphh1": guesses_phot[guess_index].pphh,
+ "gs2": zero.copy(), "ph2": guesses_phot2[guess_index].ph,
+ "pphh2": guesses_phot2[guess_index].pphh
+ }) for guess_index in np.arange(n_guesses)]
+ else:
+ final_guesses = [AmplitudeVector(**{
+ "ph": guesses_elec[guess_index].ph,
+ "gs1": zero.copy(), "ph1": guesses_phot[guess_index].ph,
+ "gs2": zero.copy(), "ph2": guesses_phot2[guess_index].ph
+ }) for guess_index in np.arange(n_guesses)]
+
+ # TODO: maybe make these values accessible by a keyword, since
+ # they can tune the performance. From my experience these work
+ # quite well though, but depending on how strong a state couples
+ # to the single photon dispersion mode and how many states one
+ # requests, adjusting these values can increase the
+ # convergence rate
+ # for stronger coupling e.g. 2
+ final_guesses[n_guesses - 2].gs1.set_from_ndarray(np.array([5]))
+ # for stronger coupling e.g. 5
+ final_guesses[n_guesses - 1].gs2.set_from_ndarray(np.array([20]))
+
+ return [vec / np.sqrt(vec @ vec) for vec in final_guesses]
+
+
def setup_solver_printing(solmethod_name, matrix, kind, default_print,
output=None):
"""
diff --git a/libadcc/pyiface/export_Tensor.cc b/libadcc/pyiface/export_Tensor.cc
index ecff8ac7..d8ef1bf0 100644
--- a/libadcc/pyiface/export_Tensor.cc
+++ b/libadcc/pyiface/export_Tensor.cc
@@ -23,6 +23,9 @@
#include
#include
#include
+#include
+#include "../AdcMemory.hh"
+#include "../tests/wrap_libtensor.hh"
namespace libadcc {
@@ -346,6 +349,17 @@ static py::object Tensor___repr__(const Tensor& self) {
return Tensor___str__(self);
}
+static ten_ptr set_lt_scalar(const py::float_ n) {
+ auto adcmem_ptr = std::shared_ptr(new AdcMemory());
+ libtensor::bispace<1> bis(1);
+ std::vector ax{{"x", 1}};
+ double vector[] = {n};
+ auto v_ptr = std::make_shared>(bis);
+ libtensor::btod_import_raw<1>(vector, bis.get_bis().get_dims()).perform(*v_ptr);
+ std::shared_ptr tensor_ptr = wrap_libtensor(adcmem_ptr, ax, v_ptr);
+ return tensor_ptr;
+}
+
//
// Operations with a scalar
//
@@ -540,6 +554,9 @@ void export_Tensor(py::module& m) {
m.def("trace", &Tensor_trace_2, "tensor"_a);
m.def("linear_combination_strict", &linear_combination_strict, "coefficients"_a,
"tensors"_a);
+ // This is necessary for smooth handling of ground state contributions
+ // in AmplitudeVectors
+ m.def("set_lt_scalar", &set_lt_scalar, "a"_a);
}
} // namespace libadcc