Transition dipole moments for asymmetric operators (#158)

* fix tdm for asymmetric operators

* generalize mtms for asymmetric operators

* small changes bmatrix

* update tests

* fix cvs tdms

* flake8 code cleanup

* pytest.skip and rename operators

* update dump adcc

* make use of State2States in test_IsrMatrix

* flake8

* add adcc_reference data

* use cached properties

* add tests against adcc reference data

* reviews

* reviews: reuse lower-order functions

* new test data

* remove np pins

* try newer macos runner

---------

Co-authored-by: Maximilian Scheurer <max.scheurer@me.com>

.github/workflows/ci.yaml

@@ -19,7 +19,7 @@ jobs:
         include:
           - {version: '3.7', os: ubuntu-latest, documentation: True}
           - {version: '3.9', os: ubuntu-latest, documentation: False}
-          - {version: '3.7', os: macOS-10.15  , documentation: False}
+          - {version: '3.7', os: macos-11 , documentation: False}
     steps:
       - uses: actions/checkout@v2
       - uses: actions/setup-python@v2

adcc/IsrMatrix.py

@@ -42,7 +42,7 @@ class IsrMatrix(AdcMatrixlike):
         "adc3":  dict(ph_ph=3, ph_pphh=2,    pphh_ph=2,    pphh_pphh=1),     # noqa: E501
     }
 
