diff --git a/README.md b/README.md
index e75af6a..ac9820b 100644
--- a/README.md
+++ b/README.md
@@ -3,15 +3,30 @@ ResponseFun
 
 Fun with Response Functions in the Algebraic Diagrammatic Construction Framework.
 
+
+### Installation
+
+You can install `responsefun` with Conda using the conda-forge channel:
+
+```bash
+conda install responsefun -c conda-forge
+```
+
+
 ### Development
 
 A suitable `conda` environment for development can be created from the `ci_env.yml` file.
+After cloning the repository and navigating to the responsefun directory, the `responsefun` python package can be installed with the following command:
+`pip install -e .`
 
 Tests can be run with `pytest responsefun`, or `pytest --pyargs responsefun` if the package is installed.
-To exclude the slowest test, run `pytest -k "not slow" responsefun`.
+To exclude the slowest test, run `pytest -m "not slow" responsefun`.
 
 Code style is enforced through `black` (formatting), `isort` (sorting import statements), and `ruff` (linting).
 
+### Citation
+
+If you use `responsefun`, please cite [our article in JCTC.](https://doi.org/10.1021/acs.jctc.3c00456)
 ### Copyright
 
 Copyright (c) 2023, The `responsefun` Developers
diff --git a/ci_env.yml b/ci_env.yml
index 4b72a61..22d198f 100644
--- a/ci_env.yml
+++ b/ci_env.yml
@@ -3,10 +3,10 @@ channels:
   - conda-forge
 dependencies:
   - libblas =*=*mkl
-  - adcc >=0.15.16,<0.16.0
+  - adcc >=0.16.0
   - respondo >=0.0.5
   - numpy >=1.14
   - isort
   - ruff
   - pip:
-    - pyscf
\ No newline at end of file
+    - pyscf
diff --git a/examples/alpha.py b/examples/alpha.py
index 72451c5..6f60326 100644
--- a/examples/alpha.py
+++ b/examples/alpha.py
@@ -6,7 +6,7 @@
 from pyscf import gto, scf
 
 from responsefun import evaluate_property_isr, TransitionMoment
-from responsefun.symbols_and_labels import O, gamma, n, op_a, op_b, w, w_n
+from responsefun.symbols_and_labels import O, gamma, n, mu_a, mu_b, w, w_n
 
 # run SCF in PySCF
 mol = gto.M(
@@ -26,8 +26,8 @@
 
 # define symbolic SOS expression
 alpha_sos_expr = (
-        TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w - 1j*gamma)
-        + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, O) / (w_n + w + 1j*gamma)
+        TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w - 1j*gamma)
+        + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n + w + 1j*gamma)
 )
 # compute polarizability
 alpha_tens = evaluate_property_isr(
diff --git a/examples/beta_complex.py b/examples/beta_complex.py
index 2442906..65c240c 100644
--- a/examples/beta_complex.py
+++ b/examples/beta_complex.py
@@ -9,9 +9,9 @@
     O,
     gamma,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     p,
     w_1,
     w_2,
@@ -38,19 +38,16 @@
 
 # compute the complex beta tensor
 beta_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n - w_o - 1j*gamma) * (w_p - w_2 - 1j*gamma))
-    + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n + w_1 + 1j*gamma) * (w_p - w_2 - 1j*gamma))
-    + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_c, p, shifted=True)
-    * TransitionMoment(p, op_a, O) / ((w_n + w_1 + 1j*gamma) * (w_p + w_o + 1j*gamma))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n - w_o - 1j*gamma) * (w_p - w_2 - 1j*gamma))
+    + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n + w_1 + 1j*gamma) * (w_p - w_2 - 1j*gamma))
+    + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_c, p, shifted=True)
+    * TransitionMoment(p, mu_a, O) / ((w_n + w_1 + 1j*gamma) * (w_p + w_o + 1j*gamma))
 )
-# the minus sign is needed, because the negative charge is not yet included
-# in the operator definitions
-# TODO: remove minus after adc-connect/adcc#190 is merged
-beta_tens = -1.0 * evaluate_property_isr(
+beta_tens = evaluate_property_isr(
     state, beta_term, [n, p],
-    perm_pairs=[(op_b, w_1), (op_c, w_2)], excluded_states=O,
+    perm_pairs=[(mu_b, w_1), (mu_c, w_2)], excluded_states=O,
     freqs_in=[(w_1, 0.5), (w_2, 0.5)], freqs_out=(w_o, w_1+w_2),
     damping=0.01, conv_tol=1e-4,
 )
diff --git a/examples/beta_sos.py b/examples/beta_sos.py
index 65da119..799856b 100644
--- a/examples/beta_sos.py
+++ b/examples/beta_sos.py
@@ -6,9 +6,9 @@
 from responsefun.symbols_and_labels import (
     O,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     p,
     w_1,
     w_2,
@@ -18,8 +18,8 @@
 )
 
 beta_sos_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n - w_o) * (w_p - w_2))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n - w_o) * (w_p - w_2))
 )
 
 beta_sos = SumOverStates(
@@ -27,7 +27,7 @@
     [n, p],  # indices of summation
     freqs_in=[w_1, w_2],  # frequencies of incident photons
     freqs_out=w_o,  # frequency of resulting photon
-    perm_pairs=[(op_a, -w_o), (op_b, w_1), (op_c, w_2)],  # tuples to be permuted
+    perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],  # tuples to be permuted
     excluded_states=O  # states excluded from the summations
 )
 
diff --git a/examples/esp.py b/examples/esp.py
index 11abe22..000739c 100644
--- a/examples/esp.py
+++ b/examples/esp.py
@@ -6,7 +6,7 @@
 from pyscf import gto, scf
 
 from responsefun import evaluate_property_isr, TransitionMoment
