Active space approaches

FEATURES.md

@@ -16,6 +16,7 @@ The following table summarises the available methods and routines for the ansatz
 | CCSDT       |  RUG  |  RUG  |  RUG  |       |       |       |  RUG  |  RUG  |   -   |
 | CCSDTQ      |   g   |   g   |       |       |       |       |       |       |   -   |
 | CCSD(T)     |  RuG  |  RuG  |       |       |       |       |       |       |   -   |
+| CCSDt'      |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
 | CC2         |  RUG  |  RUG  |  RUG  |       |       |       |  RUG  |  RUG  |   -   |
 | CC3         |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
 | QCISD       |  RUG  |  RUG  |       |       |       |       |       |       |   -   |

ebcc/ansatz.py

@@ -18,6 +18,7 @@
     "CC2": ("CC2", "", 0, 0),
     "CC3": ("CC3", "", 0, 0),
     "QCISD": ("QCISD", "", 0, 0),
+    "CCSDt'": ("CCSDt'", "", 0, 0),
     "CCSD-S-1-1": ("CCSD", "S", 1, 1),
     "CCSD-SD-1-1": ("CCSD", "SD", 1, 1),
     "CCSD-SD-1-2": ("CCSD", "SD", 1, 2),
@@ -37,6 +38,7 @@ def name_to_identifier(name):
     iden = name.replace("(", "x").replace(")", "x")
     iden = iden.replace("[", "y").replace("]", "y")
     iden = iden.replace("-", "_")
+    iden = iden.replace("'", "p")
 
     return iden
 
@@ -55,6 +57,7 @@ def identifity_to_name(iden):
         name = name.replace("x", "(", 1).replace("x", ")", 1)
     while "y" in name:
         name = name.replace("y", "(", 1).replace("y", ")", 1)
+    name = name.replace("p", "'")
 
     return name
 
@@ -173,124 +176,137 @@ def is_one_shot(self):
             for ansatz in (self.fermion_ansatz, self.boson_ansatz)
         )
 
-    @property
-    def correlated_cluster_ranks(self):
-        """Get a list of cluster operator rank numbers for each of
-        the fermionic, bosonic, and coupling ansatzes, for the
-        correlated space (see space.py).
+    def fermionic_cluster_ranks(self, spin_type="G"):
+        """Get a list of cluster operator ranks for the fermionic
+        space.
 
         Returns
         -------
