Additional operator integrals and selection of gauge origin

adcc/DataHfProvider.py

@@ -203,21 +203,39 @@ def __init__(self, data):
                                  + str((3, nb, nb)) + " not "
                                  + str(mmp["elec_1"].shape))
             opprov.electric_dipole = np.asarray(mmp["elec_1"])
+        if "elec_2" in mmp:
+            def get_integral_quad(gauge_origin):
+                integral_string = "elec_2"
+                if gauge_origin != 'origin':
+                    integral_string = f'{integral_string}_{gauge_origin}'
+                return np.asarray(mmp[integral_string])
+            if mmp["elec_2"].shape != (9, nb, nb):
+                raise ValueError("multipoles/elec_2 is expected to have "
+                                 "shape " + str((9, nb, nb)) + " not "
+                                 + str(mmp["elec_2"].shape))
+            opprov.electric_quadrupole = get_integral_quad
         magm = data.get("magnetic_moments", {})
         if "mag_1" in magm:
+            def get_integral_magm(gauge_origin):
+                integral_string = "mag_1"
+                if gauge_origin != 'origin':
+                    integral_string = f'{integral_string}_{gauge_origin}'
+                return np.asarray(magm[integral_string])
             if magm["mag_1"].shape != (3, nb, nb):
                 raise ValueError("magnetic_moments/mag_1 is expected to have "
                                  "shape " + str((3, nb, nb)) + " not "
                                  + str(magm["mag_1"].shape))
-            opprov.magnetic_dipole = np.asarray(magm["mag_1"])
+            opprov.magnetic_dipole = get_integral_magm
         derivs = data.get("derivatives", {})
         if "nabla" in derivs:
+            def get_integral_deriv(gauge_origin):
+                return np.asarray(derivs["nabla"])
             if derivs["nabla"].shape != (3, nb, nb):
                 raise ValueError("derivatives/nabla is expected to "
                                  "have shape "
                                  + str((3, nb, nb)) + " not "
                                  + str(derivs["nabla"].shape))
-            opprov.nabla = np.asarray(derivs["nabla"])
+            opprov.nabla = get_integral_deriv
         self.operator_integral_provider = opprov
 
     #
@@ -260,7 +278,7 @@ def get_backend(self):
     def get_energy_scf(self):
         return get_scalar_value(self.data, "energy_scf", 0.0)
 
-    def get_nuclear_multipole(self, order):
+    def get_nuclear_multipole(self, order, gauge_origin):
         if order == 0:  # The function interface needs an np.array on return
             nuc_0 = get_scalar_value(self.data, "multipoles/nuclear_0", 0.0)
             return np.array([nuc_0])

adcc/ElectronicTransition.py

@@ -165,35 +165,47 @@ def transition_dipole_moment(self):
             for tdm in self.transition_dm
         ])
 
-    @cached_property
     @mark_excitation_property()
     @timed_member_call(timer="_property_timer")
-    def transition_dipole_moment_velocity(self):
+    def transition_dipole_moment_velocity(self, gauge_origin=[0.0, 0.0, 0.0]):
         """List of transition dipole moments in the
         velocity gauge of all computed states"""
         if self.property_method.level == 0:
             warnings.warn("ADC(0) transition velocity dipole moments "
                           "are known to be faulty in some cases.")
-        dipole_integrals = self.operators.nabla
+        dipole_integrals = self.operators.nabla(gauge_origin)
         return np.array([
             [product_trace(comp, tdm) for comp in dipole_integrals]
             for tdm in self.transition_dm
         ])
 
-    @cached_property
     @mark_excitation_property()
     @timed_member_call(timer="_property_timer")
-    def transition_magnetic_dipole_moment(self):
+    def transition_magnetic_dipole_moment(self, gauge_origin=[0.0, 0.0, 0.0]):
         """List of transition magnetic dipole moments of all computed states"""
         if self.property_method.level == 0:
             warnings.warn("ADC(0) transition magnetic dipole moments "
                           "are known to be faulty in some cases.")