-    def __init__(self, method, hf_or_mp, operators, block_orders=None):
+    def __init__(self, method, hf_or_mp, operator, block_orders=None):
         """
         Initialise an ISR matrix of a given one-particle operator
         for the provided ADC method.
@@ -53,7 +53,7 @@ def __init__(self, method, hf_or_mp, operators, block_orders=None):
             Method to use.
         hf_or_mp : adcc.ReferenceState or adcc.LazyMp
             HF reference or MP ground state.
-        operators : adcc.OneParticleOperator or list of adcc.OneParticleOperator
+        operator : adcc.OneParticleOperator or list of adcc.OneParticleOperator
                     objects
             One-particle matrix elements associated with a one-particle operator.
         block_orders : optional
@@ -71,12 +71,12 @@ def __init__(self, method, hf_or_mp, operators, block_orders=None):
         if not isinstance(method, AdcMethod):
             method = AdcMethod(method)
 
-        if not isinstance(operators, list):
-            self.operators = [operators]
+        if not isinstance(operator, list):
+            self.operator = [operator]
         else:
-            self.operators = operators.copy()
-        if not all(isinstance(op, OneParticleOperator) for op in self.operators):
-            raise TypeError("operators is not a valid object. It needs to be "
+            self.operator = operator.copy()
+        if not all(isinstance(op, OneParticleOperator) for op in self.operator):
+            raise TypeError("operator is not a valid object. It needs to be "
                             "either an OneParticleOperator or a list of "
                             "OneParticleOperator objects.")
 
@@ -115,11 +115,11 @@ def __init__(self, method, hf_or_mp, operators, block_orders=None):
             if self.is_core_valence_separated:
                 variant = "cvs"
             blocks = [{
-                block: ppbmatrix.block(self.ground_state, operator,
+                block: ppbmatrix.block(self.ground_state, op,
                                        block.split("_"), order=order,
                                        variant=variant)
                 for block, order in self.block_orders.items() if order is not None
-            } for operator in self.operators]
+            } for op in self.operator]
             # TODO Rename to self.block in 0.16.0
             self.blocks_ph = [{
                 b: bl[b].apply for b in bl
@@ -145,10 +145,10 @@ def matvec(self, v):
 
     def rmatvec(self, v):
         # Hermitian operators
-        if all(op.is_symmetric for op in self.operators):
+        if all(op.is_symmetric for op in self.operator):
             return self.matvec(v)
         else:
-            diffv = [op.ov + op.vo.transpose((1, 0)) for op in self.operators]
+            diffv = [op.ov + op.vo.transpose((1, 0)) for op in self.operator]
             # anti-Hermitian operators
             if all(dv.dot(dv) < 1e-12 for dv in diffv):
                 return [

adcc/adc_pp/bmatrix.py

@@ -122,9 +122,9 @@ def block_pphh_ph_0(ground_state, op):
 #
 def block_ph_pphh_1(ground_state, op):
     if op.is_symmetric:
-        op_vo = op.ov.transpose((1, 0))
+        op_vo = op.ov.transpose()
     else:
-        op_vo = op.vo.copy()
+        op_vo = op.vo
     t2 = ground_state.t2(b.oovv)
 
     def apply(ampl):
@@ -141,9 +141,9 @@ def apply(ampl):
 
 def block_pphh_ph_1(ground_state, op):
     if op.is_symmetric:
-        op_vo = op.ov.transpose((1, 0))
+        op_vo = op.ov.transpose()
     else:
-        op_vo = op.vo.copy()
+        op_vo = op.vo
     t2 = ground_state.t2(b.oovv)
 
     def apply(ampl):
@@ -175,9 +175,9 @@ def apply(ampl):
 #
 def block_ph_ph_2(ground_state, op):
     if op.is_symmetric:
-        op_vo = op.ov.transpose((1, 0))
+        op_vo = op.ov.transpose()
     else:
-        op_vo = op.vo.copy()
+        op_vo = op.vo
     p0 = ground_state.mp2_diffdm
     t2 = ground_state.t2(b.oovv)
 

adcc/adc_pp/modified_transition_moments.py

@@ -30,48 +30,55 @@
 from adcc.AmplitudeVector import AmplitudeVector
 
 
-def mtm_adc0(mp, dipop, intermediates):
-    return AmplitudeVector(ph=dipop.ov)
+def mtm_adc0(mp, op, intermediates):
+    f1 = op.ov if op.is_symmetric else op.vo.transpose()
+    return AmplitudeVector(ph=f1)
 
 
-def mtm_adc1(mp, dipop, intermediates):
-    f1 = dipop.ov - einsum("ijab,jb->ia", mp.t2(b.oovv), dipop.ov)
-    return AmplitudeVector(ph=f1)
+def mtm_adc1(mp, op, intermediates):
+    ampl = mtm_adc0(mp, op, intermediates)
+    f1 = - 1.0 * einsum("ijab,jb->ia", mp.t2(b.oovv), op.ov)
+    return ampl + AmplitudeVector(ph=f1)
 
 
-def mtm_adc2(mp, dipop, intermediates):
+def mtm_adc2(mp, op, intermediates):
     t2 = mp.t2(b.oovv)
     p0 = mp.mp2_diffdm
 
+    op_vo = op.ov.transpose() if op.is_symmetric else op.vo
+
+    ampl = mtm_adc1(mp, op, intermediates)
     f1 = (
-        + dipop.ov
-        - einsum("ijab,jb->ia", t2,
-                 + dipop.ov - 0.5 * einsum("jkbc,kc->jb", t2, dipop.ov))
-        + 0.5 * einsum("ij,ja->ia", p0.oo, dipop.ov)
-        - 0.5 * einsum("ib,ab->ia", dipop.ov, p0.vv)
-        + einsum("ib,ab->ia", p0.ov, dipop.vv)
-        - einsum("ij,ja->ia", dipop.oo, p0.ov)
-        - einsum("ijab,jb->ia", mp.td2(b.oovv), dipop.ov)
+        + 0.5 * einsum("ijab,jkbc,ck->ia", t2, t2, op_vo)
+        + 0.5 * einsum("ij,aj->ia", p0.oo, op_vo)
+        - 0.5 * einsum("bi,ab->ia", op_vo, p0.vv)
+        + 1.0 * einsum("ib,ab->ia", p0.ov, op.vv)
+        - 1.0 * einsum("ji,ja->ia", op.oo, p0.ov)
+        - 1.0 * einsum("ijab,jb->ia", mp.td2(b.oovv), op.ov)
     )
     f2 = (
-        + einsum("ijac,bc->ijab", t2, dipop.vv).antisymmetrise(2, 3)
-        - einsum("ik,kjab->ijab", dipop.oo, t2).antisymmetrise(0, 1)
+        + 1.0 * einsum("ijac,bc->ijab", t2, op.vv).antisymmetrise(2, 3)
+        + 1.0 * einsum("ki,jkab->ijab", op.oo, t2).antisymmetrise(0, 1)
     )
-    return AmplitudeVector(ph=f1, pphh=f2)
+    return ampl + AmplitudeVector(ph=f1, pphh=f2)
 
 
-def mtm_cvs_adc0(mp, dipop, intermediates):
-    return AmplitudeVector(ph=dipop.cv)
+def mtm_cvs_adc0(mp, op, intermediates):
+    f1 = op.cv if op.is_symmetric else op.vc.transpose()
+    return AmplitudeVector(ph=f1)
+
 
+def mtm_cvs_adc2(mp, op, intermediates):
+    op_vc = op.cv.transpose() if op.is_symmetric else op.vc
+    op_oc = op.co.transpose() if op.is_symmetric else op.oc
 
-def mtm_cvs_adc2(mp, dipop, intermediates):
+    ampl = mtm_cvs_adc0(mp, op, intermediates)
     f1 = (
-        + dipop.cv
-        - einsum("Ib,ba->Ia", dipop.cv, intermediates.cvs_p0.vv)
-        - einsum("Ij,ja->Ia", dipop.co, intermediates.cvs_p0.ov)
+        - 0.5 * einsum("bI,ab->Ia", op_vc, intermediates.cvs_p0.vv)
+        - 1.0 * einsum("jI,ja->Ia", op_oc, intermediates.cvs_p0.ov)
     )
-    f2 = (1 / sqrt(2)) * einsum("Ik,kjab->jIab", dipop.co, mp.t2(b.oovv))
-    return AmplitudeVector(ph=f1, pphh=f2)
+    f2 = (1 / sqrt(2)) * einsum("kI,kjab->jIab", op_oc, mp.t2(b.oovv))
+    return ampl + AmplitudeVector(ph=f1, pphh=f2)
 
 
 DISPATCH = {
@@ -84,7 +91,7 @@ def mtm_cvs_adc2(mp, dipop, intermediates):
 }
 
 
-def modified_transition_moments(method, ground_state, dipole_operator=None,
+def modified_transition_moments(method, ground_state, operator=None,
                                 intermediates=None):
     """Compute the modified transition moments (MTM) for the provided
     ADC method with reference to the passed ground state.