-        ranks : tuple of tuple of int
-            Cluster operator ranks for the fermionic, bosonic, and
-            coupling ansatzes, for the correlated space.
+        ranks : list of tuples
+            List of cluster operator ranks, each element is a tuple
+            containing the name, the slices and the rank.
         """
 
         ranks = []
+        if not self.fermion_ansatz:
+            return ranks
 
         notations = {
-            "S": [1],
-            "D": [2],
-            "T": [3],
-            "Q": [4],
-            "2": [1, 2],
-            "3": [1, 2, 3],
-            "4": [1, 2, 3, 4],
+            "S": [("t1", "ov", 1)],
+            "D": [("t2", "oovv", 2)],
+            "T": [("t3", "ooovvv", 3)],
+            "t'": [("t3", "OOOVVV", 3)],
         }
-
-        for i, op in enumerate([self.fermion_ansatz, self.boson_ansatz]):
-            # Remove any perturbative corrections
-            while "(" in op:
-                start = op.index("(")
-                end = op.index(")")
-                op = op[:start]
-                if (end + 1) < len(op):
-                    op += op[end + 1 :]
-
-            # Check in order of longest -> shortest string in case
-            # one method name starts with a substring equal to the
-            # name of another method
-            if i == 0:
-                for method_type in sorted(METHOD_TYPES, key=len)[::-1]:
-                    if op.startswith(method_type):
-                        op = op.lstrip(method_type)
-                        break
-
-            # If it's Moller-Plesset perturbation theory, we only
-            # need to initialise second-order amplitudes
-            if method_type == "MP":
-                op = "D"
-
-            # Remove any lower case characters, as these correspond
-            # to active space
-            op = "".join([char for char in op if char.isupper() or char.isnumeric()])
-
-            # Determine the ranks
-            ranks_entry = set()
-            for char in op:
-                for rank in notations[char]:
-                    ranks_entry.add(rank)
-            ranks.append(tuple(sorted(list(ranks_entry))))
-
-        # Get the coupling ranks
-        for op in [self.fermion_coupling_rank, self.boson_coupling_rank]:
-            ranks.append(tuple(range(1, op + 1)))
-
-        return tuple(ranks)
-
-    @property
-    def active_cluster_ranks(self):
-        """Get a list of cluster operator rank numbers for each of
-        the fermionic, bosonic, and coupling ansatzes, for the
-        active space (see space.py).
+        if spin_type == "R":
+            notations["Q"] = [("t4a", "oooovvvv", 4), ("t4b", "oooovvvv", 4)]
+        else:
+            notations["Q"] = [("t4", "oooovvvv", 4)]
+        notations["2"] = notations["S"] + notations["D"]
+        notations["3"] = notations["2"] + notations["T"]
+        notations["4"] = notations["3"] + notations["Q"]
+
+        # Remove any perturbative corrections
+        op = self.fermion_ansatz
+        while "(" in op:
+            start = op.index("(")
+            end = op.index(")")
+            op = op[:start]
+            if (end + 1) < len(op):
+                op += op[end + 1 :]
+
+        # Check in order of longest to shortest string in case one
+        # method name starts with a substring equal to the name of
+        # another method
+        for method_type in sorted(METHOD_TYPES, key=len)[::-1]:
+            if op.startswith(method_type):
+                op = op.lstrip(method_type)
+                break
+
+        # If it's MP we only ever need to initialise second-order
+        # amplitudes
+        if method_type == "MP":
+            op = "D"
+
+        # Determine the ranks
+        for key in sorted(notations.keys(), key=len)[::-1]:
+            if key in op:
+                ranks += notations[key]
+                op = op.replace(key, "")
+
+        # Check there are no duplicates
+        if len(ranks) != len(set(ranks)):
+            raise util.ModelNotImplemented("Duplicate ranks in %s" % self.fermion_ansatz)
+
+        # Sort the ranks by the cluster operator dimension
+        ranks = sorted(ranks, key=lambda x: x[2])
+
+        return ranks
+
+    def bosonic_cluster_ranks(self, spin_type="G"):
+        """Get a list of cluster operator ranks for the bosonic
+        space.
 
         Returns
         -------
-        ranks : tuple of tuple of int
-            Cluster operator ranks for the fermionic, bosonic, and
-            coupling ansatzes, for the active space.
+        ranks : list of tuples
+            List of cluster operator ranks, each element is a tuple
+            containing the name, the slices and the rank.
         """
 
         ranks = []
+        if not self.boson_ansatz:
+            return ranks
 
         notations = {
-            "s": [1],
-            "d": [2],
-            "t": [3],
-            "q": [4],
+            "S": [("s1", "b", 1)],
+            "D": [("s2", "bb", 2)],
+            "T": [("s3", "bbb", 3)],
         }
+        notations["2"] = notations["S"] + notations["D"]
+        notations["3"] = notations["2"] + notations["T"]
+
+        # Remove any perturbative corrections
+        op = self.boson_ansatz
+        while "(" in op:
+            start = op.index("(")
+            end = op.index(")")
+            op = op[:start]
+            if (end + 1) < len(op):
+                op += op[end + 1 :]
+
+        # Determine the ranks
+        for key in sorted(notations.keys(), key=len)[::-1]:
+            if key in op:
+                ranks += notations[key]
+                op = op.replace(key, "")
+
+        # Check there are no duplicates
+        if len(ranks) != len(set(ranks)):
+            raise util.ModelNotImplemented("Duplicate ranks in %s" % self.boson_ansatz)
+
+        # Sort the ranks by the cluster operator dimension
+        ranks = sorted(ranks, key=lambda x: x[2])
+
+        return ranks
+
+    def coupling_cluster_ranks(self, spin_type="G"):
+        """Get a list of cluster operator ranks for the coupling
+        between fermionic and bosonic spaces.
+
+        Returns
+        -------
+        ranks : list of tuple
+            List of cluster operator ranks, each element is a tuple
+            containing the name, slice, fermionic rank and bosonic
+            rank.
+        """
+
+        ranks = []
+
+        for fermion_rank in range(1, self.fermion_coupling_rank + 1):
+            for boson_rank in range(1, self.boson_coupling_rank + 1):
+                name = f"u{fermion_rank}{boson_rank}"
+                key = "b" * boson_rank + "o" * fermion_rank + "v" * fermion_rank
+                ranks.append((name, key, fermion_rank, boson_rank))
 