-        mag_dipole_integrals = self.operators.magnetic_dipole
+        mag_dipole_integrals = self.operators.magnetic_dipole(gauge_origin)
         return np.array([
             [product_trace(comp, tdm) for comp in mag_dipole_integrals]
             for tdm in self.transition_dm
         ])
 
+    @mark_excitation_property()
+    @timed_member_call(timer="_property_timer")
+    def transition_quadrupole_moment(self, gauge_origin=[0.0, 0.0, 0.0]):
+        """List of transition quadrupole moments of all computed states"""
+        if self.property_method.level == 0:
+            warnings.warn("ADC(0) transition quadrupole moments are known to be "
+                          "faulty in some cases.")
+        quadrupole_integrals = self.operators.electric_quadrupole(gauge_origin)
+        return np.array([
+            [[product_trace(quad1, tdm)
+                for quad1 in quad] for quad in quadrupole_integrals]
+            for tdm in self.transition_dm
+        ])
+
     @cached_property
     @mark_excitation_property()
     def oscillator_strength(self):
@@ -204,26 +216,25 @@ def oscillator_strength(self):
                                self.excitation_energy)
         ])
 
-    @cached_property
     @mark_excitation_property()
-    def oscillator_strength_velocity(self):
+    def oscillator_strength_velocity(self, gauge_origin=[0.0, 0.0, 0.0]):
         """List of oscillator strengths in
         velocity gauge of all computed states"""
         return 2. / 3. * np.array([
             np.linalg.norm(tdm)**2 / np.abs(ev)
-            for tdm, ev in zip(self.transition_dipole_moment_velocity,
+            for tdm, ev in zip(self.transition_dipole_moment_velocity(gauge_origin),
                                self.excitation_energy)
         ])
 
-    @cached_property
     @mark_excitation_property()
-    def rotatory_strength(self):
+    def rotatory_strength(self, gauge_origin=[0.0, 0.0, 0.0]):
         """List of rotatory strengths of all computed states"""
         return np.array([
             np.dot(tdm, magmom) / ee
-            for tdm, magmom, ee in zip(self.transition_dipole_moment_velocity,
-                                       self.transition_magnetic_dipole_moment,
-                                       self.excitation_energy)
+            for tdm, magmom, ee in zip(
+                self.transition_dipole_moment_velocity(gauge_origin),
+                self.transition_magnetic_dipole_moment(gauge_origin),
+                self.excitation_energy)
         ])
 
     @property

adcc/OperatorIntegrals.py

@@ -102,6 +102,9 @@ def available(self):
             "electric_dipole",
             "magnetic_dipole",
             "nabla",
+            "electric_quadrupole_traceless",
+            "electric_quadrupole",
+            "diamagnetic_magnetizability",
             "pe_induction_elec",
             "pcm_potential_elec",
         )
@@ -125,6 +128,38 @@ def import_dipole_like_operator(self, integral, is_symmetric=True):
             dipoles.append(dip_ff)
         return dipoles
 
+    def __import_gauge_dependent_dipole_like(self, callback, is_symmetric=True):
+        """Returns a function that imports a gauge-dependent dipole like
+        operator.
+        The returned function imports the operator and transforms it to the
+        molecular orbital basis.
+
+        Parameters
+        ----------
+        callback : callable
+            Function that computes the operator in atomic orbitals using
+            the gauge-origin (str or list) as single argument
+        is_symmetric : bool, optional
+            if the imported operator is symmetric, by default True
+        """
+        if not callable(callback):
+            raise TypeError
+
+        def process_operator(gauge_origin=[0.0, 0.0, 0.0], callback=callback,
+                             is_symmetric=is_symmetric):
+            dipoles = []
+            for i, component in enumerate(["x", "y", "z"]):
+                dip_backend = callback(gauge_origin)[i]
+                dip_bb = replicate_ao_block(self.mospaces, dip_backend,
+                                            is_symmetric=is_symmetric)
+                dip_ff = OneParticleOperator(self.mospaces,
+                                             is_symmetric=is_symmetric)
+                transform_operator_ao2mo(dip_bb, dip_ff, self.__coefficients,
+                                         self.__conv_tol)
+                dipoles.append(dip_ff)
+            return dipoles
+        return process_operator
+
     @property
     @timed_member_call("_import_timer")
     def electric_dipole(self):
