From f9650fccf2e1a63d87ca1d3ecaa41e695fe50b0b Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Fri, 14 Mar 2025 13:10:56 +0100 Subject: [PATCH 01/13] new operators --- responsefun/AdccProperties.py | 172 ++++++++++++++++++++++++++---- responsefun/symbols_and_labels.py | 30 +++++- 2 files changed, 182 insertions(+), 20 deletions(-) diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 5d1de6e..69b413a 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -64,6 +64,12 @@ class Operator: dim=1, is_imag=False, ), + Operator( + name="electric_dipole_velocity", + symbol="mu_p", + symmetry=Symmetry.ANTIHERMITIAN, + dim=1 + ), Operator( name="magnetic_dipole", symbol="m", @@ -71,6 +77,18 @@ class Operator: dim=1, is_imag=True, ), + Operator( + name="electric_quadrupole", + symbol="q", + symmetry=Symmetry.HERMITIAN, + dim=2 + ), + Operator( + name="electric_quadrupole_velocity", + symbol="q_p", + symmetry=Symmetry.ANTIHERMITIAN, + dim=2 + ), Operator( name="diamagnetic_magnetizability", symbol="xi", @@ -142,12 +160,40 @@ 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. + """ + # Note that the charge adcc/PR#190 must already be contained in the integrals + # passed to this function. + + ref_mom = np.array( + [product_trace(ints, state.ground_state.reference_state.density) + for ints in integrals] + ) + if pm_level == 1: + gs_mom = ref_mom + elif pm_level == 2: + mp2corr = -1.0 * np.array([product_trace(ints, state.ground_state.mp2_diffdm) + for ints in integrals]) + 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]] = "mass_center"): self._state = state self._state_size = len(state.excitation_energy_uncorrected) self._property_method = self._state.property_method @@ -340,27 +386,14 @@ def _operator(self) -> Operator: @property def integrals(self) -> list[adcc.OneParticleOperator]: - return self._state.reference_state.operators.magnetic_dipole + 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 - else: - raise NotImplementedError( - "Only magnetic dipole moments for level 1 and 2 are implemented." - ) + # Note: This next line is needed because of missing minus sign in the definition + # of the operators, TODO: remove if adcc/PR#190 is merged + ints = [-1.0 * ints for ints in self.integrals] + return compute_ground_state_moment(self._state, ints, self._pm_level) def _transition_moment(self) -> np.ndarray: return self._state.transition_magnetic_dipole_moment @@ -368,5 +401,106 @@ def _transition_moment(self) -> np.ndarray: def _state_to_state_transition_moment(self) -> np.ndarray: if isinstance(self._state, MockExcitedStates): return self._state.transition_magnetic_moment_s2s + else: + return compute_state_to_state_transition_moments(self._state, self.integrals) + + +class ElectricDipoleVelocity(AdccProperties): + @property + def _operator(self) -> Operator: + return get_operator_by_name("electric_dipole_velocity") + + @property + def integrals(self) -> list[adcc.OneParticleOperator]: + return self._state.reference_state.operators.electric_dipole_velocity(self._gauge_origin) + + @property + def gs_moment(self) -> np.ndarray: + raise NotImplementedError("Ground-state moments not implemented for electric dipole " + "operator in velocity gauge.") + + def _transition_moment(self) -> np.ndarray: + return self._state.transition_dipole_moment_velocity + + def _state_to_state_transition_moment(self) -> np.ndarray: + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError() + else: + return compute_state_to_state_transition_moments(self._state, self.integrals) + + +class ElectricQuadrupole(AdccProperties): + @property + def _operator(self) -> Operator: + return get_operator_by_name("electric_quadrupole") + + @property + def integrals(self) -> list[adcc.OneParticleOperator]: + return self._state.reference_state.operators.magnetic_dipole + + @property + def gs_moment(self) -> np.ndarray: + nuclear_contribution = self._state.reference_state.nuclear_quadrupole(self._gauge_origin) + # Note: This next line is needed because of missing minus sign in the definition + # of the operators, TODO: remove if adcc/PR#190 is merged + ints = [-1.0 * ints for ints in self.integrals] + return compute_ground_state_moment(self._state, ints, self._pm_level, + nuclear_contribution=nuclear_contribution) + + + def _transition_moment(self) -> np.ndarray: + return self._state.transition_dipole_moment_velocity + + def _state_to_state_transition_moment(self) -> np.ndarray: + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError() + else: + return compute_state_to_state_transition_moments(self._state, self.integrals) + + +class ElectricQuadrupoleVelocity(AdccProperties): + @property + def _operator(self) -> Operator: + return get_operator_by_name("electric_quadrupole_velocity") + + @property + def integrals(self) -> list[adcc.OneParticleOperator]: + 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(self) -> np.ndarray: + return self._state.transition_quadrupole_moment_velocity(self._gauge_origin) + + def _state_to_state_transition_moment(self) -> np.ndarray: + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError() + else: + return compute_state_to_state_transition_moments(self._state, self.integrals) + + +class DiamagneticMagnetizability(AdccProperties): + @property + def _operator(self) -> Operator: + return get_operator_by_name("diamagentic_magnetizability") + + @property + def integrals(self) -> list[adcc.OneParticleOperator]: + return self._state.reference_state.operators.diamagentic_magnetizability(self._gauge_origin) + + @property + def gs_moment(self) -> np.ndarray: + return compute_ground_state_moment(self._state, self.integrals, self._pm_level) + + def _transition_moment(self) -> np.ndarray: + return compute_transition_moments(self._state, self.integrals) + + def _state_to_state_transition_moment(self) -> np.ndarray: + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError() else: return compute_state_to_state_transition_moments(self._state, self.integrals) \ No newline at end of file diff --git a/responsefun/symbols_and_labels.py b/responsefun/symbols_and_labels.py index 2e61d17..a583337 100644 --- a/responsefun/symbols_and_labels.py +++ b/responsefun/symbols_and_labels.py @@ -30,9 +30,37 @@ op_d = OneParticleOperator("D", "electric_dipole", False) op_e = OneParticleOperator("E", "electric_dipole", False) +# electric dipole operators in velocity gauge +op_vel_a = OneParticleOperator("A", "electric_dipole_velocity", False) +op_vel_b = OneParticleOperator("B", "electric_dipole_velocity", False) +op_vel_c = OneParticleOperator("C", "electric_dipole_velocity", False) +op_vel_d = OneParticleOperator("D", "electric_dipole_velocity", False) +op_vel_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 +opm_e = OneParticleOperator("E", "magnetic_dipole", False) + +# electric quadrupole operators +opq_ab = OneParticleOperator("AB", "electric_quadrupole", False) +opq_bc = OneParticleOperator("BC", "electric_quadrupole", False) +opq_cd = OneParticleOperator("CD", "electric_quadrupole", False) +opq_de = OneParticleOperator("DE", "electric_quadrupole", False) +opq_ef = OneParticleOperator("EF", "electric_quadrupole", False) + +# electric quadrupole operators in velocity gauge +opq_vel_ab = OneParticleOperator("AB", "electric_quadrupole_velocity", False) +opq_vel_bc = OneParticleOperator("BC", "electric_quadrupole_velocity", False) +opq_vel_cd = OneParticleOperator("CD", "electric_quadrupole_velocity", False) +opq_vel_de = OneParticleOperator("DE", "electric_quadrupole_velocity", False) +opq_vel_ef = OneParticleOperator("EF", "electric_quadrupole_velocity", False) + +# diamagnetic magnetizability operator +opxi_ab = OneParticleOperator("AB", "diamagnetic_magnetizability", False) +opxi_bc = OneParticleOperator("BC", "diamagnetic_magnetizability", False) +opxi_cd = OneParticleOperator("CD", "diamagnetic_magnetizability", False) +opxi_de = OneParticleOperator("DE", "diamagnetic_magnetizability", False) +opxi_ef = OneParticleOperator("EF", "diamagnetic_magnetizability", False) \ No newline at end of file From 4a1f9a8ec01c1e5c2a594355e2d4aa11616d6eff Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Tue, 18 Mar 2025 13:40:35 +0100 Subject: [PATCH 02/13] gauge origin selection new operators, some minor testing --- responsefun/AdccProperties.py | 123 +++++++++++++----- responsefun/SumOverStates.py | 4 +- responsefun/evaluate_property.py | 12 +- responsefun/test_twodimensional_operators.py | 116 +++++++++++++++++ .../testdata/dump_full_diagonalization.py | 9 +- 5 files changed, 222 insertions(+), 42 deletions(-) create mode 100644 responsefun/test_twodimensional_operators.py diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 69b413a..5be708c 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -161,22 +161,28 @@ def compute_state_to_state_transition_moments(state, integrals, initial_state=No def compute_ground_state_moment(state, integrals, pm_level, - nuclear_contribution = None) -> np.ndarray: + nuclear_contribution=None) -> np.ndarray: """ Computes the ground state moments. """ - # Note that the charge adcc/PR#190 must already be contained in the integrals + # Note that the charge (see adcc/PR#190) must already be contained in the integrals # passed to this function. - ref_mom = np.array( - [product_trace(ints, state.ground_state.reference_state.density) - for ints in integrals] - ) + # 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 == 1: gs_mom = ref_mom elif pm_level == 2: - mp2corr = -1.0 * np.array([product_trace(ints, state.ground_state.mp2_diffdm) - for ints in integrals]) + 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( @@ -193,7 +199,7 @@ class AdccProperties(ABC): from adcc for a given operator.""" def __init__(self, state: Union[adcc.ExcitedStates, MockExcitedStates], - gauge_origin: Union[str, tuple[float, float, float]] = "mass_center"): + 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 @@ -288,14 +294,18 @@ def _state_to_state_transition_moment(self) -> np.ndarray: def modified_transition_moments( self, comp: Union[int, None] = None - ) -> Union[adcc.AmplitudeVector, list[adcc.AmplitudeVector]]: + ) -> Union[adcc.AmplitudeVector, np.ndarray[adcc.AmplitudeVector]]: if comp is None: op = 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] + if isinstance(op, adcc.OneParticleOperator): + mtms = modified_transition_moments( + self._property_method, self._state.ground_state, op + ) + elif isinstance(op, list): + mtms = np.array([(modified_transition_moments(self._property_method, + self._state.ground_state, op_ints)) for op_ints in self.integrals]) return mtms def modified_transition_moments_reverse( @@ -347,8 +357,16 @@ def build_adcc_properties( ) -> 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 @@ -360,7 +378,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: @@ -386,17 +407,23 @@ def _operator(self) -> Operator: @property def integrals(self) -> list[adcc.OneParticleOperator]: - return self._state.reference_state.operators.magnetic_dipole(self._gauge_origin) + 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: - # Note: This next line is needed because of missing minus sign in the definition - # of the operators, TODO: remove if adcc/PR#190 is merged - ints = [-1.0 * ints for ints in self.integrals] - return compute_ground_state_moment(self._state, ints, self._pm_level) + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError + else: + return compute_ground_state_moment(self._state, self.integrals, self._pm_level) def _transition_moment(self) -> np.ndarray: - return self._state.transition_magnetic_dipole_moment + if isinstance(self._state, MockExcitedStates): + return self._state.transition_magnetic_dipole_moment + else: + return self._state.transition_magnetic_dipole_moment(self._gauge_origin) def _state_to_state_transition_moment(self) -> np.ndarray: if isinstance(self._state, MockExcitedStates): @@ -412,7 +439,10 @@ def _operator(self) -> Operator: @property def integrals(self) -> list[adcc.OneParticleOperator]: - return self._state.reference_state.operators.electric_dipole_velocity(self._gauge_origin) + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError + else: + return self._state.reference_state.operators.electric_dipole_velocity @property def gs_moment(self) -> np.ndarray: @@ -436,20 +466,27 @@ 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.electric_quadrupole(self._gauge_origin) @property def gs_moment(self) -> np.ndarray: - nuclear_contribution = self._state.reference_state.nuclear_quadrupole(self._gauge_origin) - # Note: This next line is needed because of missing minus sign in the definition - # of the operators, TODO: remove if adcc/PR#190 is merged - ints = [-1.0 * ints for ints in self.integrals] - return compute_ground_state_moment(self._state, ints, self._pm_level, + 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(self) -> np.ndarray: - return self._state.transition_dipole_moment_velocity + if isinstance(self._state, MockExcitedStates): + return self._state.transition_quadrupole_moment + else: + return self._state.transition_quadrupole_moment(self._gauge_origin) def _state_to_state_transition_moment(self) -> np.ndarray: if isinstance(self._state, MockExcitedStates): @@ -465,8 +502,11 @@ def _operator(self) -> Operator: @property def integrals(self) -> list[adcc.OneParticleOperator]: - return self._state.reference_state.operators.\ - electric_quadrupole_velocity(self._gauge_origin) + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError + else: + return self._state.reference_state.operators.\ + electric_quadrupole_velocity(self._gauge_origin) @property def gs_moment(self) -> np.ndarray: @@ -474,7 +514,10 @@ def gs_moment(self) -> np.ndarray: "quadrupole operator in velocity gauge.") def _transition_moment(self) -> np.ndarray: - return self._state.transition_quadrupole_moment_velocity(self._gauge_origin) + if isinstance(self._state, MockExcitedStates): + return self._state.transition_quadrupole_moment_velocity + else: + return self._state.transition_quadrupole_moment_velocity(self._gauge_origin) def _state_to_state_transition_moment(self) -> np.ndarray: if isinstance(self._state, MockExcitedStates): @@ -490,14 +533,24 @@ def _operator(self) -> Operator: @property def integrals(self) -> list[adcc.OneParticleOperator]: - return self._state.reference_state.operators.diamagentic_magnetizability(self._gauge_origin) + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError + else: + return self._state.reference_state.operators.\ + diamagentic_magnetizability(self._gauge_origin) @property def gs_moment(self) -> np.ndarray: - return compute_ground_state_moment(self._state, self.integrals, self._pm_level) + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError + else: + return compute_ground_state_moment(self._state, self.integrals, self._pm_level) def _transition_moment(self) -> np.ndarray: - return compute_transition_moments(self._state, self.integrals) + if isinstance(self._state, MockExcitedStates): + raise NotImplementedError + else: + return compute_transition_moments(self._state, self.integrals) def _state_to_state_transition_moment(self) -> np.ndarray: if isinstance(self._state, MockExcitedStates): diff --git a/responsefun/SumOverStates.py b/responsefun/SumOverStates.py index 3e9a127..4db6e4c 100644 --- a/responsefun/SumOverStates.py +++ b/responsefun/SumOverStates.py @@ -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..637ac3f 100644 --- a/responsefun/evaluate_property.py +++ b/responsefun/evaluate_property.py @@ -309,7 +309,7 @@ def determine_rvecs(rvecs_dict_list, input_subs, adcc_prop, rhs = adcop.modified_transition_moments() if key[3] == 0.0: for c in components: - # list indices must be integers (1-D operators) + # # list indices must be integers (1-D operators) c = c[0] if len(c) == 1 else c response[c] = solve_response( matrix, rhs[c], -key[2], gamma=0.0, projection=projection, **solver_args @@ -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. @@ -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. @@ -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. @@ -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 = [] @@ -1312,6 +1315,7 @@ def evaluate_property_sos_fast( f"Created string of subscript labels that is used by np.einsum for term {it+1}:\n", einsum_string, ) + print(array_list[-1].shape) res_tens += factor * np.einsum(einsum_string, *array_list) res_tens = process_complex_factor(sos, res_tens) diff --git a/responsefun/test_twodimensional_operators.py b/responsefun/test_twodimensional_operators.py new file mode 100644 index 0000000..e93084f --- /dev/null +++ b/responsefun/test_twodimensional_operators.py @@ -0,0 +1,116 @@ +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, + n, + op_a, + op_b, + op_c, + op_d, + opq_ab, + opq_bc, + opq_cd, + opq_de, + opq_ef, + 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, opq_ab, n) * TransitionMoment(n, op_c, O) / (w_n - w) + + TransitionMoment(O, op_c, n) * TransitionMoment(n, opq_ab, O) / (w_n + w) + ), + "bc": ( + TransitionMoment(O, op_a, n) * TransitionMoment(n, opq_bc, O) / (w_n - w) + + TransitionMoment(O, opq_bc, n) * TransitionMoment(n, op_a, O) / (w_n + w) + ), + "abcd": ( + TransitionMoment(O, opq_ab, n) * TransitionMoment(n, opq_cd, O) / (w_n - w) + + TransitionMoment(O, opq_cd, n) * TransitionMoment(n, opq_ab, O) / (w_n + w) + ), +} + + +SOS_beta_like = { + "ab": ( + TransitionMoment(O, opq_ab, n) + * TransitionMoment(n, op_b, k) + * TransitionMoment(k, op_c, O) + / ((w_n - w_o) * (w_k - w_2)) + ), + "cd": ( + TransitionMoment(O, op_a, n) + * TransitionMoment(n, op_b, k) + * TransitionMoment(k, opq_cd, O) + / ((w_n - w_o) * (w_k - w_2)) + ), + "abde": ( + TransitionMoment(O, opq_ab, n) + * TransitionMoment(n, op_c, k) + * TransitionMoment(k, opq_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/dump_full_diagonalization.py b/responsefun/testdata/dump_full_diagonalization.py index a961db2..b1e6d56 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("_") @@ -24,7 +25,7 @@ def main(): ) state = adcc.run_adc(method=method, data_or_matrix=scfres, n_singlets=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 +56,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): From 2b6b950461b76b976ecf3796e86a13f019c2a596 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Tue, 18 Mar 2025 17:56:49 +0100 Subject: [PATCH 03/13] rename operators --- examples/alpha.py | 6 +- examples/beta_complex.py | 20 +- examples/beta_sos.py | 12 +- examples/esp.py | 6 +- examples/mcd.py | 8 +- examples/rixs.py | 8 +- examples/second_order_nlo.py | 12 +- examples/third_order_nlo.py | 18 +- examples/threepa.py | 12 +- examples/tpa.py | 8 +- responsefun/SumOverStates.py | 4 +- responsefun/evaluate_property.py | 6 +- responsefun/symbols_and_labels.py | 60 +++--- responsefun/test_extra_terms.py | 22 +- responsefun/test_magnetic_properties.py | 200 +++++++++---------- responsefun/test_property.py | 72 +++---- responsefun/test_twodimensional_operators.py | 48 ++--- 17 files changed, 261 insertions(+), 261 deletions(-) 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..330564d 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,19 @@ # 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( 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..aa60130 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 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..72b0a53 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,8 +56,8 @@ 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), @@ -69,7 +69,7 @@ def compute_beta_parallel(beta_tens, dip_mom): # TODO: remove minus after adc-connect/adcc#190 is merged beta_tens = -1.0 * 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..c565e30 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,8 +45,8 @@ 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): @@ -56,7 +56,7 @@ def threepa_average(tens): # TODO: remove minus after adc-connect/adcc#190 is merged threepa_tens = -1.0 * 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/responsefun/SumOverStates.py b/responsefun/SumOverStates.py index 4db6e4c..9480e02 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): (, ), - 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_, w_2)]. Returns ---------- @@ -325,7 +325,7 @@ def __init__( List of (op, freq) pairs whose permutation yields the full SOS expression; (op, freq): (, ), - 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_, w_2)]. excluded_states: list of or int, optional List of states that are excluded from the summation. diff --git a/responsefun/evaluate_property.py b/responsefun/evaluate_property.py index 637ac3f..b4666d0 100644 --- a/responsefun/evaluate_property.py +++ b/responsefun/evaluate_property.py @@ -498,7 +498,7 @@ def evaluate_property_isr( List of (op, freq) pairs whose permutation yields the full SOS expression; (op, freq): (, ), - 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 or int, optional List of states that are excluded from the summation. @@ -798,7 +798,7 @@ def evaluate_property_sos( List of (op, freq) pairs whose permutation yields the full SOS expression; (op, freq): (, ), - 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 or int, optional List of states that are excluded from the summation. @@ -1053,7 +1053,7 @@ def evaluate_property_sos_fast( List of (op, freq) pairs whose permutation yields the full SOS expression; (op, freq): (, ), - 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 or int, optional List of states that are excluded from the summation. diff --git a/responsefun/symbols_and_labels.py b/responsefun/symbols_and_labels.py index a583337..609810e 100644 --- a/responsefun/symbols_and_labels.py +++ b/responsefun/symbols_and_labels.py @@ -24,43 +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 -op_vel_a = OneParticleOperator("A", "electric_dipole_velocity", False) -op_vel_b = OneParticleOperator("B", "electric_dipole_velocity", False) -op_vel_c = OneParticleOperator("C", "electric_dipole_velocity", False) -op_vel_d = OneParticleOperator("D", "electric_dipole_velocity", False) -op_vel_e = OneParticleOperator("E", "electric_dipole_velocity", False) +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) +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 -opq_ab = OneParticleOperator("AB", "electric_quadrupole", False) -opq_bc = OneParticleOperator("BC", "electric_quadrupole", False) -opq_cd = OneParticleOperator("CD", "electric_quadrupole", False) -opq_de = OneParticleOperator("DE", "electric_quadrupole", False) -opq_ef = OneParticleOperator("EF", "electric_quadrupole", False) +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 -opq_vel_ab = OneParticleOperator("AB", "electric_quadrupole_velocity", False) -opq_vel_bc = OneParticleOperator("BC", "electric_quadrupole_velocity", False) -opq_vel_cd = OneParticleOperator("CD", "electric_quadrupole_velocity", False) -opq_vel_de = OneParticleOperator("DE", "electric_quadrupole_velocity", False) -opq_vel_ef = OneParticleOperator("EF", "electric_quadrupole_velocity", False) +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 operator -opxi_ab = OneParticleOperator("AB", "diamagnetic_magnetizability", False) -opxi_bc = OneParticleOperator("BC", "diamagnetic_magnetizability", False) -opxi_cd = OneParticleOperator("CD", "diamagnetic_magnetizability", False) -opxi_de = OneParticleOperator("DE", "diamagnetic_magnetizability", False) -opxi_ef = OneParticleOperator("EF", "diamagnetic_magnetizability", False) \ No newline at end of file +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..64279e4 100644 --- a/responsefun/test_magnetic_properties.py +++ b/responsefun/test_magnetic_properties.py @@ -12,16 +12,16 @@ k, m, n, - op_a, - op_b, - op_c, - op_d, - op_e, - opm_a, - opm_b, - opm_c, - opm_d, - opm_e, + mu_a, + mu_b, + mu_c, + mu_d, + mu_e, + m_a, + m_b, + m_c, + m_d, + m_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..0272e36 100644 --- a/responsefun/test_property.py +++ b/responsefun/test_property.py @@ -25,10 +25,10 @@ k, m, n, - op_a, - op_b, - op_c, - op_d, + mu_a, + mu_b, + mu_c, + mu_d, p, w, w_prime, @@ -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)], ), } diff --git a/responsefun/test_twodimensional_operators.py b/responsefun/test_twodimensional_operators.py index e93084f..f2e5e9e 100644 --- a/responsefun/test_twodimensional_operators.py +++ b/responsefun/test_twodimensional_operators.py @@ -11,15 +11,15 @@ O, k, n, - op_a, - op_b, - op_c, - op_d, - opq_ab, - opq_bc, - opq_cd, - opq_de, - opq_ef, + mu_a, + mu_b, + mu_c, + mu_d, + q_ab, + q_bc, + q_cd, + q_de, + q_ef, w, w_1, w_2, @@ -41,37 +41,37 @@ def run_scf(molecule, basis, backend="pyscf"): SOS_alpha_like = { "ab": ( - TransitionMoment(O, opq_ab, n) * TransitionMoment(n, op_c, O) / (w_n - w) - + TransitionMoment(O, op_c, n) * TransitionMoment(n, opq_ab, O) / (w_n + w) + 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, op_a, n) * TransitionMoment(n, opq_bc, O) / (w_n - w) - + TransitionMoment(O, opq_bc, n) * TransitionMoment(n, op_a, O) / (w_n + w) + 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, opq_ab, n) * TransitionMoment(n, opq_cd, O) / (w_n - w) - + TransitionMoment(O, opq_cd, n) * TransitionMoment(n, opq_ab, O) / (w_n + w) + 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, opq_ab, n) - * TransitionMoment(n, op_b, k) - * TransitionMoment(k, op_c, O) + 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, op_a, n) - * TransitionMoment(n, op_b, k) - * TransitionMoment(k, opq_cd, O) + 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, opq_ab, n) - * TransitionMoment(n, op_c, k) - * TransitionMoment(k, opq_de, O) + TransitionMoment(O, q_ab, n) + * TransitionMoment(n, mu_c, k) + * TransitionMoment(k, q_de, O) / ((w_n - w_o) * (w_k - w_2)) ), } From 3135c78f5aea0c695382b8a7c5af0dae95e092be Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Fri, 21 Mar 2025 11:26:37 +0100 Subject: [PATCH 04/13] changes to adapt to new adcc version --- responsefun/testdata/dump_full_diagonalization.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/responsefun/testdata/dump_full_diagonalization.py b/responsefun/testdata/dump_full_diagonalization.py index b1e6d56..f3f67ff 100644 --- a/responsefun/testdata/dump_full_diagonalization.py +++ b/responsefun/testdata/dump_full_diagonalization.py @@ -23,7 +23,8 @@ 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(gauge_origin) From 3ab92e0a07422a846c9c26cf8cb39b2532a1ec34 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Fri, 21 Mar 2025 15:25:50 +0100 Subject: [PATCH 05/13] some typos, tolerance criteria --- responsefun/AdccProperties.py | 2 +- responsefun/test_property.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 5be708c..0810d98 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -353,7 +353,7 @@ 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) diff --git a/responsefun/test_property.py b/responsefun/test_property.py index 0272e36..feb3739 100644 --- a/responsefun/test_property.py +++ b/responsefun/test_property.py @@ -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("_") From 1a320dc0944894c2575478a7a2577b991023c31b Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Fri, 21 Mar 2025 16:05:11 +0100 Subject: [PATCH 06/13] pmlevel in gs moments --- responsefun/AdccProperties.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 0810d98..1f43815 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -177,7 +177,7 @@ def compute_ground_state_moment(state, integrals, pm_level, for c in components: ref_mom[c] = product_trace(integrals[c], state.ground_state.reference_state.density) - if pm_level == 1: + if pm_level in [0, 1]: gs_mom = ref_mom elif pm_level == 2: mp2corr = np.zeros(op_shape) From 1e293999b1b040ed731ad37102ddb61c62334d23 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Fri, 21 Mar 2025 17:09:34 +0100 Subject: [PATCH 07/13] code style & is_imag Adcc Properties --- responsefun/AdccProperties.py | 9 ++++++--- responsefun/test_magnetic_properties.py | 12 ++++++------ responsefun/test_property.py | 4 ++-- responsefun/test_twodimensional_operators.py | 5 ++--- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 1f43815..2ff21cd 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -68,7 +68,8 @@ class Operator: name="electric_dipole_velocity", symbol="mu_p", symmetry=Symmetry.ANTIHERMITIAN, - dim=1 + dim=1, + is_imag=True, ), Operator( name="magnetic_dipole", @@ -81,13 +82,15 @@ class Operator: name="electric_quadrupole", symbol="q", symmetry=Symmetry.HERMITIAN, - dim=2 + dim=2, + is_imag=False, ), Operator( name="electric_quadrupole_velocity", symbol="q_p", symmetry=Symmetry.ANTIHERMITIAN, - dim=2 + dim=2, + is_imag=True, ), Operator( name="diamagnetic_magnetizability", diff --git a/responsefun/test_magnetic_properties.py b/responsefun/test_magnetic_properties.py index 64279e4..6b9ed34 100644 --- a/responsefun/test_magnetic_properties.py +++ b/responsefun/test_magnetic_properties.py @@ -11,17 +11,17 @@ O, k, m, - n, - mu_a, - mu_b, - mu_c, - mu_d, - mu_e, m_a, m_b, m_c, m_d, m_e, + mu_a, + mu_b, + mu_c, + mu_d, + mu_e, + n, p, w, w_1, diff --git a/responsefun/test_property.py b/responsefun/test_property.py index feb3739..2ad7b2d 100644 --- a/responsefun/test_property.py +++ b/responsefun/test_property.py @@ -24,14 +24,13 @@ gamma, k, m, - n, mu_a, mu_b, mu_c, mu_d, + n, 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 diff --git a/responsefun/test_twodimensional_operators.py b/responsefun/test_twodimensional_operators.py index f2e5e9e..9d2e701 100644 --- a/responsefun/test_twodimensional_operators.py +++ b/responsefun/test_twodimensional_operators.py @@ -10,16 +10,14 @@ from responsefun.symbols_and_labels import ( O, k, - n, mu_a, mu_b, mu_c, - mu_d, + n, q_ab, q_bc, q_cd, q_de, - q_ef, w, w_1, w_2, @@ -30,6 +28,7 @@ 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, From d9800400ce3e120cff5b02f5bb02ae013e602c87 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Fri, 28 Mar 2025 13:52:31 +0100 Subject: [PATCH 08/13] typo and pin new adcc version --- ci_env.yml | 4 ++-- pyproject.toml | 2 +- responsefun/AdccProperties.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ci_env.yml b/ci_env.yml index 8dfd54c..22d198f 100644 --- a/ci_env.yml +++ b/ci_env.yml @@ -3,10 +3,10 @@ channels: - conda-forge dependencies: - libblas =*=*mkl - - adcc >=0.15.16 + - 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/pyproject.toml b/pyproject.toml index 2e3f5a7..a126f4d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,7 +21,7 @@ classifiers = [ requires-python = ">=3.7" # Declare any run-time dependencies that should be installed with the package. dependencies = [ - "adcc>=0.15.16", + "adcc>=0.16.0", "respondo>=0.0.5", "numpy>=1.14", "sympy", diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 2ff21cd..e53280d 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -540,7 +540,7 @@ def integrals(self) -> list[adcc.OneParticleOperator]: raise NotImplementedError else: return self._state.reference_state.operators.\ - diamagentic_magnetizability(self._gauge_origin) + diamagnetic_magnetizability(self._gauge_origin) @property def gs_moment(self) -> np.ndarray: From 066a511e3dd865b191a4224c44f6bcc26356374a Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Mon, 31 Mar 2025 09:36:08 +0200 Subject: [PATCH 09/13] remove minus sign (adcc/PR#190) --- examples/beta_complex.py | 5 +---- examples/mcd.py | 5 +---- examples/second_order_nlo.py | 5 +---- examples/threepa.py | 5 +---- 4 files changed, 4 insertions(+), 16 deletions(-) diff --git a/examples/beta_complex.py b/examples/beta_complex.py index 330564d..65c240c 100644 --- a/examples/beta_complex.py +++ b/examples/beta_complex.py @@ -45,10 +45,7 @@ + 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=[(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), diff --git a/examples/mcd.py b/examples/mcd.py index aa60130..7509dd4 100644 --- a/examples/mcd.py +++ b/examples/mcd.py @@ -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/second_order_nlo.py b/examples/second_order_nlo.py index 72b0a53..7286d2c 100644 --- a/examples/second_order_nlo.py +++ b/examples/second_order_nlo.py @@ -64,10 +64,7 @@ def compute_beta_parallel(beta_tens, dip_mom): "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=[(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)], freqs_in=[(w_1, freqs[0]), (w_2, freqs[1])], diff --git a/examples/threepa.py b/examples/threepa.py index c565e30..e6d7730 100644 --- a/examples/threepa.py +++ b/examples/threepa.py @@ -51,10 +51,7 @@ def threepa_average(tens): 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=[(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)], From fd5db9ea75a35c5e83c69ccd4e13810fc92681a1 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Sun, 25 May 2025 11:46:08 +0200 Subject: [PATCH 10/13] update testdata --- responsefun/testdata/0_download_testdata.sh | 4 ++-- responsefun/testdata/SHA256SUMS | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) 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 From 1f6de39c5cdbb819d3066d5441c202f13246f392 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Sun, 25 May 2025 11:54:36 +0200 Subject: [PATCH 11/13] update README.md --- README.md | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) 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 From 26426cfe222d8a81ad9797fde92f04704794c0c9 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Mon, 26 May 2025 02:46:10 +0200 Subject: [PATCH 12/13] pin sympy version, cleanup --- pyproject.toml | 2 +- responsefun/AdccProperties.py | 3 --- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index dc17854..36de855 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ dependencies = [ "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 e53280d..94cdf98 100644 --- a/responsefun/AdccProperties.py +++ b/responsefun/AdccProperties.py @@ -168,9 +168,6 @@ def compute_ground_state_moment(state, integrals, pm_level, """ Computes the ground state moments. """ - # Note that the charge (see adcc/PR#190) must already be contained in the integrals - # passed to this function. - # convert integrals to np.array (tuples are not allowed as indices for lists) integrals = np.array(integrals) op_shape = np.shape(integrals) From 75a5f425c80b0b466ba5a1b9c0e42466427d1354 Mon Sep 17 00:00:00 2001 From: Friederike Schneider Date: Tue, 27 May 2025 10:52:08 +0200 Subject: [PATCH 13/13] github comments --- responsefun/AdccProperties.py | 123 ++++++++++++------------------ responsefun/SumOverStates.py | 4 +- responsefun/evaluate_property.py | 3 +- responsefun/symbols_and_labels.py | 2 +- 4 files changed, 53 insertions(+), 79 deletions(-) diff --git a/responsefun/AdccProperties.py b/responsefun/AdccProperties.py index 94cdf98..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 . # - import warnings from abc import ABC, abstractmethod, abstractproperty from dataclasses import dataclass @@ -210,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) @@ -247,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: @@ -255,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: @@ -285,35 +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 + self, comp: Union[int, tuple[int], None] = None ) -> Union[adcc.AmplitudeVector, np.ndarray[adcc.AmplitudeVector]]: if comp is None: - op = self.integrals - else: - op = np.array(self.integrals)[comp] - if isinstance(op, adcc.OneParticleOperator): - mtms = modified_transition_moments( - self._property_method, self._state.ground_state, op - ) - elif isinstance(op, list): mtms = np.array([(modified_transition_moments(self._property_method, self._state.ground_state, op_ints)) for op_ints in self.integrals]) + else: + 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: @@ -324,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 @@ -390,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): @@ -419,17 +422,11 @@ def gs_moment(self) -> np.ndarray: else: return compute_ground_state_moment(self._state, self.integrals, self._pm_level) - def _transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - return self._state.transition_magnetic_dipole_moment - else: - return self._state.transition_magnetic_dipole_moment(self._gauge_origin) + def _transition_moment_mock(self) -> np.ndarray: + return self._state.transition_magnetic_dipole_moment - def _state_to_state_transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - return self._state.transition_magnetic_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_magnetic_moment_s2s class ElectricDipoleVelocity(AdccProperties): @@ -449,14 +446,11 @@ def gs_moment(self) -> np.ndarray: raise NotImplementedError("Ground-state moments not implemented for electric dipole " "operator in velocity gauge.") - def _transition_moment(self) -> np.ndarray: + def _transition_moment_mock(self) -> np.ndarray: return self._state.transition_dipole_moment_velocity - def _state_to_state_transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - raise NotImplementedError() - else: - return compute_state_to_state_transition_moments(self._state, self.integrals) + def _state_to_state_transition_moment_mock(self) -> np.ndarray: + raise NotImplementedError class ElectricQuadrupole(AdccProperties): @@ -465,7 +459,7 @@ def _operator(self) -> Operator: return get_operator_by_name("electric_quadrupole") @property - def integrals(self) -> list[adcc.OneParticleOperator]: + def integrals(self) -> list[list[adcc.OneParticleOperator]]: if isinstance(self._state, MockExcitedStates): raise NotImplementedError else: @@ -481,18 +475,11 @@ def gs_moment(self) -> np.ndarray: 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 _transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - return self._state.transition_quadrupole_moment - else: - return self._state.transition_quadrupole_moment(self._gauge_origin) - - def _state_to_state_transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - raise NotImplementedError() - else: - return compute_state_to_state_transition_moments(self._state, self.integrals) + def _state_to_state_transition_moment_mock(self) -> np.ndarray: + raise NotImplementedError class ElectricQuadrupoleVelocity(AdccProperties): @@ -501,7 +488,7 @@ def _operator(self) -> Operator: return get_operator_by_name("electric_quadrupole_velocity") @property - def integrals(self) -> list[adcc.OneParticleOperator]: + def integrals(self) -> list[list[adcc.OneParticleOperator]]: if isinstance(self._state, MockExcitedStates): raise NotImplementedError else: @@ -513,17 +500,11 @@ def gs_moment(self) -> np.ndarray: raise NotImplementedError("Ground-state moments not implemented for electric " "quadrupole operator in velocity gauge.") - def _transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - return self._state.transition_quadrupole_moment_velocity - else: - return self._state.transition_quadrupole_moment_velocity(self._gauge_origin) + def _transition_moment_mock(self) -> np.ndarray: + return self._state.transition_quadrupole_moment_velocity - def _state_to_state_transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - raise NotImplementedError() - else: - return compute_state_to_state_transition_moments(self._state, self.integrals) + def _state_to_state_transition_moment_mock(self) -> np.ndarray: + raise NotImplementedError class DiamagneticMagnetizability(AdccProperties): @@ -532,7 +513,7 @@ def _operator(self) -> Operator: return get_operator_by_name("diamagentic_magnetizability") @property - def integrals(self) -> list[adcc.OneParticleOperator]: + def integrals(self) -> list[list[adcc.OneParticleOperator]]: if isinstance(self._state, MockExcitedStates): raise NotImplementedError else: @@ -546,14 +527,8 @@ def gs_moment(self) -> np.ndarray: else: return compute_ground_state_moment(self._state, self.integrals, self._pm_level) - def _transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - raise NotImplementedError - else: - return compute_transition_moments(self._state, self.integrals) + def _transition_moment_mock(self) -> np.ndarray: + raise NotImplementedError - def _state_to_state_transition_moment(self) -> np.ndarray: - if isinstance(self._state, MockExcitedStates): - raise NotImplementedError() - else: - return compute_state_to_state_transition_moments(self._state, self.integrals) \ No newline at end of file + 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 9480e02..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): (, ), - e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_, 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): (, ), - e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_, w_2)]. + e.g., [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)]. excluded_states: list of or int, optional List of states that are excluded from the summation. diff --git a/responsefun/evaluate_property.py b/responsefun/evaluate_property.py index b4666d0..3067a74 100644 --- a/responsefun/evaluate_property.py +++ b/responsefun/evaluate_property.py @@ -309,7 +309,7 @@ def determine_rvecs(rvecs_dict_list, input_subs, adcc_prop, rhs = adcop.modified_transition_moments() if key[3] == 0.0: for c in components: - # # list indices must be integers (1-D operators) + # list indices must be integers (1-D operators) c = c[0] if len(c) == 1 else c response[c] = solve_response( matrix, rhs[c], -key[2], gamma=0.0, projection=projection, **solver_args @@ -1315,7 +1315,6 @@ def evaluate_property_sos_fast( f"Created string of subscript labels that is used by np.einsum for term {it+1}:\n", einsum_string, ) - print(array_list[-1].shape) res_tens += factor * np.einsum(einsum_string, *array_list) res_tens = process_complex_factor(sos, res_tens) diff --git a/responsefun/symbols_and_labels.py b/responsefun/symbols_and_labels.py index 609810e..97b500f 100644 --- a/responsefun/symbols_and_labels.py +++ b/responsefun/symbols_and_labels.py @@ -58,7 +58,7 @@ qp_de = OneParticleOperator("DE", "electric_quadrupole_velocity", False) qp_ef = OneParticleOperator("EF", "electric_quadrupole_velocity", False) -# diamagnetic magnetizability operator +# diamagnetic magnetizability operators xi_ab = OneParticleOperator("AB", "diamagnetic_magnetizability", False) xi_bc = OneParticleOperator("BC", "diamagnetic_magnetizability", False) xi_cd = OneParticleOperator("CD", "diamagnetic_magnetizability", False)