Fix MCD (#15)

* fix mcd_bterm

* pin zarr 2.0

* fix sos_mcd_bterm

* update testdata with fixed magnetic tdms

* add todo for transition_polarizability

* fix ruff

* fix ci

Author: apapapostolou

ci_env.yml

@@ -9,6 +9,6 @@ dependencies:
   - pytest
   - pytest-cov
   - coveralls
-  - zarr
+  - zarr>=2,<3
   - ruff
   - isort

pyproject.toml

@@ -42,7 +42,7 @@ test = [
   "pytest-cov",
   "pytest-runner",
   "pyscf",
-  "zarr",
+  "zarr>=2,<3",
 ]
 
 [tool.setuptools]

respondo/mcd.py

@@ -1,43 +1,12 @@
 import numpy as np
-from adcc import AmplitudeVector
 from adcc.adc_pp.modified_transition_moments import modified_transition_moments
 from adcc.Excitation import Excitation
-from adcc.functions import einsum
 from adcc.workflow import construct_adcmatrix
 
 from .misc import select_property_method
 from .solve_response import solve_response, transition_polarizability
 
 
-def compute_adc1_f1_mag(magdip, ground_state):
-    mtm = magdip.ov + einsum("ijab,jb->ia", ground_state.t2("o1o1v1v1"), magdip.ov)
-    return AmplitudeVector(ph=mtm)
-
-
-def compute_adc2_f1_mag(magdip, ground_state):
-    t2 = ground_state.t2("o1o1v1v1")
-    td2 = ground_state.td2("o1o1v1v1")
-    p0 = ground_state.mp2_diffdm
-    d = magdip
-    return (
-        d.ov
-        + 1.0 * einsum("ijab,jb->ia", t2, d.ov + 0.5 * einsum("jkbc,kc->jb", t2, d.ov))
-        + 0.5 * (einsum("ij,ja->ia", p0.oo, d.ov) - 1.0 * einsum("ib,ab->ia", d.ov, p0.vv))
-        - 1.0 * einsum("ib,ab->ia", p0.ov, d.vv)
-        - 1.0 * einsum("ij,ja->ia", d.oo, p0.ov)
-        + 1.0 * einsum("ijab,jb->ia", td2, d.ov)
-    )
-
-
-def compute_adc2_f2_mag(magdip, ground_state):
-    t2 = ground_state.t2("o1o1v1v1")
-    term1 = -1.0 * einsum("ijac,bc->ijab", t2, magdip.vv)
-    term2 = -1.0 * einsum("ik,kjab->ijab", magdip.oo, t2)
-    term1 = term1.antisymmetrise(2, 3)
-    term2 = term2.antisymmetrise(0, 1)
-    return term1 - term2
-
-
 def mcd_bterm(state, property_method=None, **solver_args):
     if not isinstance(state, Excitation):
         raise TypeError()
@@ -49,30 +18,12 @@ def mcd_bterm(state, property_method=None, **solver_args):
 
     dips_el = hf.operators.electric_dipole
     dips_mag = hf.operators.magnetic_dipole
-    # Electric dipole right-hand-side
     rhss_el = modified_transition_moments(property_method, mp, dips_el)
+    rhss_mag = modified_transition_moments(property_method, mp, dips_mag)
 
-    # Magnetic dipole right-hand-side
-    # TODO: temporary hack, add to adcc...
-    if property_method.name == "adc0":
-        rhss_mag = modified_transition_moments(property_method, mp, dips_mag)
-        rhss_mag = [-1.0 * rhs_mag for rhs_mag in rhss_mag]
-    elif property_method.name == "adc1":
-        rhss_mag = [-1.0 * compute_adc1_f1_mag(mag, mp) for mag in dips_mag]
-    elif property_method.name == "adc2":
-        rhss_mag = [
-            -1.0
-            * AmplitudeVector(
-                ph=compute_adc2_f1_mag(mag, mp),
-                pphh=compute_adc2_f2_mag(mag, mp),
-            )
-            for mag in dips_mag
-        ]
-    else:
-        raise NotImplementedError("")
-
+    # the minus sign is required due to the anti-hermiticity of the magnetic dipole operator
     response_mag = [
-        solve_response(matrix, rhs_mag, omega=0.0, gamma=0.0, **solver_args) for rhs_mag in rhss_mag
+        -1.0 * solve_response(matrix, rhs_mag, omega=0.0, gamma=0.0, **solver_args) for rhs_mag in rhss_mag
     ]
 
     v_f = state.excitation_vector
@@ -93,14 +44,16 @@ def projection(X, bl=None):
         for rhs_el in rhss_el
     ]
 