-from responsefun.symbols_and_labels import f, gamma, n, op_a, op_b, w, w_f, w_n
+from responsefun.symbols_and_labels import f, gamma, n, mu_a, mu_b, w, w_f, w_n
 
 # run SCF in PySCF
 mol = gto.M(
@@ -27,8 +27,8 @@
 
 # compute esp tensor
 sos_expr = (
-    TransitionMoment(f, op_a, n) * TransitionMoment(n, op_b, f) / (w_n - w_f - w - 1j*gamma)
-    + TransitionMoment(f, op_b, n) * TransitionMoment(n, op_a, f) / (w_n - w_f + w + 1j*gamma)
+    TransitionMoment(f, mu_a, n) * TransitionMoment(n, mu_b, f) / (w_n - w_f - w - 1j*gamma)
+    + TransitionMoment(f, mu_b, n) * TransitionMoment(n, mu_a, f) / (w_n - w_f + w + 1j*gamma)
 )
 tens = evaluate_property_isr(
     state,  # ExcitedStates object returned by the adcc calculation
diff --git a/examples/mcd.py b/examples/mcd.py
index bf84122..7509dd4 100644
--- a/examples/mcd.py
+++ b/examples/mcd.py
@@ -8,7 +8,7 @@
 
 from responsefun import evaluate_property_isr, TransitionMoment
 from responsefun.misc import epsilon
-from responsefun.symbols_and_labels import O, j, k, op_a, op_b, opm_c, w_j, w_k
+from responsefun.symbols_and_labels import O, j, k, mu_a, mu_b, m_c, w_j, w_k
 
 # run SCF in PySCF
 mol = gto.M(
@@ -28,11 +28,11 @@
 
 # define symbolic SOS expressions
 mcd_sos_expr1 = (
-        TransitionMoment(O, opm_c, k) * TransitionMoment(k, op_b, j, shifted=True)
-        * TransitionMoment(j, op_a, O) / w_k
+        TransitionMoment(O, m_c, k) * TransitionMoment(k, mu_b, j, shifted=True)
+        * TransitionMoment(j, mu_a, O) / w_k
 )
 mcd_sos_expr2 = (
-        TransitionMoment(O, op_b, k) * TransitionMoment(k, opm_c, j) * TransitionMoment(j, op_a, O)
+        TransitionMoment(O, mu_b, k) * TransitionMoment(k, m_c, j) * TransitionMoment(j, mu_a, O)
         / (w_k - w_j)
 )
 # compute MCD B term for the first excited state
@@ -48,8 +48,5 @@
     conv_tol=1e-4,
 )
 
-# the minus sign is needed, because the negative charge is not yet included
-# in the operator definitions
-# TODO: remove minus after adc-connect/adcc#190 is merged
-bterm = -1.0 * np.einsum("abc,abc->", epsilon, mcd_tens1 + mcd_tens2)
+bterm = np.einsum("abc,abc->", epsilon, mcd_tens1 + mcd_tens2)
 print(f"The MCD Faraday B term for excited state {final_state} is {bterm:.2f} (a.u.).")
\ No newline at end of file
diff --git a/examples/rixs.py b/examples/rixs.py
index 036c89c..a5230d4 100644
--- a/examples/rixs.py
+++ b/examples/rixs.py
@@ -12,8 +12,8 @@
     f,
     gamma,
     n,
-    op_a,
-    op_b,
+    mu_a,
+    mu_b,
     w,
     w_f,
     w_n,
@@ -38,7 +38,7 @@
 state = adcc.adc2(scfres, n_singlets=5)
 
 # compute RIXS tensor within the rotating-wave approximation
-rixs_sos_rwa = TransitionMoment(f, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w - 1j * gamma)
+rixs_sos_rwa = TransitionMoment(f, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w - 1j * gamma)
 rixs_rwa = evaluate_property_isr(
     state,
     rixs_sos_rwa,
@@ -55,7 +55,7 @@
     state,
     rixs_sos_rwa,
     [n],
-    perm_pairs=[(op_a, w + 1j * gamma), (op_b, -w_prime - 1j * gamma)],
+    perm_pairs=[(mu_a, w + 1j * gamma), (mu_b, -w_prime - 1j * gamma)],
     freqs_in=(w, ev2au(534.74)),
     freqs_out=(w_prime, w-w_f),
     damping=ev2au(0.124),
diff --git a/examples/second_order_nlo.py b/examples/second_order_nlo.py
index f9425d0..7286d2c 100644
--- a/examples/second_order_nlo.py
+++ b/examples/second_order_nlo.py
@@ -11,9 +11,9 @@
 from responsefun.symbols_and_labels import (
     O,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     p,
     w_1,
     w_2,
@@ -56,20 +56,17 @@ def compute_beta_parallel(beta_tens, dip_mom):
 
 # compute the first hyperpolarizability tensor
 beta_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n - w_o) * (w_p - w_2))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n - w_o) * (w_p - w_2))
 )
 processes = {
     "static": (0.0, 0.0), "OR": (w_ruby, -w_ruby),
     "EOPE": (w_ruby, 0.0), "SHG": (w_ruby, w_ruby)
 }
 for process, freqs in processes.items():
-    # the minus sign is needed, because the negative charge is not yet included
-    # in the operator definitions
-    # TODO: remove minus after adc-connect/adcc#190 is merged
-    beta_tens = -1.0 * evaluate_property_isr(
+    beta_tens = evaluate_property_isr(
         state, beta_term, [n, p],
-        perm_pairs=[(op_a, -w_o), (op_b, w_1), (op_c, w_2)],
+        perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],
         freqs_in=[(w_1, freqs[0]), (w_2, freqs[1])],
         freqs_out=(w_o, w_1+w_2),
         excluded_states=O,
diff --git a/examples/third_order_nlo.py b/examples/third_order_nlo.py
index 0105fa4..eab7c85 100644
--- a/examples/third_order_nlo.py
+++ b/examples/third_order_nlo.py
@@ -12,10 +12,10 @@
     O,
     m,
     n,
-    op_a,
-    op_b,
-    op_c,
-    op_d,
+    mu_a,
+    mu_b,
+    mu_c,
+    mu_d,
     p,
     w_1,
     w_2,
@@ -55,16 +55,16 @@ def compute_gamma_average(gamma_tens):
 state = adcc.adc2(scfres, n_singlets=5)
 # compute the second hyperpolarizability tensor
 gamma_term_I = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, m, shifted=True)
-    * TransitionMoment(m, op_c, p, shifted=True) * TransitionMoment(p, op_d, O)
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, m, shifted=True)
+    * TransitionMoment(m, mu_c, p, shifted=True) * TransitionMoment(p, mu_d, O)
     / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
 )
 gamma_term_II = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, O)
-    * TransitionMoment(O, op_c, m) * TransitionMoment(m, op_d, O)
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O)
+    * TransitionMoment(O, mu_c, m) * TransitionMoment(m, mu_d, O)
     / ((w_n - w_o) * (w_m - w_3) * (w_m + w_2))
 )
-perm_pairs = [(op_a, -w_o), (op_b, w_1), (op_c, w_2), (op_d, w_3)]
+perm_pairs = [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2), (mu_d, w_3)]
 processes = {
     "static": (0.0, 0.0, 0.0), "dcOR": (w_ruby, -w_ruby, 0.0),
     "EOKE": (w_ruby, 0.0, 0.0), "IDRI": (w_ruby, -w_ruby, w_ruby),
diff --git a/examples/threepa.py b/examples/threepa.py
index 1697c5b..e6d7730 100644
--- a/examples/threepa.py
+++ b/examples/threepa.py
@@ -11,9 +11,9 @@
     f,
     m,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     w_1,
     w_2,
     w_3,
@@ -45,18 +45,15 @@ def threepa_average(tens):
 print(state.describe())
 
 threepa_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, m)
-    * TransitionMoment(m, op_c, f) / ((w_n - w_1) * (w_m - w_1 - w_2))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, m)
+    * TransitionMoment(m, mu_c, f) / ((w_n - w_1) * (w_m - w_1 - w_2))
 )
 
 for es in range(5):
     print(f"===== State {es} ===== ")
