ADC Nuclear Gradients

adcc/ElectronicStates.py

@@ -160,6 +160,19 @@ def excitation_energy(self):
         """Excitation energies including all corrections in atomic units"""
         return self._excitation_energy
 
+    @property
+    def total_energy(self):
+        """Total energies (ground state + excitation) of all computed states"""
+        return np.array([self._total_energy(i) for i in range(self.size)])
+
+    def _total_energy(self, state_n: int) -> float:
+        """Computes the total energy for a single state"""
+        if self.method.level == MethodLevel.ZERO:
+            return (self._excitation_energy[state_n]
+                    + self.reference_state.energy_scf)
+        return (self._excitation_energy[state_n]
+                + self.ground_state.energy(self.method.level.to_int()))
+
     @property
     def excitation_energy_uncorrected(self):
         """Excitation energies without any corrections in atomic units"""

adcc/Excitation.py

@@ -94,6 +94,21 @@ def transition_quadrupole_moment_velocity(self,
             state_n=self.index, gauge_origin=gauge_origin
         )
 
+    @property
+    def ground_state(self):
+        """The ground state (LazyMp) of the parent ElectronicStates"""
+        return self._parent_state.ground_state
+
+    @property
+    def reference_state(self):
+        """The HF reference state of the parent ElectronicStates"""
+        return self._parent_state.reference_state
+
+    @property
+    def method(self):
+        """The ADC method of the parent ElectronicStates"""
+        return self._parent_state.method
+
     @property
     def oscillator_strength(self) -> np.float64:
         """The oscillator strength"""

adcc/ReferenceState.py

@@ -159,7 +159,9 @@ def __init__(self, hfdata, core_orbitals=None, frozen_core=None,
         )
 
         self.environment = None  # no environment attached by default
-        for name in ["excitation_energy_corrections", "environment"]:
+        for name in ["excitation_energy_corrections",
+                     "environment",
+                     "gradient_provider"]:
             if hasattr(hfdata, name):
                 setattr(self, name, getattr(hfdata, name))
 

adcc/StateView.py

@@ -39,6 +39,11 @@ def excitation_energy(self):
         """Excitation energy including all corrections in atomic units"""
         return self._parent_state.excitation_energy[self.index]
 
+    @property
+    def total_energy(self) -> float:
+        """Total energy (ground state + excitation) of this state"""
+        return self._parent_state._total_energy(self.index)
+
     @property
     def excitation_energy_uncorrected(self):
         """Excitation energy without any corrections in atomic units"""

adcc/__init__.py

@@ -47,6 +47,7 @@
 from .TwoParticleOperator import TwoParticleOperator
 from .TwoParticleDensity import TwoParticleDensity
 from .opt_einsum_integration import register_with_opt_einsum
+from .gradients import NuclearGradientScanner, nuclear_gradient
 
 # This has to be the last set of import
 from .guess import (guess_symmetries, guess_zero, guesses_any, guesses_singlet,
@@ -68,6 +69,7 @@
            "guess_symmetries", "guesses_spin_flip", "guess_zero", "LazyMp",
            "adc0", "cis", "adc1", "adc2", "adc2x", "adc3",
            "cvs_adc0", "cvs_adc1", "cvs_adc2", "cvs_adc2x", "cvs_adc3",
+           "nuclear_gradient", "NuclearGradientScanner",
            "banner"]
 
 __version__ = "0.18.0"

adcc/backends/psi4.py

@@ -24,6 +24,8 @@
 
 from libadcc import HartreeFockProvider
 
+from adcc.gradients import GradientComponents
+
 import psi4
 
 from .EriBuilder import EriBuilder
@@ -32,6 +34,72 @@
 from ..OneParticleOperator import OneParticleOperator
 
 
+class Psi4GradientProvider:
+    def __init__(self, wfn):
+        self.wfn = wfn
+        self.mol = self.wfn.molecule()
+        self.backend = "psi4"
+        self.mints = psi4.core.MintsHelper(self.wfn)
+
+    def correlated_gradient(self, g1_ao, w_ao, g2_ao_1, g2_ao_2):
+        """
+        g1_ao: relaxed one-particle density matrix
+        w_ao:  energy-weighted density matrix
+        g2_ao_1, g2_ao_2: relaxed two-particle density matrices
+        """
+        natoms = self.mol.natom()
+        Gradient = {}
+        Gradient["N"] = np.zeros((natoms, 3))
+        Gradient["S"] = np.zeros((natoms, 3))
+        Gradient["T"] = np.zeros((natoms, 3))
+        Gradient["V"] = np.zeros((natoms, 3))
+        Gradient["OEI"] = np.zeros((natoms, 3))
+        Gradient["ERI"] = np.zeros((natoms, 3))
+        Gradient["Total"] = np.zeros((natoms, 3))
+
+        # 1st Derivative of Nuclear Repulsion
+        Gradient["N"] = psi4.core.Matrix.to_array(
+            self.mol.nuclear_repulsion_energy_deriv1([0, 0, 0])
+        )
+        # Build Integral Derivatives
+        cart = ['_X', '_Y', '_Z']
+        oei_dict = {"S": "OVERLAP", "T": "KINETIC", "V": "POTENTIAL"}
+
+        deriv1_mat = {}
+        deriv1_np = {}
+        # 1st Derivative of OEIs
+        for atom in range(natoms):
+            for key in oei_dict:
+                string = key + str(atom)
+                deriv1_mat[string] = self.mints.ao_oei_deriv1(oei_dict[key], atom)
+                for p in range(3):
+                    map_key = string + cart[p]
+                    deriv1_np[map_key] = np.asarray(deriv1_mat[string][p])
+                    if key == "S":
+                        Gradient["S"][atom, p] = np.sum(w_ao * deriv1_np[map_key])
+                    else:
+                        Gradient[key][atom, p] = np.sum(g1_ao * deriv1_np[map_key])
+        # Build Total OEI Gradient
+        Gradient["OEI"] = Gradient["T"] + Gradient["V"] + Gradient["S"]
+
+        # Build TE contributions
+        for atom in range(natoms):
+            deriv2 = self.mints.ao_tei_deriv1(atom)
+            for p in range(3):
+                deriv2_np = np.asarray(deriv2[p])
+                Gradient["ERI"][atom, p] += np.einsum(
+                    'pqrs,prqs->', g2_ao_1, deriv2_np, optimize=True
+                )
+                Gradient["ERI"][atom, p] -= np.einsum(
+                    'pqrs,psqr->', g2_ao_2, deriv2_np, optimize=True
+                )
+        ret = GradientComponents(
+            natoms, Gradient["N"], Gradient["S"],
+            Gradient["T"] + Gradient["V"], Gradient["ERI"]
+        )
+        return ret
+
+
 class Psi4OperatorIntegralProvider:
     available: tuple[str, ...] = (
         "overlap", "electric_dipole", "electric_dipole_velocity", "magnetic_dipole",
@@ -169,6 +237,7 @@ def __init__(self, wfn):
                                           wfn.nalpha(), wfn.nbeta(),
                                           self.restricted)
         self.operator_integral_provider = Psi4OperatorIntegralProvider(self.wfn)
+        self.gradient_provider = Psi4GradientProvider(self.wfn)
 
         self.environment = None
         self.environment_implementation = None