-    term1 = -transition_polarizability(property_method, mp, response_mag, dips_el, v_f)
-    term2 = -transition_polarizability(property_method, mp, v_f, dips_mag, response_el)
+    term1 = transition_polarizability(property_method, mp, v_f, dips_el, response_mag)
+    term2 = transition_polarizability(property_method, mp, v_f, dips_mag, response_el)
 
     # Levi-Civita tensor
     epsilon = np.zeros((3, 3, 3))
     epsilon[0, 1, 2] = epsilon[1, 2, 0] = epsilon[2, 0, 1] = 1
     epsilon[2, 1, 0] = epsilon[0, 2, 1] = epsilon[1, 0, 2] = -1
 
     tdip_f = state.transition_dipole_moment
-    B = np.einsum("abc,a,bc->", epsilon, tdip_f, term1 + term2)
+    # the minus sign accounts for the negative charge, since it is not included in the operators
+    # TODO as soon as PR #190 in adcc is merged: remove minus
+    B = -1.0 * np.einsum("abc,a,bc->", epsilon, tdip_f, term1 + np.transpose(term2))
     return B

respondo/solve_response.py

@@ -89,6 +89,8 @@ def solve_response(
         )
 
 
+# TODO: fix transition_polarizability and transition_polarizability_complex
+# by swapping from_vecs and to_vecs (for ADC it is defined the other way around)
 # from_vecs * B(ops) * to_vecs
 def transition_polarizability(method, ground_state, from_vecs, ops, to_vecs):
     if not isinstance(from_vecs, list):

respondo/sos.py

@@ -112,19 +112,19 @@ def sos_mcd_bterm(state, final_state=0):
         for A in range(3):
             for B in range(3):
                 term1[A, B] -= (
-                    state.transition_magnetic_dipole_moment[k][A]
-                    * s2s_tdm[k, final_state, B]
+                    state.transition_magnetic_dipole_moment[k][B]
+                    * s2s_tdm[final_state, k, A]
                     / (state.excitation_energy_uncorrected[k])
                 )
                 if k != final_state:
                     term2[A, B] += (
-                        state.transition_dipole_moment[k][B]
-                        * s2s_tdm_mag[k, final_state, A]
+                        state.transition_dipole_moment[k][A]
+                        * s2s_tdm_mag[final_state, k, B]
                         / (state.excitation_energy_uncorrected[k] - e_f)
                     )
 
     epsilon = np.zeros((3, 3, 3))
     epsilon[0, 1, 2] = epsilon[1, 2, 0] = epsilon[2, 0, 1] = 1
     epsilon[2, 1, 0] = epsilon[0, 2, 1] = epsilon[1, 0, 2] = -1
-    B = np.einsum("abc,a,bc->", epsilon, tdip_f, term1 + term2)
+    B = -1.0 * np.einsum("abc,a,bc->", epsilon, tdip_f, term1 + term2)
     return B

respondo/testdata/0_download_testdata.sh

@@ -2,9 +2,9 @@
 
 # taken from adcc repository, written by @mfherbst
 
-SOURCE="https://wwwagdreuw.iwr.uni-heidelberg.de/respondo_test_data/0.1.0"
+SOURCE="https://wwwagdreuw.iwr.uni-heidelberg.de/respondo_test_data/0.2.0"
 DATAFILES=(
-	data_0.0.1.tar.gz
+	data_0.2.0.tar.gz
 )
 #
 # -----

respondo/testdata/SHA256SUMS

@@ -1 +1 @@
-63522c1f857212c5622de5b18f89e9c8c1d9676f5b5621b24f9508fc48d5d6ef  data_0.0.1.tar.gz
+ed7bf0506c2c9d299529a56a5b4a6c76b9985b2bd2954b8be13abc33375342ed  data_0.2.0.tar.gz

respondo/testdata/dump_full_diagonalization.py

@@ -47,7 +47,7 @@ def main():
         z.create_group("excitation")
         exci = z["excitation"]
         propkeys = state.excitation_property_keys
-        propkeys.extend(state.excitation_energy_corrections.keys())
+        propkeys.extend([k.name for k in state._excitation_energy_corrections])
         for key in propkeys:
             try:
                 d = getattr(state, key)