-        for i, op in enumerate([self.fermion_ansatz, self.boson_ansatz]):
-            # Remove any perturbative corrections
-            while "(" in op:
-                start = op.index("(")
-                end = op.index(")")
-                op = op[:start]
-                if (end + 1) < len(op):
-                    op += op[end + 1 :]
-
-            # Check in order of longest -> shortest string in case
-            # one method name starts with a substring equal to the
-            # name of another method
-            if i == 0:
-                for method_type in sorted(METHOD_TYPES, key=len)[::-1]:
-                    if op.startswith(method_type):
-                        op = op.lstrip(method_type)
-                        break
-
-            # Remove any lower case characters, as these correspond
-            # to active space
-            op = "".join([char for char in op if char.islower()])
-
-            # Determine the ranks
-            ranks_entry = set()
-            for char in op:
-                for rank in notations[char]:
-                    ranks_entry.add(rank)
-            ranks.append(tuple(sorted(list(ranks_entry))))
-
-        # Get the coupling ranks
-        # FIXME how to handle? if it's ever supported
-        ranks.append(tuple())
-
-        return tuple(ranks)
+        return ranks

ebcc/brueckner.py

@@ -130,21 +130,21 @@ def transform_amplitudes(self, u, amplitudes=None):
         ca = u[nocc:, nocc:]
 
         # Transform T amplitudes:
-        for n in self.cc.ansatz.correlated_cluster_ranks[0]:
-            args = [self.cc.amplitudes["t%d" % n], tuple(range(n * 2))]
+        for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
+            args = [self.cc.amplitudes[name], tuple(range(n * 2))]
             for i in range(n):
                 args += [ci, (i, i + n * 2)]
             for i in range(n):
                 args += [ca, (i + n, i + n * 3)]
             args += [tuple(range(n * 2, n * 4))]
-            self.cc.amplitudes["t%d" % n] = util.einsum(*args)
+            self.cc.amplitudes[name] = util.einsum(*args)
 
         # Transform S amplitudes:
-        for n in self.cc.ansatz.correlated_cluster_ranks[1]:
+        for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
             raise util.ModelNotImplemented  # TODO
 
         # Transform U amplitudes:
-        for n in self.cc.ansatz.correlated_cluster_ranks[2]:
+        for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
             raise util.ModelNotImplemented  # TODO
 
         return self.cc.amplitudes
@@ -264,6 +264,10 @@ def kernel(self):
 
         return self.cc.e_corr
 
+    @property
+    def spin_type(self):
+        return self.cc.spin_type
+
 
 @util.inherit_docstrings
 class BruecknerUEBCC(BruecknerREBCC):
@@ -335,22 +339,22 @@ def transform_amplitudes(self, u, amplitudes=None):
         ca = {"a": u[0][nocc[0] :, nocc[0] :], "b": u[1][nocc[1] :, nocc[1] :]}
 
         # Transform T amplitudes:
-        for n in self.cc.ansatz.correlated_cluster_ranks[0]:
+        for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
             for comb in util.generate_spin_combinations(n, unique=True):
-                args = [getattr(self.cc.amplitudes["t%d" % n], comb), tuple(range(n * 2))]
+                args = [getattr(self.cc.amplitudes[name], comb), tuple(range(n * 2))]
             for i in range(n):
                 args += [ci[comb[i]], (i, i + n * 2)]
             for i in range(n):
                 args += [ca[comb[i + n]], (i + n, i + n * 3)]
             args += [tuple(range(n * 2, n * 4))]
-            setattr(self.cc.amplitudes["t%d" % n], comb, util.einsum(*args))
+            setattr(self.cc.amplitudes[name], comb, util.einsum(*args))
 
         # Transform S amplitudes:
-        for n in self.cc.ansatz.correlated_cluster_ranks[1]:
+        for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
             raise util.ModelNotImplemented  # TODO
 
         # Transform U amplitudes:
-        for n in self.cc.ansatz.correlated_cluster_ranks[2]:
+        for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
             raise util.ModelNotImplemented  # TODO
 
         return self.cc.amplitudes