@@ -135,18 +170,98 @@ def electric_dipole(self):
     @property
     @timed_member_call("_import_timer")
     def magnetic_dipole(self):
-        """Return the magnetic dipole integrals in the molecular orbital basis."""
-        return self.import_dipole_like_operator("magnetic_dipole",
-                                                is_symmetric=False)
+        """
+        Returns a function to obtain magnetic dipole intergrals
+        in the molecular orbital basis dependent on the selected gauge origin.
+        The default gauge origin is set to [0.0, 0.0, 0.0].
+        """
+        callback = self.provider_ao.magnetic_dipole
+        return self.__import_gauge_dependent_dipole_like(callback,
+                                                         is_symmetric=False)
 
     @property
     @timed_member_call("_import_timer")
     def nabla(self):
         """
-        Return the momentum (nabla operator) integrals
-        in the molecular orbital basis.
+        Returns a function to obtain nabla intergrals
+        in the molecular orbital basis dependent on the selected gauge origin.
+        The default gauge origin is set to [0.0, 0.0, 0.0].
+        """
+        callback = self.provider_ao.nabla
+        return self.__import_gauge_dependent_dipole_like(callback,
+                                                         is_symmetric=False)
+
+    def __import_quadrupole_like_operator(self, callback, is_symmetric=True):
+        """Returns a function that imports a gauge-dependent quadrupole like
+        operator.
+        The returned function imports the operator and transforms it to the
+        molecular orbital basis.
+
+        Parameters
+        ----------
+        callback : callable
+            Function that computes the operator in atomic orbitals using
+            a the gauge-origin (str or list) as single argument
+        is_symmetric : bool, optional
+            if the imported operator is symmetric, by default True
+        """
+        if not callable(callback):
+            raise TypeError
+
+        def process_operator(gauge_origin=[0.0, 0.0, 0.0], callback=callback,
+                             is_symmetric=is_symmetric):
+            quad = []
+            quadrupoles = []
+            for i, component in enumerate(["xx", "xy", "xz",
+                                           "yx", "yy", "yz",
+                                           "zx", "zy", "zz"]):
+                quad_backend = callback(gauge_origin)[i]
+                quad_bb = replicate_ao_block(self.mospaces, quad_backend,
+                                             is_symmetric=is_symmetric)
+                quad_ff = OneParticleOperator(self.mospaces,
+                                              is_symmetric=is_symmetric)
+                transform_operator_ao2mo(quad_bb, quad_ff,
+                                         self.__coefficients,
+                                         self.__conv_tol)
+                quad.append(quad_ff)
+            quadrupoles.append(quad[:3])
+            quadrupoles.append(quad[3:6])
+            quadrupoles.append(quad[6:])
+            return quadrupoles
+        return process_operator
+
+    @property
+    @timed_member_call("_import_timer")
+    def electric_quadrupole_traceless(self):
+        """
+        Returns a function to obtain traceless electric quadrupole integrals
+        in the molecular orbital basis dependent on the selected gauge origin.
+        The default gauge origin is set to [0.0, 0.0, 0.0].
+        """
+        callback = self.provider_ao.electric_quadrupole_traceless
+        return self.__import_quadrupole_like_operator(callback, is_symmetric=False)
+
+    @property
+    @timed_member_call("_import_timer")
+    def electric_quadrupole(self):
+        """
+        Returns a function to obtain electric quadrupole integrals
+        in the molecular orbital basis dependent on the selected gauge origin.
+        The default gauge origin is set to [0.0, 0.0, 0.0].
+        """
+        callback = self.provider_ao.electric_quadrupole
+        return self.__import_quadrupole_like_operator(callback, is_symmetric=False)
+
+    @property
+    @timed_member_call("_import_timer")
+    def diamagnetic_magnetizability(self):
+        """
+        Returns a function to obtain diamagnetic magnetizability integrals
+        in the molecular orbital basis dependent on the selected gauge origin.
+        The default gauge origin is set to [0.0, 0.0, 0.0].
         """