-    # the minus sign is needed, because the negative charge is not yet included
-    # in the operator definitions
-    # TODO: remove minus after adc-connect/adcc#190 is merged
-    threepa_tens = -1.0 * evaluate_property_isr(
+    threepa_tens = evaluate_property_isr(
         state, threepa_term, [n, m],
-        perm_pairs=[(op_a, w_1), (op_b, w_2), (op_c, w_3)],
+        perm_pairs=[(mu_a, w_1), (mu_b, w_2), (mu_c, w_3)],
         freqs_in=[(w_1, w_f/3), (w_2, w_f/3), (w_3, w_f/3)],
         excited_state=es, conv_tol=1e-5,
     )
diff --git a/examples/tpa.py b/examples/tpa.py
index 4ecde7f..b2948c1 100644
--- a/examples/tpa.py
+++ b/examples/tpa.py
@@ -10,8 +10,8 @@
     O,
     f,
     n,
-    op_a,
-    op_b,
+    mu_a,
+    mu_b,
     w_1,
     w_2,
     w_n,
@@ -42,14 +42,14 @@ def tpa_average(tens):
 print(state.describe())
 
 tpa_term = (
-    TransitionMoment(f, op_b, n) * TransitionMoment(n, op_a, O) / (w_n - w_1)
+    TransitionMoment(f, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n - w_1)
 )
 
 for es in range(1):
     print(f"===== State {es} ===== ")
     tpa_tens = evaluate_property_isr(
         state, tpa_term, [n],
-        perm_pairs=[(op_a, w_1), (op_b, w_2)],
+        perm_pairs=[(mu_a, w_1), (mu_b, w_2)],
         freqs_in=[(w_1, w_f/2), (w_2, w_f/2)],
         excited_state=es, conv_tol=1e-4,
     )
diff --git a/pyproject.toml b/pyproject.toml
index f44a78f..36de855 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,10 +21,10 @@ classifiers = [
 requires-python = ">=3.7"
 # Declare any run-time dependencies that should be installed with the package.
 dependencies = [
-    "adcc>=0.15.16,<0.16.0",
+    "adcc>=0.16.0",
     "respondo>=0.0.5",
     "numpy>=1.14",
-    "sympy",
+    "sympy<=1.13",
     "tqdm",
     "anytree",
 ]
diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py
index 5d1de6e..4b7b099 100644
--- a/responsefun/AdccProperties.py
+++ b/responsefun/AdccProperties.py
@@ -15,7 +15,6 @@
 #  You should have received a copy of the GNU Lesser General Public License
 #  along with responsefun. If not, see <http:www.gnu.org/licenses/>.
 #
-
 import warnings
 from abc import ABC, abstractmethod, abstractproperty
 from dataclasses import dataclass
@@ -64,6 +63,13 @@ class Operator:
         dim=1,
         is_imag=False,
     ),
+    Operator(
+        name="electric_dipole_velocity",
+        symbol="mu_p",
+        symmetry=Symmetry.ANTIHERMITIAN,
+        dim=1,
+        is_imag=True,
+    ),
     Operator(
         name="magnetic_dipole",
         symbol="m",
@@ -71,6 +77,20 @@ class Operator:
         dim=1,
         is_imag=True,
     ),
+    Operator(
+        name="electric_quadrupole",
+        symbol="q",
+        symmetry=Symmetry.HERMITIAN,
+        dim=2,
+        is_imag=False,
+    ),
+    Operator(
+        name="electric_quadrupole_velocity",
+        symbol="q_p",
+        symmetry=Symmetry.ANTIHERMITIAN,
+        dim=2,
+        is_imag=True,
+    ),
     Operator(
         name="diamagnetic_magnetizability",
         symbol="xi",
@@ -142,12 +162,43 @@ def compute_state_to_state_transition_moments(state, integrals, initial_state=No
     return np.squeeze(s2s_tm)
 
 
+def compute_ground_state_moment(state, integrals, pm_level,
+                                nuclear_contribution=None) -> np.ndarray:
+    """
+    Computes the ground state moments. 
+    """
+    # convert integrals to np.array (tuples are not allowed as indices for lists)
+    integrals = np.array(integrals)
+    op_shape = np.shape(integrals)
+    iterables = [list(range(shape)) for shape in op_shape]
+    components = list(product(*iterables))
+    ref_mom = np.zeros(op_shape)
+    for c in components:
+        ref_mom[c] = product_trace(integrals[c], state.ground_state.reference_state.density)
+
+    if pm_level in [0, 1]:
+        gs_mom = ref_mom
+    elif pm_level == 2:
+        mp2corr = np.zeros(op_shape)
+        for c in components:
+            mp2corr[c] = product_trace(integrals[c], state.ground_state.mp2_diffdm)
+        gs_mom = ref_mom + mp2corr
+    else:
+        raise NotImplementedError(
+            "Only ground state moments for level 1 and 2 are implemented."
+        )
+    if nuclear_contribution is not None:
+        gs_mom = gs_mom + nuclear_contribution
+
+    return gs_mom
+
+    
 class AdccProperties(ABC):
     """Abstract base class encompassing all properties that can be obtained
     from adcc for a given operator."""
 
     def __init__(self, state: Union[adcc.ExcitedStates, MockExcitedStates],
-                 gauge_origin: Union[str, tuple[float, float, float], None] = None):
+                 gauge_origin: Union[str, tuple[float, float, float]] = "origin"):
         self._state = state
         self._state_size = len(state.excitation_energy_uncorrected)
         self._property_method = self._state.property_method
@@ -158,6 +209,11 @@ def __init__(self, state: Union[adcc.ExcitedStates, MockExcitedStates],
 
         self._gauge_origin = gauge_origin
 
+        if isinstance(self._state, MockExcitedStates) and \
+            self._gauge_origin != "origin":
+            raise NotImplementedError("MockExcitedStates class is only available "
+                                      "with gauge origin selected as origin.")
+
         # to make things faster if not all state-to-state transition moments are needed
         # but only from or to a specific state
         self._s2s_tm_i = np.empty((self._state_size), dtype=object)
@@ -195,7 +251,10 @@ def revert_transition_moment(self, moment: Any) -> Any:
 
     @cached_property
     def transition_moment(self) -> np.ndarray:
-        return self._transition_moment()
+        if isinstance(self._state, MockExcitedStates):
+            return self._transition_moment_mock()
+        else:
+            return compute_transition_moments(self._state, self.integrals)
 
     @property
     def transition_moment_reverse(self) -> np.ndarray:
@@ -203,7 +262,10 @@ def transition_moment_reverse(self) -> np.ndarray:
 
     @cached_property
     def state_to_state_transition_moment(self) -> np.ndarray:
-        return self._state_to_state_transition_moment()
+        if isinstance(self._state, MockExcitedStates):
+            return self._state_to_state_transition_moment_mock()
+        else:
+            return compute_state_to_state_transition_moments(self._state, self.integrals)
     
     def s2s_tm_view(self, initial_state=None, final_state=None):
         if initial_state is None and final_state is None:
@@ -233,31 +295,31 @@ def s2s_tm_view(self, initial_state=None, final_state=None):
             return s2s_tm
 
     @abstractmethod
-    def _transition_moment(self) -> np.ndarray:
+    def _transition_moment_mock(self) -> np.ndarray:
         pass
 
     @abstractmethod
-    def _state_to_state_transition_moment(self) -> np.ndarray:
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
         pass
 
     def modified_transition_moments(
-        self, comp: Union[int, None] = None
-    ) -> Union[adcc.AmplitudeVector, list[adcc.AmplitudeVector]]:
+        self, comp: Union[int, tuple[int], None] = None
+    ) -> Union[adcc.AmplitudeVector, np.ndarray[adcc.AmplitudeVector]]:
         if comp is None:
-            op = self.integrals
+            mtms = np.array([(modified_transition_moments(self._property_method,
+                             self._state.ground_state, op_ints)) for op_ints in self.integrals])
         else:
-            op = self.integrals[comp]
-        mtms = modified_transition_moments(
-            self._property_method, self._state.ground_state, op
-        )
+            op = np.array(self.integrals)[comp]
+            mtms = modified_transition_moments(self._property_method,
+                                               self._state.ground_state, op)
         return mtms
 
     def modified_transition_moments_reverse(
-        self, comp: Union[int, None] = None
-    ) -> Union[adcc.AmplitudeVector, list[adcc.AmplitudeVector]]:
+        self, comp: Union[int, tuple[int], None] = None
+    ) -> Union[adcc.AmplitudeVector, np.ndarray[adcc.AmplitudeVector]]:
         return self.revert_transition_moment(self.modified_transition_moments(comp))
 
-    def isr_matrix(self, comp: Union[int, None] = None) -> adcc.IsrMatrix:
+    def isr_matrix(self, comp: Union[int, tuple[int], None] = None) -> adcc.IsrMatrix:
         if comp is None:
             op = self.integrals
         else:
@@ -268,7 +330,7 @@ def transition_polarizability(
         self,
         to_vec: Union[adcc.AmplitudeVector, RV],
         from_vec: Union[adcc.AmplitudeVector, RV],
-        comp: Union[int, None] = None
+        comp: Union[int, tuple[int], None] = None
         ) -> np.ndarray:
         if comp is None:
             op = self.integrals
@@ -297,12 +359,20 @@ def transition_polarizability(
 def build_adcc_properties(
     state: Union[adcc.ExcitedStates, MockExcitedStates],
     op_type: str,
-    gauge_origin: Union[str, tuple[float, float, float], None] = None
+    gauge_origin: Union[str, tuple[float, float, float]] = "origin"
 ) -> AdccProperties:
     if op_type == "electric_dipole":
         return ElectricDipole(state, gauge_origin)
+    elif op_type == "electric_dipole_velocity":
+        return ElectricDipoleVelocity(state, gauge_origin)
     elif op_type == "magnetic_dipole":
         return MagneticDipole(state, gauge_origin)
+    elif op_type == "electric_quadrupole":
+        return ElectricQuadrupole(state, gauge_origin)
+    elif op_type == "electric_quadrupole_velocity":
+        return ElectricQuadrupoleVelocity(state, gauge_origin)
+    elif op_type == "diamagnetic_magnetizability":
+        return DiamagneticMagnetizability(state, gauge_origin)
     else:
         raise NotImplementedError
 
@@ -314,7 +384,10 @@ def _operator(self) -> Operator:
 
     @property
     def integrals(self) -> list[adcc.OneParticleOperator]:
-        return self._state.reference_state.operators.electric_dipole
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
+        else:
+            return self._state.reference_state.operators.electric_dipole
 
     @property
     def gs_moment(self) -> np.ndarray:
@@ -323,14 +396,11 @@ def gs_moment(self) -> np.ndarray:
         else:
             return self._state.ground_state.dipole_moment(self._pm_level)
 
-    def _transition_moment(self) -> np.ndarray:
+    def _transition_moment_mock(self) -> np.ndarray:
         return self._state.transition_dipole_moment
 
-    def _state_to_state_transition_moment(self) -> np.ndarray:
-        if isinstance(self._state, MockExcitedStates):
-            return self._state.transition_dipole_moment_s2s
-        else:
-            return compute_state_to_state_transition_moments(self._state, self.integrals)
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
+        return self._state.transition_dipole_moment_s2s
 
 
 class MagneticDipole(AdccProperties):
@@ -340,33 +410,125 @@ def _operator(self) -> Operator:
 
     @property
     def integrals(self) -> list[adcc.OneParticleOperator]:
-        return self._state.reference_state.operators.magnetic_dipole
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
+        else:
+            return self._state.reference_state.operators.magnetic_dipole(self._gauge_origin)
 
     @property
     def gs_moment(self) -> np.ndarray:
-        # the minus sign is needed, because the negative charge is not yet included
-        # in the operator definitions
-        # TODO: remove minus after adc-connect/adcc#190 is merged
-        ref_dipmom = -1.0 * np.array(
-            [product_trace(dip, self._state.ground_state.reference_state.density)
-             for dip in self.integrals]
-        )
-        if self._pm_level in [0, 1]:
-            return ref_dipmom
-        elif self._pm_level == 2:
-            mp2corr = -1.0 * np.array([product_trace(dip, self._state.ground_state.mp2_diffdm)
-                                       for dip in self.integrals])
-            return ref_dipmom + mp2corr
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
         else:
-            raise NotImplementedError(
-                "Only magnetic dipole moments for level 1 and 2 are implemented."
-            )
+            return compute_ground_state_moment(self._state, self.integrals, self._pm_level)
 
-    def _transition_moment(self) -> np.ndarray:
+    def _transition_moment_mock(self) -> np.ndarray:
         return self._state.transition_magnetic_dipole_moment
 
-    def _state_to_state_transition_moment(self) -> np.ndarray:
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
+        return self._state.transition_magnetic_moment_s2s
+
+
+class ElectricDipoleVelocity(AdccProperties):
+    @property
+    def _operator(self) -> Operator:
+        return get_operator_by_name("electric_dipole_velocity")
+
+    @property
+    def integrals(self) -> list[adcc.OneParticleOperator]:
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
+        else:
+            return self._state.reference_state.operators.electric_dipole_velocity
+
+    @property
+    def gs_moment(self) -> np.ndarray:
+        raise NotImplementedError("Ground-state moments not implemented for electric dipole "
+                                  "operator in velocity gauge.")
+
+    def _transition_moment_mock(self) -> np.ndarray:
+        return self._state.transition_dipole_moment_velocity
+
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
+        raise NotImplementedError
+
+
+class ElectricQuadrupole(AdccProperties):
+    @property
+    def _operator(self) -> Operator:
+        return get_operator_by_name("electric_quadrupole")
+
+    @property
+    def integrals(self) -> list[list[adcc.OneParticleOperator]]:
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
+        else:
+            return self._state.reference_state.operators.electric_quadrupole(self._gauge_origin)
+
+    @property
+    def gs_moment(self) -> np.ndarray:
+        nuclear_contribution = np.zeros((3,3))
+        nuc_electric_quadrupole = self._state.reference_state.nuclear_quadrupole(self._gauge_origin)
+        nuclear_contribution[np.triu_indices(3)] = nuc_electric_quadrupole
+        nuclear_contribution = nuclear_contribution + nuclear_contribution.T \
+            - np.diag(nuclear_contribution)
+        return compute_ground_state_moment(self._state, self.integrals, self._pm_level,
+                                           nuclear_contribution=nuclear_contribution)
+
+    def _transition_moment_mock(self) -> np.ndarray:
+        return self._state.transition_quadrupole_moment
+
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
+        raise NotImplementedError
+
+
+class ElectricQuadrupoleVelocity(AdccProperties):
+    @property
+    def _operator(self) -> Operator:
+        return get_operator_by_name("electric_quadrupole_velocity")
+
+    @property
+    def integrals(self) -> list[list[adcc.OneParticleOperator]]:
         if isinstance(self._state, MockExcitedStates):
-            return self._state.transition_magnetic_moment_s2s
+            raise NotImplementedError
         else:
-            return compute_state_to_state_transition_moments(self._state, self.integrals)
\ No newline at end of file
+            return self._state.reference_state.operators.\
+                electric_quadrupole_velocity(self._gauge_origin)
+
+    @property
+    def gs_moment(self) -> np.ndarray:
+        raise NotImplementedError("Ground-state moments not implemented for electric "
+                                  "quadrupole operator in velocity gauge.")
+
+    def _transition_moment_mock(self) -> np.ndarray:
+        return self._state.transition_quadrupole_moment_velocity
+
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
+        raise NotImplementedError
+
+
+class DiamagneticMagnetizability(AdccProperties):
+    @property
+    def _operator(self) -> Operator:
+        return get_operator_by_name("diamagentic_magnetizability")
+
+    @property
+    def integrals(self) -> list[list[adcc.OneParticleOperator]]:
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
+        else:
+            return self._state.reference_state.operators.\
+                diamagnetic_magnetizability(self._gauge_origin)
+
+    @property
+    def gs_moment(self) -> np.ndarray:
+        if isinstance(self._state, MockExcitedStates):
+            raise NotImplementedError
+        else:
+            return compute_ground_state_moment(self._state, self.integrals, self._pm_level)
+
+    def _transition_moment_mock(self) -> np.ndarray:
+        raise NotImplementedError
+
+    def _state_to_state_transition_moment_mock(self) -> np.ndarray:
+        raise NotImplementedError
\ No newline at end of file
diff --git a/responsefun/SumOverStates.py b/responsefun/SumOverStates.py
index 3e9a127..28801e1 100644
--- a/responsefun/SumOverStates.py
+++ b/responsefun/SumOverStates.py
@@ -180,7 +180,7 @@ def _build_sos_via_permutation(term, perm_pairs):
         List of (op, freq) pairs whose permutation yields the full SOS expression;
         (op, freq): (<class 'responsefun.operators.OneParticleOperator'>,
         <class 'sympy.core.symbol.Symbol'>),
-        e.g., [(op_a, -w_o), (op_b, w_1), (op_c, w_2)].
+        e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)].
 
     Returns
     ----------