@@ -95,9 +102,9 @@ def modified_transition_moments(method, ground_state, dipole_operator=None,
         Provide a method at which to compute the MTMs
     ground_state : adcc.LazyMp
         The MP ground state
-    dipole_operator : adcc.OneParticleOperator or list, optional
-        Only required if different dipole operators than the standard
-        dipole operators in the MO basis should be used.
+    operator : adcc.OneParticleOperator or list, optional
+        Only required if different operators than the standard
+        electric dipole operators in the MO basis should be used.
     intermediates : adcc.Intermediates
         Intermediates from the ADC calculation to reuse
 
@@ -113,17 +120,17 @@ def modified_transition_moments(method, ground_state, dipole_operator=None,
         intermediates = Intermediates(ground_state)
 
     unpack = False
-    if dipole_operator is None:
-        dipole_operator = ground_state.reference_state.operators.electric_dipole
-    elif not isinstance(dipole_operator, list):
+    if operator is None:
+        operator = ground_state.reference_state.operators.electric_dipole
+    elif not isinstance(operator, list):
         unpack = True
-        dipole_operator = [dipole_operator]
+        operator = [operator]
     if method.name not in DISPATCH:
         raise NotImplementedError("modified_transition_moments is not "
                                   f"implemented for {method.name}.")
 
-    ret = [DISPATCH[method.name](ground_state, dipop, intermediates)
-           for dipop in dipole_operator]
+    ret = [DISPATCH[method.name](ground_state, op, intermediates)
+           for op in operator]
     if unpack:
         assert len(ret) == 1
         ret = ret[0]

adcc/adc_pp/test_modified_transition_moments.py

@@ -21,59 +21,72 @@
 ##
 ## ---------------------------------------------------------------------
 import unittest
+import itertools
 import numpy as np
 
 from .modified_transition_moments import modified_transition_moments
 
 from adcc.misc import expand_test_templates
 from adcc.testdata.cache import cache
 
-from pytest import approx
+from pytest import skip, approx
 
 # The methods to test
 basemethods = ["adc0", "adc1", "adc2"]
 methods = [m for bm in basemethods for m in [bm, "cvs-" + bm]]
 
+operator_kinds = ["electric", "magnetic"]
 
-@expand_test_templates(methods)
+
+@expand_test_templates(list(itertools.product(methods, operator_kinds)))
 class TestModifiedTransitionMoments(unittest.TestCase):
-    def base_test(self, system, method, kind):
-        state = cache.adc_states[system][method][kind]
-        ref = cache.reference_data[system][method][kind]
+    def base_test(self, system, method, kind, op_kind):
+        state = cache.adcc_states[system][method][kind]
+        ref = cache.adcc_reference_data[system][method][kind]
         n_ref = len(state.excitation_vector)
 
-        mtms = modified_transition_moments(method, state.ground_state)
+        if op_kind == "electric":
+            dips = state.reference_state.operators.electric_dipole
+            ref_tdm = ref["transition_dipole_moments"]
+        elif op_kind == "magnetic":
+            dips = state.reference_state.operators.magnetic_dipole
+            ref_tdm = ref["transition_magnetic_dipole_moments"]
+        else:
+            skip("Tests are only implemented for electric "
+                 "and magnetic dipole operators.")
+
+        mtms = modified_transition_moments(method, state.ground_state, dips)
+
         for i in range(n_ref):
-            # computing the scalar product of the eigenvector
+            # Computing the scalar product of the eigenvector
             # and the modified transition moments yields
             # the transition dipole moment (doi.org/10.1063/1.1752875)
             excivec = state.excitation_vector[i]
-            res_tdm = -1.0 * np.array([excivec @ mtms[i] for i in range(3)])
-            ref_tdm = ref["transition_dipole_moments"][i]
+            res_tdm = np.array([excivec @ mtms[i] for i in range(3)])
 
             # Test norm and actual values
             res_tdm_norm = np.sum(res_tdm * res_tdm)
-            ref_tdm_norm = np.sum(ref_tdm * ref_tdm)
+            ref_tdm_norm = np.sum(ref_tdm[i] * ref_tdm[i])
             assert res_tdm_norm == approx(ref_tdm_norm, abs=1e-8)
-            np.testing.assert_allclose(res_tdm, ref_tdm, atol=1e-8)
+            np.testing.assert_allclose(res_tdm, ref_tdm[i], atol=1e-8)
 
     #
     # General
     #
-    def template_h2o_sto3g_singlets(self, method):
-        self.base_test("h2o_sto3g", method, "singlet")
+    def template_h2o_sto3g_singlets(self, method, op_kind):
+        self.base_test("h2o_sto3g", method, "singlet", op_kind)
 
-    def template_h2o_def2tzvp_singlets(self, method):
-        self.base_test("h2o_def2tzvp", method, "singlet")
+    def template_h2o_def2tzvp_singlets(self, method, op_kind):
+        self.base_test("h2o_def2tzvp", method, "singlet", op_kind)
 
-    def template_h2o_sto3g_triplets(self, method):
-        self.base_test("h2o_sto3g", method, "triplet")
+    def template_h2o_sto3g_triplets(self, method, op_kind):
+        self.base_test("h2o_sto3g", method, "triplet", op_kind)
 
-    def template_h2o_def2tzvp_triplets(self, method):
-        self.base_test("h2o_def2tzvp", method, "triplet")
+    def template_h2o_def2tzvp_triplets(self, method, op_kind):
+        self.base_test("h2o_def2tzvp", method, "triplet", op_kind)
 
-    def template_cn_sto3g(self, method):
-        self.base_test("cn_sto3g", method, "state")
+    def template_cn_sto3g(self, method, op_kind):
+        self.base_test("cn_sto3g", method, "any", op_kind)
 
-    def template_cn_ccpvdz(self, method):
-        self.base_test("cn_ccpvdz", method, "state")
+    def template_cn_ccpvdz(self, method, op_kind):
+        self.base_test("cn_ccpvdz", method, "any", op_kind)

adcc/adc_pp/transition_dm.py

@@ -69,7 +69,7 @@ def tdm_cvs_adc2(mp, amplitude, intermediates):
     )
 
     # cvs_adc2_dp0_vc
-    dm.vc = u1.transpose() - einsum("ab,Ib->aI", p0.vv, u1)
+    dm.vc -= 0.5 * einsum("ab,Ib->aI", p0.vv, u1)
     return dm
 
 
@@ -86,7 +86,7 @@ def tdm_adc2(mp, amplitude, intermediates):
     # Compute ADC(2) tdm
     dm.oo = (  # adc2_dp0_oo
         - einsum("ia,ja->ij", p0.ov, u1)
-        - einsum("ikab,jkab->ij", u2, t2)
+        - einsum("ikab,jkab->ji", u2, t2)
     )
     dm.vv = (  # adc2_dp0_vv
         + einsum("ia,ib->ab", u1, p0.ov)