Davidson works when doubles folding in ADC(2)

adcc/AdcMatrix.py

@@ -30,7 +30,7 @@
 from .functions import ones_like
 from .Intermediates import Intermediates
 from .AmplitudeVector import AmplitudeVector
-
+# from .solver.explicit_symmetrisation import IndexSymmetrisation
 
 class AdcExtraTerm:
     def __init__(self, matrix, blocks):
@@ -650,3 +650,56 @@ def block_view(self, block):
                                   "projected ADC matrices.")
         # TODO The way to implement this is to ask the inner matrix to
         #      a block_view and then wrap that in an AdcMatrixProjected.
+
+
+class AdcMatrixFolded(AdcMatrix):
+    def __init__(self, matrix):
+        """
+        Initialise a ADC(2) matrix when using doubles-folding.
+
+        Parameters
+        ----------
+        matrix : AdcMatrix
+           ADC(2) matrix
+        """
+        super().__init__(matrix.method, matrix.ground_state,
+                         block_orders=matrix.block_orders,
+                         intermediates=matrix.intermediates)
+        if matrix.method.name != "adc2":
+            raise TypeError("The matrix should be ADC2 Matrix for doubles-folding.")
+
+        from .solver.explicit_symmetrisation import IndexSymmetrisation
+        self.isymm = IndexSymmetrisation(matrix)
+
+    def diagonal(self):
+        """
+        Return the approximate diagonal of the ADC(2) matrix with doubles-folding.
+        """
+        return AmplitudeVector(ph=super().diagonal().ph)
+
+    def u2_fold(self, other):
+        """
+        Return the doubles component of u in terms of its singles component.
+        """
+        diag = super().diagonal().pphh
+        e = diag.ones_like()
+        u2 = self.block_apply("pphh_ph", other.ph) / (e * self.omega - diag)
+        return self.isymm.symmetrise(AmplitudeVector(pphh=u2))
+
+    def matvec(self, other):
+        """
+        Compute the doubles-folded matrix-vector product of the singles vector with
+        an effective ADC matrix which depends on the eigenvalue(omega).
+        """
+        u2 = self.u2_fold(other)
+        return AmplitudeVector(ph=self.block_apply("ph_ph", other.ph)
+                               + self.block_apply("ph_pphh", u2.pphh))
+
+    def unfold(self, other):
+        """
+        Recompute the doubles component and return the complete vector.
+        """
+        other_symm = self.isymm.symmetrise(other)
+        u2 = self.u2_fold(other_symm)
+
+        return AmplitudeVector(ph=other_symm.ph, pphh=u2.pphh)