@@ -325,7 +325,7 @@ def __init__(
             List of (op, freq) pairs whose permutation yields the full SOS expression;
             (op, freq): (<class 'responsefun.operators.OneParticleOperator'>,
             <class 'sympy.core.symbol.Symbol'>),
-            e.g., [(op_a, -w_o), (op_b, w_1), (op_c, w_2)].
+            e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)].
 
         excluded_states: list of <class 'sympy.core.symbol.Symbol'> or int, optional
             List of states that are excluded from the summation.
@@ -395,12 +395,12 @@ def __init__(
         validate_summation_indices(self.expr, self.summation_indices)
 
         self._operators, self._operators_unshifted = extract_operators_from_sos(self.expr)
-        self._components = {op.comp for op in self._operators}
+        self._components = {op_char for op in self._operators for op_char in op.comp}
         self._order = len(self._components)
         if self._components.difference(ABC[: self._order]):
             raise ValueError(
                 f"It is important that the Cartesian components of an order {self._order} tensor "
-                f"be specified as {ABC[:self._order]}."
+                f"be specified as {ABC[:self._order]} and not {self._components}."
             )
         self._initial_state, self._final_state, self._excited_state = \
             extract_initial_final_excited_from_sos(self.expr, self.summation_indices)
diff --git a/responsefun/evaluate_property.py b/responsefun/evaluate_property.py
index 6d2ac0f..3067a74 100644
--- a/responsefun/evaluate_property.py
+++ b/responsefun/evaluate_property.py
@@ -476,6 +476,7 @@ def evaluate_property_isr(
     omegas=None,
     gamma_val=None,
     final_state=None,
+    gauge_origin="origin",
     **solver_args,
 ):
     """Compute a molecular property with the ADC/ISR approach from its SOS expression.
@@ -497,7 +498,7 @@ def evaluate_property_isr(
         List of (op, freq) pairs whose permutation yields the full SOS expression;
         (op, freq): (<class 'responsefun.operators.OneParticleOperator'>,
         <class 'sympy.core.symbol.Symbol'>),
-        e.g., [(op_a, -w_o), (op_b, w_1), (op_c, w_2)].
+        e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)].
 
     excluded_states: list of <class 'sympy.core.symbol.Symbol'> or int, optional
         List of states that are excluded from the summation.
@@ -622,7 +623,7 @@ def projection(X, bl=None):
     # store adcc properties for the required operators in a dict
     adcc_prop = {}
     for op_type in sos.operator_types:
-        adcc_prop[op_type] = build_adcc_properties(state, op_type)
+        adcc_prop[op_type] = build_adcc_properties(state, op_type, gauge_origin=gauge_origin)
 
     rvecs_dict_tot, rvecs_solution, rvecs_mapping = determine_rvecs(
         rvecs_dict_list, input_subs, adcc_prop, state, projection, **solver_args
@@ -776,6 +777,7 @@ def evaluate_property_sos(
     omegas=None,
     gamma_val=None,
     final_state=None,
+    gauge_origin="origin",
 ):
     """Compute a molecular property from its SOS expression.
 
@@ -796,7 +798,7 @@ def evaluate_property_sos(
         List of (op, freq) pairs whose permutation yields the full SOS expression;
         (op, freq): (<class 'responsefun.operators.OneParticleOperator'>,
         <class 'sympy.core.symbol.Symbol'>),
-        e.g., [(op_a, -w_o), (op_b, w_1), (op_c, w_2)].
+        e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)].
 
     excluded_states: list of <class 'sympy.core.symbol.Symbol'> or int, optional
         List of states that are excluded from the summation.