-        return self.import_dipole_like_operator("nabla", is_symmetric=False)
+        callback = self.provider_ao.diamagnetic_magnetizability
+        return self.__import_quadrupole_like_operator(callback, is_symmetric=False)
 
     def __import_density_dependent_operator(self, ao_callback, is_symmetric=True):
         """Returns a function that imports a density-dependent operator.

adcc/ReferenceState.py

@@ -225,4 +225,9 @@ def dipole_moment(self):
         return self.nuclear_dipole - np.array([product_trace(comp, self.density)
                                                for comp in dipole_integrals])
 
+    def nuclear_quadrupole(self, gauge_origin="origin"):
+        if isinstance(gauge_origin, str):
+            gauge_origin = self.determine_gauge_origin(gauge_origin)
+        return super().nuclear_quadrupole(np.array(gauge_origin))
+
 # TODO some nice describe method

adcc/adc_pp/modified_transition_moments.py

@@ -85,6 +85,7 @@ def mtm_cvs_adc2(mp, op, intermediates):
     "adc0": mtm_adc0,
     "adc1": mtm_adc1,
     "adc2": mtm_adc2,
+    "adc2x": mtm_adc2,  # Identical to ADC(2)
     "cvs-adc0": mtm_cvs_adc0,
     "cvs-adc1": mtm_cvs_adc0,  # Identical to CVS-ADC(0)
     "cvs-adc2": mtm_cvs_adc2,

adcc/adc_pp/test_modified_transition_moments.py

@@ -36,6 +36,7 @@
 methods = [m for bm in basemethods for m in [bm, "cvs-" + bm]]
 
 operator_kinds = ["electric", "magnetic"]
+gauge_origins = ["origin", "mass_center", "charge_center"]
 
 
 @expand_test_templates(list(itertools.product(methods, operator_kinds)))
@@ -49,8 +50,9 @@ def base_test(self, system, method, kind, op_kind):
             dips = state.reference_state.operators.electric_dipole
             ref_tdm = ref["transition_dipole_moments"]
         elif op_kind == "magnetic":
-            dips = state.reference_state.operators.magnetic_dipole
-            ref_tdm = ref["transition_magnetic_dipole_moments"]
+            for gauge_origin in gauge_origins:
+                dips = state.reference_state.operators.magnetic_dipole(gauge_origin)
+                ref_tdm = ref["transition_magnetic_dipole_moments"][gauge_origin]
         else:
             skip("Tests are only implemented for electric "
                  "and magnetic dipole operators.")

adcc/backends/molsturm.py

@@ -21,14 +21,17 @@
 ##
 ## ---------------------------------------------------------------------
 import numpy as np
-
+import warnings
 from molsturm.State import State
 from adcc.DataHfProvider import DataHfProvider
 
 
 def convert_scf_to_dict(scfres):
     if not isinstance(scfres, State):
         raise TypeError("Unsupported type for backends.molsturm.import_scf.")
+    warnings.warn("Gauge origin selection not available "
+                  "in Molsturm. "
+                  "The gauge origin is selected as [0, 0, 0].")
 
     n_alpha = scfres["n_alpha"]
     n_beta = scfres["n_beta"]