diff --git a/ci_env.yml b/ci_env.yml index d2a9d37..0bd5cfb 100644 --- a/ci_env.yml +++ b/ci_env.yml @@ -9,6 +9,6 @@ dependencies: - pytest - pytest-cov - coveralls - - zarr + - zarr>=2,<3 - ruff - isort \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index eee3015..3a8652f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,7 +42,7 @@ test = [ "pytest-cov", "pytest-runner", "pyscf", - "zarr", + "zarr>=2,<3", ] [tool.setuptools] diff --git a/respondo/mcd.py b/respondo/mcd.py index 308b65d..ada46e1 100644 --- a/respondo/mcd.py +++ b/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,8 +44,8 @@ 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)) @@ -102,5 +53,7 @@ def projection(X, bl=None): 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 diff --git a/respondo/solve_response.py b/respondo/solve_response.py index 6902018..9726b32 100644 --- a/respondo/solve_response.py +++ 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): diff --git a/respondo/sos.py b/respondo/sos.py index cd08c4b..9dc551e 100644 --- a/respondo/sos.py +++ b/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 diff --git a/respondo/testdata/0_download_testdata.sh b/respondo/testdata/0_download_testdata.sh index 7abf607..0cff9ab 100755 --- a/respondo/testdata/0_download_testdata.sh +++ 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 ) # # ----- diff --git a/respondo/testdata/SHA256SUMS b/respondo/testdata/SHA256SUMS index ef3ec4a..bdca93a 100644 --- a/respondo/testdata/SHA256SUMS +++ b/respondo/testdata/SHA256SUMS @@ -1 +1 @@ -63522c1f857212c5622de5b18f89e9c8c1d9676f5b5621b24f9508fc48d5d6ef data_0.0.1.tar.gz +ed7bf0506c2c9d299529a56a5b4a6c76b9985b2bd2954b8be13abc33375342ed data_0.2.0.tar.gz \ No newline at end of file diff --git a/respondo/testdata/dump_full_diagonalization.py b/respondo/testdata/dump_full_diagonalization.py index d9cb6ab..ad52b4d 100644 --- a/respondo/testdata/dump_full_diagonalization.py +++ b/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)