@@ -944,7 +946,7 @@ def evaluate_property_sos(
     # store adcc properties for the required operators in a dict
     adcc_prop = {}
     for op_type in sos.operator_types:
-        adcc_prop[op_type] = build_adcc_properties(state, op_type)
+        adcc_prop[op_type] = build_adcc_properties(state, op_type, gauge_origin=gauge_origin)
 
     print(f"Summing over {len(state.excitation_energy_uncorrected)} excited states ...")
     for term_dict in tqdm(term_list):
@@ -1030,6 +1032,7 @@ def evaluate_property_sos_fast(
     omegas=None,
     gamma_val=None,
     final_state=None,
+    gauge_origin="origin",
 ):
     """Compute a molecular property from its SOS expression using the Einstein summation convention.
 
@@ -1050,7 +1053,7 @@ def evaluate_property_sos_fast(
         List of (op, freq) pairs whose permutation yields the full SOS expression;
         (op, freq): (<class 'responsefun.operators.OneParticleOperator'>,
         <class 'sympy.core.symbol.Symbol'>),
-        e.g., [(op_a, -w_o), (op_b, w_1), (op_c, w_2)].
+        e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)].
 
     excluded_states: list of <class 'sympy.core.symbol.Symbol'> or int, optional
         List of states that are excluded from the summation.
@@ -1173,7 +1176,7 @@ def evaluate_property_sos_fast(
     # store adcc properties for the required operators in a dict
     adcc_prop = {}
     for op_type in sos.operator_types:
-        adcc_prop[op_type] = build_adcc_properties(state, op_type)
+        adcc_prop[op_type] = build_adcc_properties(state, op_type, gauge_origin=gauge_origin)
 
     for it, term in enumerate(term_list):
         einsum_list = []
diff --git a/responsefun/symbols_and_labels.py b/responsefun/symbols_and_labels.py
index 2e61d17..97b500f 100644
--- a/responsefun/symbols_and_labels.py
+++ b/responsefun/symbols_and_labels.py
@@ -24,15 +24,43 @@
 w_k = TransitionFrequency(k, real=True)
 
 # electric dipole operators
-op_a = OneParticleOperator("A", "electric_dipole", False)
-op_b = OneParticleOperator("B", "electric_dipole", False)
-op_c = OneParticleOperator("C", "electric_dipole", False)
-op_d = OneParticleOperator("D", "electric_dipole", False)
-op_e = OneParticleOperator("E", "electric_dipole", False)
+mu_a = OneParticleOperator("A", "electric_dipole", False)
+mu_b = OneParticleOperator("B", "electric_dipole", False)
+mu_c = OneParticleOperator("C", "electric_dipole", False)
+mu_d = OneParticleOperator("D", "electric_dipole", False)
+mu_e = OneParticleOperator("E", "electric_dipole", False)
+
+# electric dipole operators in velocity gauge
+mup_a = OneParticleOperator("A", "electric_dipole_velocity", False)
+mup_b = OneParticleOperator("B", "electric_dipole_velocity", False)
+mup_c = OneParticleOperator("C", "electric_dipole_velocity", False)
+mup_d = OneParticleOperator("D", "electric_dipole_velocity", False)
+mup_e = OneParticleOperator("E", "electric_dipole_velocity", False)
 
 # magnetic dipole operators
-opm_a = OneParticleOperator("A", "magnetic_dipole", False)
-opm_b = OneParticleOperator("B", "magnetic_dipole", False)
-opm_c = OneParticleOperator("C", "magnetic_dipole", False)
-opm_d = OneParticleOperator("D", "magnetic_dipole", False)
-opm_e = OneParticleOperator("E", "magnetic_dipole", False)
\ No newline at end of file
+m_a = OneParticleOperator("A", "magnetic_dipole", False)
+m_b = OneParticleOperator("B", "magnetic_dipole", False)
+m_c = OneParticleOperator("C", "magnetic_dipole", False)
+m_d = OneParticleOperator("D", "magnetic_dipole", False)
+m_e = OneParticleOperator("E", "magnetic_dipole", False)
+
+# electric quadrupole operators
+q_ab = OneParticleOperator("AB", "electric_quadrupole", False)
+q_bc = OneParticleOperator("BC", "electric_quadrupole", False)
+q_cd = OneParticleOperator("CD", "electric_quadrupole", False)
+q_de = OneParticleOperator("DE", "electric_quadrupole", False)
+q_ef = OneParticleOperator("EF", "electric_quadrupole", False)
+
+# electric quadrupole operators in velocity gauge
+qp_ab = OneParticleOperator("AB", "electric_quadrupole_velocity", False)
+qp_bc = OneParticleOperator("BC", "electric_quadrupole_velocity", False)
+qp_cd = OneParticleOperator("CD", "electric_quadrupole_velocity", False)
+qp_de = OneParticleOperator("DE", "electric_quadrupole_velocity", False)
+qp_ef = OneParticleOperator("EF", "electric_quadrupole_velocity", False)
+
+# diamagnetic magnetizability operators
+xi_ab = OneParticleOperator("AB", "diamagnetic_magnetizability", False)
+xi_bc = OneParticleOperator("BC", "diamagnetic_magnetizability", False)
+xi_cd = OneParticleOperator("CD", "diamagnetic_magnetizability", False)
+xi_de = OneParticleOperator("DE", "diamagnetic_magnetizability", False)
+xi_ef = OneParticleOperator("EF", "diamagnetic_magnetizability", False)
\ No newline at end of file
diff --git a/responsefun/test_extra_terms.py b/responsefun/test_extra_terms.py
index 3150f15..3979036 100644
--- a/responsefun/test_extra_terms.py
+++ b/responsefun/test_extra_terms.py
@@ -8,9 +8,9 @@
     O,
     k,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     w_1,
     w_2,
     w_k,
@@ -34,21 +34,21 @@ def run_scf(molecule, basis, backend="pyscf"):
 SOS_expressions = {
     "beta1": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, k)
-            * TransitionMoment(k, op_c, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, k)
+            * TransitionMoment(k, mu_c, O)
             / ((w_n - w_o) * (w_k - w_2))
         ),
-        [(op_a, -w_o), (op_b, w_1), (op_c, w_2)],
+        [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],
     ),
     "beta2": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, k, shifted=True)
-            * TransitionMoment(k, op_c, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, k, shifted=True)
+            * TransitionMoment(k, mu_c, O)
             / ((w_n - w_o) * (w_k - w_2))
         ),
-        [(op_a, -w_o), (op_b, w_1), (op_c, w_2)],
+        [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],
     ),    
 }
 
diff --git a/responsefun/test_magnetic_properties.py b/responsefun/test_magnetic_properties.py
index afe0e7f..6b9ed34 100644
--- a/responsefun/test_magnetic_properties.py
+++ b/responsefun/test_magnetic_properties.py
@@ -11,17 +11,17 @@
     O,
     k,
     m,
+    m_a,
+    m_b,
+    m_c,
+    m_d,
+    m_e,
+    mu_a,
+    mu_b,
+    mu_c,
+    mu_d,
+    mu_e,
     n,
-    op_a,
-    op_b,
-    op_c,
-    op_d,
-    op_e,
-    opm_a,
-    opm_b,
-    opm_c,
-    opm_d,
-    opm_e,
     p,
     w,
     w_1,
@@ -49,55 +49,55 @@ def run_scf(molecule, basis, backend="pyscf"):
 
 SOS_alpha_like = {
     "a": (
-        TransitionMoment(O, opm_a, n) * TransitionMoment(n, op_b, O) / (w_n - w)
-        + TransitionMoment(O, op_b, n) * TransitionMoment(n, opm_a, O) / (w_n + w)
+        TransitionMoment(O, m_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w)
+        + TransitionMoment(O, mu_b, n) * TransitionMoment(n, m_a, O) / (w_n + w)
     ),
     "b": (
-        TransitionMoment(O, op_a, n) * TransitionMoment(n, opm_b, O) / (w_n - w)
-        + TransitionMoment(O, opm_b, n) * TransitionMoment(n, op_a, O) / (w_n + w)
+        TransitionMoment(O, mu_a, n) * TransitionMoment(n, m_b, O) / (w_n - w)
+        + TransitionMoment(O, m_b, n) * TransitionMoment(n, mu_a, O) / (w_n + w)
     ),
     "ab": (
-        TransitionMoment(O, opm_a, n) * TransitionMoment(n, opm_b, O) / (w_n - w)
-        + TransitionMoment(O, opm_b, n) * TransitionMoment(n, opm_a, O) / (w_n + w)
+        TransitionMoment(O, m_a, n) * TransitionMoment(n, m_b, O) / (w_n - w)
+        + TransitionMoment(O, m_b, n) * TransitionMoment(n, m_a, O) / (w_n + w)
     ),
 }
 
 
 SOS_beta_like = {
     "a": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, k)
-        * TransitionMoment(k, op_c, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, k)
+        * TransitionMoment(k, mu_c, O)
         / ((w_n - w_o) * (w_k - w_2))
     ),
     "b": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, opm_b, k)
-        * TransitionMoment(k, op_c, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, m_b, k)
+        * TransitionMoment(k, mu_c, O)
         / ((w_n - w_o) * (w_k - w_2))
     ),
     "c": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, op_b, k)
-        * TransitionMoment(k, opm_c, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, mu_b, k)
+        * TransitionMoment(k, m_c, O)
         / ((w_n - w_o) * (w_k - w_2))
     ),
     "ac": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, k)
-        * TransitionMoment(k, opm_c, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, k)
+        * TransitionMoment(k, m_c, O)
         / ((w_n - w_o) * (w_k - w_2))
     ),
     "bc": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, opm_b, k)
-        * TransitionMoment(k, opm_c, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, m_b, k)
+        * TransitionMoment(k, m_c, O)
         / ((w_n - w_o) * (w_k - w_2))
     ),
     "abc": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, opm_b, k)
-        * TransitionMoment(k, opm_c, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, m_b, k)
+        * TransitionMoment(k, m_c, O)
         / ((w_n - w_o) * (w_k - w_2))
     ),
 }
@@ -105,66 +105,66 @@ def run_scf(molecule, basis, backend="pyscf"):
 
 SOS_gamma_like = {
     "a": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, op_d, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, mu_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "b": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, opm_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, op_d, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, m_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, mu_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "d": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, opm_d, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, m_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "ac": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, opm_c, p)
-        * TransitionMoment(p, op_d, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, m_c, p)
+        * TransitionMoment(p, mu_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "ad": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, opm_d, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, m_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "cd": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, opm_c, p)
-        * TransitionMoment(p, opm_d, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, m_c, p)
+        * TransitionMoment(p, m_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "abd": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, opm_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, opm_d, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, m_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, m_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "bcd": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, opm_b, m)
-        * TransitionMoment(m, opm_c, p)
-        * TransitionMoment(p, opm_d, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, m_b, m)
+        * TransitionMoment(m, m_c, p)
+        * TransitionMoment(p, m_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
     "abcd": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, opm_b, m)
-        * TransitionMoment(m, opm_c, p)
-        * TransitionMoment(p, opm_d, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, m_b, m)
+        * TransitionMoment(m, m_c, p)
+        * TransitionMoment(p, m_d, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
     ),
 }
@@ -172,43 +172,43 @@ def run_scf(molecule, basis, backend="pyscf"):
 
 SOS_delta_like = {
     "a": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, op_d, k)
-        * TransitionMoment(k, op_e, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, mu_d, k)
+        * TransitionMoment(k, mu_e, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3) * (w_k - w_2))
     ),
     "b": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, opm_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, op_d, k)
-        * TransitionMoment(k, op_e, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, m_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, mu_d, k)
+        * TransitionMoment(k, mu_e, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3) * (w_k - w_2))
     ),
     "e": (
-        TransitionMoment(O, op_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, op_d, k)
-        * TransitionMoment(k, opm_e, O)
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, mu_d, k)
+        * TransitionMoment(k, m_e, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3) * (w_k - w_2))
     ),
     "ae": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, op_b, m)
-        * TransitionMoment(m, op_c, p)
-        * TransitionMoment(p, op_d, k)
-        * TransitionMoment(k, opm_e, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, mu_b, m)
+        * TransitionMoment(m, mu_c, p)
+        * TransitionMoment(p, mu_d, k)
+        * TransitionMoment(k, m_e, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3) * (w_k - w_2))
     ),
     "abcde": (
-        TransitionMoment(O, opm_a, n)
-        * TransitionMoment(n, opm_b, m)
-        * TransitionMoment(m, opm_c, p)
-        * TransitionMoment(p, opm_d, k)
-        * TransitionMoment(k, opm_e, O)
+        TransitionMoment(O, m_a, n)
+        * TransitionMoment(n, m_b, m)
+        * TransitionMoment(m, m_c, p)
+        * TransitionMoment(p, m_d, k)
+        * TransitionMoment(k, m_e, O)
         / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3) * (w_k - w_2))
     ),
 }
@@ -314,14 +314,14 @@ def test_h2o_sto3g_adc2(self):
         scfres = run_scf(molecule, basis)
         refstate = adcc.ReferenceState(scfres)
         term = (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, m)
-            * TransitionMoment(m, opm_c, p)
-            * TransitionMoment(p, opm_d, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, m)
+            * TransitionMoment(m, m_c, p)
+            * TransitionMoment(p, m_d, O)
             / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
         )
         sos = SumOverStates(
-            term, [n, m, p], perm_pairs=[(op_a, -w_o), (op_b, w_1), (opm_c, w_2), (opm_d, w_3)]
+            term, [n, m, p], perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (m_c, w_2), (m_d, w_3)]
         )
         mock_state = cache.data_fulldiag[case]
         state = adcc.run_adc(refstate, method=method, n_singlets=5)
diff --git a/responsefun/test_property.py b/responsefun/test_property.py
index bea7a2c..2ad7b2d 100644
--- a/responsefun/test_property.py
+++ b/responsefun/test_property.py
@@ -24,14 +24,13 @@
     gamma,
     k,
     m,
+    mu_a,
+    mu_b,
+    mu_c,
+    mu_d,
     n,
-    op_a,
-    op_b,
-    op_c,
-    op_d,
     p,
     w,
-    w_prime,
     w_1,
     w_2,
     w_3,
@@ -41,6 +40,7 @@
     w_n,
     w_o,
     w_p,
+    w_prime,
 )
 from responsefun.testdata import cache
 from responsefun.testdata.static_data import xyz
@@ -59,85 +59,85 @@ def run_scf(molecule, basis, backend="pyscf"):
 SOS_expressions = {
     "alpha": (
         (
-            TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w)
-            + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, O) / (w_n + w)
+            TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w)
+            + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n + w)
         ),
         None,
     ),
     "alpha_complex": (
         (
-            TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w - 1j * gamma)
-            + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, O) / (w_n + w + 1j * gamma)
+            TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w - 1j * gamma)
+            + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n + w + 1j * gamma)
         ),
         None,
     ),
     "rixs_short": (
-        (TransitionMoment(f, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w - 1j * gamma)),
+        (TransitionMoment(f, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w - 1j * gamma)),
         None,
     ),
     "rixs": (
         (
-            (TransitionMoment(f, op_a, n) * TransitionMoment(n, op_b, O)
+            (TransitionMoment(f, mu_a, n) * TransitionMoment(n, mu_b, O)
             / (w_n - w - 1j * gamma))
-            + (TransitionMoment(f, op_b, n) * TransitionMoment(n, op_a, O)
+            + (TransitionMoment(f, mu_b, n) * TransitionMoment(n, mu_a, O)
             / (w_n + w - w_f + 1j * gamma))
         ),
         None,
     ),
     "tpa_resonant": (
         (
-            TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, f) / (w_n - (w_f / 2))
-            + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, f) / (w_n - (w_f / 2))
+            TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, f) / (w_n - (w_f / 2))
+            + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, f) / (w_n - (w_f / 2))
         ),
         None,
     ),
     "beta": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, k)
-            * TransitionMoment(k, op_c, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, k)
+            * TransitionMoment(k, mu_c, O)
             / ((w_n - w_o) * (w_k - w_2))
         ),
-        [(op_a, -w_o), (op_b, w_1), (op_c, w_2)],
+        [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],
     ),
     "beta_complex": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, k, shifted=True)
-            * TransitionMoment(k, op_c, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, k, shifted=True)
+            * TransitionMoment(k, mu_c, O)
             / ((w_n - w_o - 1j * gamma) * (w_k - w_2 - 1j * gamma))
         ),
-        [(op_a, -w_o - 1j * gamma), (op_b, w_1 + 1j * gamma), (op_c, w_2 + 1j * gamma)],
+        [(mu_a, -w_o - 1j * gamma), (mu_b, w_1 + 1j * gamma), (mu_c, w_2 + 1j * gamma)],
     ),
     "gamma": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, m)
-            * TransitionMoment(m, op_c, p)
-            * TransitionMoment(p, op_d, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, m)
+            * TransitionMoment(m, mu_c, p)
+            * TransitionMoment(p, mu_d, O)
             / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
         ),
-        [(op_a, -w_o), (op_b, w_1), (op_c, w_2), (op_d, w_3)],
+        [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2), (mu_d, w_3)],
     ),
     "gamma_extra_terms_1": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, O)
-            * TransitionMoment(O, op_c, m)
-            * TransitionMoment(m, op_d, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, O)
+            * TransitionMoment(O, mu_c, m)
+            * TransitionMoment(m, mu_d, O)
             / ((w_n - w_o) * (w_m - w_3) * (w_m + w_2))
         ),
-        [(op_a, -w_o), (op_b, w_1), (op_c, w_2), (op_d, w_3)],
+        [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2), (mu_d, w_3)],
     ),
     "gamma_extra_terms_2": (
         (
-            TransitionMoment(O, op_a, n)
-            * TransitionMoment(n, op_b, O)
-            * TransitionMoment(O, op_c, m)
-            * TransitionMoment(m, op_d, O)
+            TransitionMoment(O, mu_a, n)
+            * TransitionMoment(n, mu_b, O)
+            * TransitionMoment(O, mu_c, m)
+            * TransitionMoment(m, mu_d, O)
             / ((w_n - w_o) * (-w_2 - w_3) * (w_m - w_3))
         ),
-        [(op_a, -w_o), (op_b, w_1), (op_c, w_2), (op_d, w_3)],
+        [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2), (mu_d, w_3)],
     ),
 }
 
@@ -320,7 +320,7 @@ def test_rixs_short(self, case):
                     err_msg = "w = {}, gamma = {}, final_state = {}".format(
                         tup[0][1], tup[1], final_state
                     )
-                    assert_allclose_signfix(rixs_isr, rixs_sos, atol=1e-7, err_msg=err_msg)
+                    assert_allclose_signfix(rixs_isr, rixs_sos, atol=5e-7, err_msg=err_msg)
 
     def test_rixs(self, case):
         molecule, basis, method = case.split("_")
@@ -457,7 +457,7 @@ def test_rixs_short(self, case):
                     err_msg = "w = {}, gamma = {}, final_state = {}".format(
                         tup[0][1], tup[1], final_state
                     )
-                    assert_allclose_signfix(rixs_isr, rixs_sos, atol=1e-7, err_msg=err_msg)
+                    assert_allclose_signfix(rixs_isr, rixs_sos, atol=5e-7, err_msg=err_msg)
 
     def test_rixs(self, case):
         molecule, basis, method = case.split("_")
diff --git a/responsefun/test_twodimensional_operators.py b/responsefun/test_twodimensional_operators.py
new file mode 100644
index 0000000..9d2e701
--- /dev/null
+++ b/responsefun/test_twodimensional_operators.py
@@ -0,0 +1,115 @@
+import adcc
+import numpy as np
+import pytest
+
+from responsefun.evaluate_property import (
+    evaluate_property_isr,
+    evaluate_property_sos_fast,
+)
+from responsefun.SumOverStates import TransitionMoment
+from responsefun.symbols_and_labels import (
+    O,
+    k,
+    mu_a,
+    mu_b,
+    mu_c,
+    n,
+    q_ab,
+    q_bc,
+    q_cd,
+    q_de,
+    w,
+    w_1,
+    w_2,
+    w_k,
+    w_n,
+    w_o,
+)
+from responsefun.testdata import cache
+from responsefun.testdata.static_data import xyz
+
+
+def run_scf(molecule, basis, backend="pyscf"):
+    scfres = adcc.backends.run_hf(
+        backend,
+        xyz=xyz[molecule],
+        basis=basis,
+    )
+    return scfres
+
+
+SOS_alpha_like = {
+    "ab": (
+        TransitionMoment(O, q_ab, n) * TransitionMoment(n, mu_c, O) / (w_n - w)
+        + TransitionMoment(O, mu_c, n) * TransitionMoment(n, q_ab, O) / (w_n + w)
+    ),
+    "bc": (
+        TransitionMoment(O, mu_a, n) * TransitionMoment(n, q_bc, O) / (w_n - w)
+        + TransitionMoment(O, q_bc, n) * TransitionMoment(n, mu_a, O) / (w_n + w)
+    ),
+    "abcd": (
+        TransitionMoment(O, q_ab, n) * TransitionMoment(n, q_cd, O) / (w_n - w)
+        + TransitionMoment(O, q_cd, n) * TransitionMoment(n, q_ab, O) / (w_n + w)
+    ),
+}
+
+
+SOS_beta_like = {
+    "ab": (
+        TransitionMoment(O, q_ab, n)
+        * TransitionMoment(n, mu_b, k)
+        * TransitionMoment(k, mu_c, O)
+        / ((w_n - w_o) * (w_k - w_2))
+    ),
+    "cd": (
+        TransitionMoment(O, mu_a, n)
+        * TransitionMoment(n, mu_b, k)
+        * TransitionMoment(k, q_cd, O)
+        / ((w_n - w_o) * (w_k - w_2))
+    ),
+    "abde": (
+        TransitionMoment(O, q_ab, n)
+        * TransitionMoment(n, mu_c, k)
+        * TransitionMoment(k, q_de, O)
+        / ((w_n - w_o) * (w_k - w_2))
+    ),
+}
+
+@pytest.mark.parametrize("ops", SOS_alpha_like.keys())
+class TestAlphaLike:
+    def test_h2o_sto3g_adc2(self, ops):
+        case = "h2o_sto3g_adc2"
+        if case not in cache.data_fulldiag:
+            pytest.skip(f"{case} cache file not available.")
+        molecule, basis, method = case.split("_")
+        scfres = run_scf(molecule, basis)
+        refstate = adcc.ReferenceState(scfres)
+        expr = SOS_alpha_like[ops]
+        mock_state = cache.data_fulldiag[case]
+        state = adcc.run_adc(refstate, method=method, n_singlets=5)
+        freq = (w, 0.5)
+        alpha_sos = evaluate_property_sos_fast(mock_state, expr, [n], freqs_in=freq, freqs_out=freq)
+        alpha_isr = evaluate_property_isr(state, expr, [n], freqs_in=freq, freqs_out=freq)
+        np.testing.assert_allclose(alpha_isr, alpha_sos, atol=1e-8)
+
+
+@pytest.mark.parametrize("ops", SOS_beta_like.keys())
+class TestBetaLike:
+    def test_h2o_sto3g_adc2(self, ops):
+        case = "h2o_sto3g_adc2"
+        if case not in cache.data_fulldiag:
+            pytest.skip(f"{case} cache file not available.")
+        molecule, basis, method = case.split("_")
+        scfres = run_scf(molecule, basis)
+        refstate = adcc.ReferenceState(scfres)
+        expr = SOS_beta_like[ops]
+        mock_state = cache.data_fulldiag[case]
+        state = adcc.run_adc(refstate, method=method, n_singlets=5)
+
+        freqs_in = [(w_1, 0.5), (w_2, 0.5)]
+        freqs_out = (w_o, 1)
+        beta_sos = evaluate_property_sos_fast(mock_state, expr, [n, k], freqs_in=freqs_in,
+                                              freqs_out=freqs_out, extra_terms=False)
+        beta_isr = evaluate_property_isr(state, expr, [n, k], freqs_in=freqs_in,
+                                         freqs_out=freqs_out, extra_terms=False)
+        np.testing.assert_allclose(beta_isr, beta_sos, atol=1e-8)
\ No newline at end of file
diff --git a/responsefun/testdata/0_download_testdata.sh b/responsefun/testdata/0_download_testdata.sh
index e6cd3d2..bfd6b90 100755
--- a/responsefun/testdata/0_download_testdata.sh
+++ b/responsefun/testdata/0_download_testdata.sh
@@ -2,9 +2,9 @@
 
 # taken from respondo
 
-SOURCE="https://wwwagdreuw.iwr.uni-heidelberg.de/responsefun_test_data/0.1.0"
+SOURCE="https://wwwagdreuw.iwr.uni-heidelberg.de/responsefun_test_data/0.2.0"
 DATAFILES=(
-	data_0.1.0.tar.gz
+	data_0.2.0.tar.gz
 )
 #
 # -----
diff --git a/responsefun/testdata/SHA256SUMS b/responsefun/testdata/SHA256SUMS
index ca2f257..ab35df9 100644
--- a/responsefun/testdata/SHA256SUMS
+++ b/responsefun/testdata/SHA256SUMS
@@ -1 +1 @@
-7530b9a13c5ad9fccd25f23981e291f5baf158dfa4b939c6c20a3b2db1d2e184  data_0.1.0.tar.gz
\ No newline at end of file
+b1604e8f7aa18e0dd01ce785166bbaf51c71ba48386c565e45707f0d06babc4b  data_0.2.0.tar.gz
\ No newline at end of file
diff --git a/responsefun/testdata/dump_full_diagonalization.py b/responsefun/testdata/dump_full_diagonalization.py
index a961db2..f3f67ff 100644
--- a/responsefun/testdata/dump_full_diagonalization.py
+++ b/responsefun/testdata/dump_full_diagonalization.py
@@ -11,6 +11,7 @@
 
 
 def main():
+    gauge_origin = "origin"
     for case in cases:
         n_singlets = cases[case]
         molecule, basis, method = case.split("_")
@@ -22,9 +23,10 @@ def main():
             # multiplicity=multiplicity,
             # conv_tol_grad=conv_tol_grad,
         )
-        state = adcc.run_adc(method=method, data_or_matrix=scfres, n_singlets=n_singlets)
+        state = adcc.run_adc(method=method, data_or_matrix=scfres, n_singlets=n_singlets,
+                             n_guesses=n_singlets)
         dips = state.reference_state.operators.electric_dipole
-        mdips = state.reference_state.operators.magnetic_dipole
+        mdips = state.reference_state.operators.magnetic_dipole(gauge_origin)
 
         # state to state transition moments
         s2s_tdms = np.zeros((state.size, state.size, 3))
@@ -55,6 +57,12 @@ def main():
                 d = getattr(state, key)
             except NotImplementedError:
                 continue
+            if callable(d):
+                try:
+                    d = d(gauge_origin)
+                except NotImplementedError:
+                    # some properties are not available for every backend
+                    continue
             if not isinstance(d, np.ndarray):
                 continue
             if not np.issubdtype(d.dtype, np.number):