diff --git a/FEATURES.md b/FEATURES.md
index 983fd9e9..b1e80faa 100644
--- a/FEATURES.md
+++ b/FEATURES.md
@@ -18,17 +18,20 @@ The following table summarises the available methods and routines for the ansatz
 | CCSD        |  RUG  |  RUG  |  RUG  |  RUG  |  RUG  |   G   |  RUG  |  RUG  |   -   |
 | CCSDT       |  RUG  |  RUG  |  RUG  |       |       |       |  RUG  |  RUG  |   -   |
 | CCSDTQ      |   g   |   g   |       |       |       |       |       |       |   -   |
-| CCSD(T)     |  RuG  |  RuG  |       |       |       |       |       |       |   -   |
-| CCSDt       |  RG   |  RG   |       |       |       |       |       |       |   -   |
+| CCSD(T)     |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
 | CCSDt'      |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
 | CC2         |  RUG  |  RUG  |  RUG  |  RUG  |  RUG  |   G   |  RUG  |  RUG  |   -   |
 | CC3         |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
+| LCCD        |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
+| LCCSD       |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
 | QCISD       |  RUG  |  RUG  |       |       |       |       |       |       |   -   |
 | DCD         |  RU   |  RU   |       |       |       |       |       |       |   -   |
 | DCSD        |  RU   |  RU   |       |       |       |       |       |       |   -   |
 | DF-CCD      |  RU   |  RU   |  RU   |  RU   |  RU   |       |  RU   |  RU   |   -   |
 | DF-CCSD     |  RU   |  RU   |  RU   |  RU   |  RU   |       |  RU   |  RU   |   -   |
 | DF-CC2      |  RU   |  RU   |  RU   |  RU   |  RU   |       |  RU   |  RU   |   -   |
+| DF-LCCD     |  RU   |  RU   |       |       |       |       |       |       |   -   |
+| DF-LCCSD    |  RU   |  RU   |       |       |       |       |       |       |   -   |
 | DF-QCISD    |  RU   |  RU   |       |       |       |       |       |       |   -   |
 | DF-DCD      |  RU   |  RU   |       |       |       |       |       |       |   -   |
 | DF-DCSD     |  RU   |  RU   |       |       |       |       |       |       |   -   |
diff --git a/ebcc/codegen/GLCCD.py b/ebcc/codegen/GLCCD.py
new file mode 100644
index 00000000..7500a418
--- /dev/null
+++ b/ebcc/codegen/GLCCD.py
@@ -0,0 +1,63 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-17T16:11:07.723784
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(t2, (0, 1, 2, 3), v.oovv, (0, 1, 2, 3), ()) * 0.25
+
+    return e_cc
+
+def update_amps(f=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t2: 
+        v: 
+
+    Returns:
+        t2new: 
+    """
+
+    t2new = einsum(t2, (0, 1, 2, 3), v.oooo, (4, 5, 0, 1), (4, 5, 2, 3)) * 0.5
+    t2new += v.oovv
+    t2new += einsum(t2, (0, 1, 2, 3), v.vvvv, (4, 5, 2, 3), (0, 1, 4, 5)) * 0.5
+    tmp0 = einsum(t2, (0, 1, 2, 3), f.vv, (4, 3), (0, 1, 4, 2))
+    t2new += np.transpose(tmp0, (1, 0, 2, 3))
+    t2new += np.transpose(tmp0, (1, 0, 3, 2)) * -1
+    del tmp0
+    tmp1 = einsum(t2, (0, 1, 2, 3), f.oo, (4, 1), (4, 0, 2, 3))
+    t2new += np.transpose(tmp1, (0, 1, 3, 2)) * -1
+    t2new += np.transpose(tmp1, (1, 0, 3, 2))
+    del tmp1
+    tmp2 = einsum(t2, (0, 1, 2, 3), v.ovov, (4, 3, 1, 5), (0, 4, 2, 5))
+    t2new += tmp2 * -1
+    t2new += np.transpose(tmp2, (0, 1, 3, 2))
+    t2new += np.transpose(tmp2, (1, 0, 2, 3))
+    t2new += np.transpose(tmp2, (1, 0, 3, 2)) * -1
+    del tmp2
+
+    return {"t2new": t2new}
+
diff --git a/ebcc/codegen/GLCCSD.py b/ebcc/codegen/GLCCSD.py
new file mode 100644
index 00000000..15f5ea0b
--- /dev/null
+++ b/ebcc/codegen/GLCCSD.py
@@ -0,0 +1,85 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-17T16:19:25.157269
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(t1, (0, 1), f.ov, (0, 1), ())
+    e_cc += einsum(v.oovv, (0, 1, 2, 3), t2, (0, 1, 2, 3), ()) * 0.25
+
+    return e_cc
+
+def update_amps(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        t1new: 
+        t2new: 
+    """
+
+    t1new = einsum(t1, (0, 1), f.vv, (2, 1), (0, 2))
+    t1new += einsum(t2, (0, 1, 2, 3), v.ovoo, (4, 3, 0, 1), (4, 2)) * -0.5
+    t1new += einsum(v.ovov, (0, 1, 2, 3), t1, (2, 1), (0, 3)) * -1
+    t1new += einsum(t2, (0, 1, 2, 3), f.ov, (1, 3), (0, 2))
+    t1new += f.ov
+    t1new += einsum(f.oo, (0, 1), t1, (1, 2), (0, 2)) * -1
+    t1new += einsum(t2, (0, 1, 2, 3), v.ovvv, (1, 4, 2, 3), (0, 4)) * -0.5
+    t2new = einsum(t2, (0, 1, 2, 3), v.oooo, (4, 5, 0, 1), (4, 5, 2, 3)) * 0.5
+    t2new += v.oovv
+    t2new += einsum(v.vvvv, (0, 1, 2, 3), t2, (4, 5, 2, 3), (4, 5, 0, 1)) * 0.5
+    tmp0 = einsum(t2, (0, 1, 2, 3), v.ovov, (4, 3, 1, 5), (0, 4, 2, 5))
+    t2new += tmp0 * -1
+    t2new += np.transpose(tmp0, (0, 1, 3, 2))
+    t2new += np.transpose(tmp0, (1, 0, 2, 3))
+    t2new += np.transpose(tmp0, (1, 0, 3, 2)) * -1
+    del tmp0
+    tmp1 = einsum(v.ooov, (0, 1, 2, 3), t1, (2, 4), (0, 1, 4, 3))
+    tmp3 = np.copy(np.transpose(tmp1, (1, 0, 2, 3))) * -1
+    del tmp1
+    tmp2 = einsum(t2, (0, 1, 2, 3), f.vv, (4, 3), (0, 1, 4, 2))
+    tmp3 += np.transpose(tmp2, (1, 0, 2, 3)) * -1
+    del tmp2
+    t2new += tmp3 * -1
+    t2new += np.transpose(tmp3, (0, 1, 3, 2))
+    del tmp3
+    tmp4 = einsum(t2, (0, 1, 2, 3), f.oo, (4, 1), (4, 0, 2, 3))
+    tmp6 = np.copy(np.transpose(tmp4, (0, 1, 3, 2))) * -1
+    del tmp4
+    tmp5 = einsum(v.ovvv, (0, 1, 2, 3), t1, (4, 1), (4, 0, 2, 3))
+    tmp6 += np.transpose(tmp5, (0, 1, 3, 2))
+    del tmp5
+    t2new += tmp6
+    t2new += np.transpose(tmp6, (1, 0, 2, 3)) * -1
+    del tmp6
+
+    return {"t1new": t1new, "t2new": t2new}
+
diff --git a/ebcc/codegen/RDFLCCD.py b/ebcc/codegen/RDFLCCD.py
new file mode 100644
index 00000000..8d012cbe
--- /dev/null
+++ b/ebcc/codegen/RDFLCCD.py
@@ -0,0 +1,87 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-18T00:18:12.591121
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    tmp0 = np.copy(np.transpose(t2, (1, 0, 2, 3))) * 2
+    tmp0 += np.transpose(t2, (1, 0, 3, 2)) * -1
+    tmp1 = einsum(v.xov, (0, 1, 2), tmp0, (1, 3, 4, 2), (3, 4, 0)) * 0.5
+    del tmp0
+    e_cc = einsum(tmp1, (0, 1, 2), v.xov, (2, 0, 1), ()) * 2
+    del tmp1
+
+    return e_cc
+
+def update_amps(f=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t2: 
+        v: 
+
+    Returns:
+        t2new: 
+    """
+
+    tmp0 = einsum(v.xoo, (0, 1, 2), v.xoo, (0, 3, 4), (1, 3, 4, 2))
+    t2new = einsum(tmp0, (0, 1, 2, 3), t2, (1, 3, 4, 5), (2, 0, 4, 5))
+    del tmp0
+    t2new += einsum(v.xov, (0, 1, 2), v.xov, (0, 3, 4), (3, 1, 4, 2))
+    tmp1 = einsum(v.xvv, (0, 1, 2), v.xvv, (0, 3, 4), (1, 3, 4, 2))
+    t2new += einsum(t2, (0, 1, 2, 3), tmp1, (4, 2, 5, 3), (0, 1, 5, 4))
+    del tmp1
+    tmp2 = einsum(v.xvv, (0, 1, 2), v.xoo, (0, 3, 4), (3, 4, 1, 2))
+    tmp3 = einsum(tmp2, (0, 1, 2, 3), t2, (4, 1, 5, 3), (4, 0, 5, 2))
+    tmp7 = np.copy(tmp3)
+    del tmp3
+    tmp4 = np.copy(np.transpose(t2, (0, 1, 3, 2))) * 2
+    tmp4 += t2 * -1
+    tmp5 = einsum(v.xov, (0, 1, 2), tmp4, (1, 3, 4, 2), (3, 4, 0)) * 0.5
+    del tmp4
+    tmp6 = einsum(v.xov, (0, 1, 2), tmp5, (3, 4, 0), (1, 3, 2, 4)) * 2
+    del tmp5
+    tmp7 += np.transpose(tmp6, (1, 0, 3, 2)) * -1
+    del tmp6
+    t2new += tmp7 * -1
+    t2new += np.transpose(tmp7, (1, 0, 3, 2)) * -1
+    del tmp7
+    tmp8 = einsum(f.oo, (0, 1), t2, (2, 1, 3, 4), (0, 2, 3, 4))
+    tmp10 = np.copy(tmp8)
+    del tmp8
+    tmp9 = einsum(t2, (0, 1, 2, 3), tmp2, (4, 1, 5, 2), (0, 4, 3, 5))
+    del tmp2
+    tmp10 += tmp9
+    del tmp9
+    t2new += np.transpose(tmp10, (0, 1, 3, 2)) * -1
+    t2new += np.transpose(tmp10, (1, 0, 2, 3)) * -1
+    del tmp10
+    tmp11 = einsum(f.vv, (0, 1), t2, (2, 3, 4, 1), (2, 3, 0, 4))
+    t2new += np.transpose(tmp11, (0, 1, 3, 2))
+    t2new += np.transpose(tmp11, (1, 0, 2, 3))
+    del tmp11
+
+    return {"t2new": t2new}
+
diff --git a/ebcc/codegen/RDFLCCSD.py b/ebcc/codegen/RDFLCCSD.py
new file mode 100644
index 00000000..4edd5cc6
--- /dev/null
+++ b/ebcc/codegen/RDFLCCSD.py
@@ -0,0 +1,117 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-18T00:18:30.019400
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(f.ov, (0, 1), t1, (0, 1), ()) * 2
+    tmp0 = np.copy(np.transpose(t2, (1, 0, 2, 3))) * 2
+    tmp0 += np.transpose(t2, (1, 0, 3, 2)) * -1
+    tmp1 = einsum(tmp0, (0, 1, 2, 3), v.xov, (4, 0, 3), (1, 2, 4)) * 0.5
+    del tmp0
+    e_cc += einsum(v.xov, (0, 1, 2), tmp1, (1, 2, 0), ()) * 2
+    del tmp1
+
+    return e_cc
+
+def update_amps(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        t1new: 
+        t2new: 
+    """
+
+    t1new = einsum(f.vv, (0, 1), t1, (2, 1), (2, 0))
+    t1new += f.ov
+    tmp0 = einsum(v.xov, (0, 1, 2), t1, (1, 2), (0,))
+    t1new += einsum(v.xov, (0, 1, 2), tmp0, (0,), (1, 2)) * 2
+    del tmp0
+    t1new += einsum(t1, (0, 1), f.oo, (2, 0), (2, 1)) * -1
+    tmp1 = einsum(v.xoo, (0, 1, 2), v.xov, (0, 3, 4), (1, 2, 3, 4))
+    tmp2 = np.copy(np.transpose(t2, (0, 1, 3, 2))) * 2
+    tmp2 += t2 * -1
+    t1new += einsum(tmp1, (0, 1, 2, 3), tmp2, (2, 0, 4, 3), (1, 4)) * -1
+    del tmp1, tmp2
+    tmp3 = einsum(t1, (0, 1), v.xoo, (2, 3, 0), (3, 1, 2))
+    tmp6 = np.copy(tmp3)
+    del tmp3
+    tmp4 = np.copy(np.transpose(t2, (0, 1, 3, 2)))
+    tmp4 += t2 * -0.5
+    tmp5 = einsum(tmp4, (0, 1, 2, 3), v.xov, (4, 0, 3), (1, 2, 4)) * 2
+    tmp6 += tmp5 * -1
+    del tmp5
+    t1new += einsum(v.xvv, (0, 1, 2), tmp6, (3, 2, 0), (3, 1)) * -1
+    t1new += einsum(tmp4, (0, 1, 2, 3), f.ov, (0, 3), (1, 2)) * 2
+    del tmp4
+    tmp7 = einsum(v.xvv, (0, 1, 2), v.xvv, (0, 3, 4), (1, 3, 4, 2))
+    t2new = einsum(t2, (0, 1, 2, 3), tmp7, (4, 3, 5, 2), (0, 1, 4, 5))
+    del tmp7
+    t2new += einsum(v.xov, (0, 1, 2), v.xov, (0, 3, 4), (1, 3, 2, 4))
+    tmp8 = einsum(v.xoo, (0, 1, 2), v.xoo, (0, 3, 4), (1, 3, 4, 2))
+    t2new += einsum(tmp8, (0, 1, 2, 3), t2, (1, 3, 4, 5), (2, 0, 4, 5))
+    del tmp8
+    tmp9 = einsum(v.xoo, (0, 1, 2), v.xvv, (0, 3, 4), (1, 2, 3, 4))
+    tmp10 = einsum(tmp9, (0, 1, 2, 3), t2, (4, 1, 5, 3), (4, 0, 5, 2))
+    tmp12 = np.copy(tmp10)
+    del tmp10
+    tmp11 = einsum(v.xov, (0, 1, 2), tmp6, (3, 4, 0), (1, 3, 2, 4))
+    del tmp6
+    tmp12 += np.transpose(tmp11, (1, 0, 3, 2))
+    del tmp11
+    t2new += tmp12 * -1
+    t2new += np.transpose(tmp12, (1, 0, 3, 2)) * -1
+    del tmp12
+    tmp13 = einsum(f.oo, (0, 1), t2, (2, 1, 3, 4), (0, 2, 3, 4))
+    tmp15 = np.copy(tmp13)
+    del tmp13
+    tmp14 = einsum(t2, (0, 1, 2, 3), tmp9, (4, 1, 5, 2), (0, 4, 3, 5))
+    del tmp9
+    tmp15 += tmp14
+    del tmp14
+    t2new += np.transpose(tmp15, (0, 1, 3, 2)) * -1
+    t2new += np.transpose(tmp15, (1, 0, 2, 3)) * -1
+    del tmp15
+    tmp16 = einsum(f.vv, (0, 1), t2, (2, 3, 4, 1), (2, 3, 0, 4))
+    tmp19 = np.copy(tmp16)
+    del tmp16
+    tmp17 = einsum(t1, (0, 1), v.xvv, (2, 3, 1), (0, 3, 2))
+    tmp18 = einsum(tmp17, (0, 1, 2), v.xov, (2, 3, 4), (0, 3, 4, 1))
+    del tmp17
+    tmp19 += tmp18
+    del tmp18
+    t2new += np.transpose(tmp19, (0, 1, 3, 2))
+    t2new += np.transpose(tmp19, (1, 0, 2, 3))
+    del tmp19
+
+    return {"t1new": t1new, "t2new": t2new}
+
diff --git a/ebcc/codegen/RLCCD.py b/ebcc/codegen/RLCCD.py
new file mode 100644
index 00000000..49c3c9b7
--- /dev/null
+++ b/ebcc/codegen/RLCCD.py
@@ -0,0 +1,75 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-17T16:11:29.501322
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(t2, (0, 1, 2, 3), v.ovov, (0, 2, 1, 3), ()) * 2
+    e_cc += einsum(v.ovov, (0, 1, 2, 3), t2, (0, 2, 3, 1), ()) * -1
+
+    return e_cc
+
+def update_amps(f=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t2: 
+        v: 
+
+    Returns:
+        t2new: 
+    """
+
+    t2new = einsum(v.oooo, (0, 1, 2, 3), t2, (3, 1, 4, 5), (0, 2, 5, 4))
+    t2new += np.transpose(v.ovov, (0, 2, 1, 3))
+    t2new += einsum(t2, (0, 1, 2, 3), v.vvvv, (4, 2, 5, 3), (0, 1, 4, 5))
+    tmp0 = einsum(t2, (0, 1, 2, 3), f.vv, (4, 3), (0, 1, 4, 2))
+    t2new += np.transpose(tmp0, (0, 1, 3, 2))
+    t2new += np.transpose(tmp0, (1, 0, 2, 3))
+    del tmp0
+    tmp1 = einsum(t2, (0, 1, 2, 3), f.oo, (4, 1), (4, 0, 2, 3))
+    tmp3 = np.copy(tmp1)
+    del tmp1
+    tmp2 = einsum(t2, (0, 1, 2, 3), v.oovv, (4, 1, 5, 2), (0, 4, 3, 5))
+    tmp3 += tmp2
+    del tmp2
+    t2new += np.transpose(tmp3, (0, 1, 3, 2)) * -1
+    t2new += np.transpose(tmp3, (1, 0, 2, 3)) * -1
+    del tmp3
+    tmp4 = einsum(v.oovv, (0, 1, 2, 3), t2, (4, 1, 5, 3), (4, 0, 5, 2))
+    tmp7 = np.copy(tmp4)
+    del tmp4
+    tmp5 = np.copy(np.transpose(t2, (0, 1, 3, 2))) * 2
+    tmp5 += t2 * -1
+    tmp6 = einsum(tmp5, (0, 1, 2, 3), v.ovov, (4, 5, 0, 3), (1, 4, 2, 5))
+    del tmp5
+    tmp7 += tmp6 * -1
+    del tmp6
+    t2new += tmp7 * -1
+    t2new += np.transpose(tmp7, (1, 0, 3, 2)) * -1
+    del tmp7
+
+    return {"t2new": t2new}
+
diff --git a/ebcc/codegen/RLCCSD.py b/ebcc/codegen/RLCCSD.py
new file mode 100644
index 00000000..8f2edd2d
--- /dev/null
+++ b/ebcc/codegen/RLCCSD.py
@@ -0,0 +1,107 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-17T16:22:36.741829
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(v.ovov, (0, 1, 2, 3), t2, (0, 2, 1, 3), ()) * 2
+    e_cc += einsum(t2, (0, 1, 2, 3), v.ovov, (0, 3, 1, 2), ()) * -1
+    e_cc += einsum(t1, (0, 1), f.ov, (0, 1), ()) * 2
+
+    return e_cc
+
+def update_amps(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        t1new: 
+        t2new: 
+    """
+
+    t1new = np.copy(f.ov)
+    t1new += einsum(t1, (0, 1), f.oo, (2, 0), (2, 1)) * -1
+    t1new += einsum(f.vv, (0, 1), t1, (2, 1), (2, 0))
+    tmp0 = np.copy(v.ovvv) * -0.5
+    tmp0 += np.transpose(v.ovvv, (0, 2, 1, 3))
+    t1new += einsum(t2, (0, 1, 2, 3), tmp0, (1, 2, 3, 4), (0, 4)) * 2
+    del tmp0
+    tmp1 = np.copy(v.ooov) * -0.5
+    tmp1 += np.transpose(v.ovoo, (0, 2, 3, 1))
+    t1new += einsum(tmp1, (0, 1, 2, 3), t2, (0, 2, 3, 4), (1, 4)) * -2
+    del tmp1
+    tmp2 = np.copy(np.transpose(v.ovov, (0, 2, 1, 3)))
+    tmp2 += v.oovv * -0.5
+    t1new += einsum(t1, (0, 1), tmp2, (0, 2, 1, 3), (2, 3)) * 2
+    del tmp2
+    tmp3 = np.copy(np.transpose(t2, (0, 1, 3, 2))) * -0.5
+    tmp3 += t2
+    t1new += einsum(f.ov, (0, 1), tmp3, (0, 2, 1, 3), (2, 3)) * 2
+    del tmp3
+    t2new = einsum(v.vvvv, (0, 1, 2, 3), t2, (4, 5, 1, 3), (4, 5, 0, 2))
+    t2new += np.transpose(v.ovov, (0, 2, 1, 3))
+    t2new += einsum(t2, (0, 1, 2, 3), v.oooo, (4, 1, 5, 0), (4, 5, 3, 2))
+    tmp4 = einsum(v.ooov, (0, 1, 2, 3), t1, (1, 4), (0, 2, 4, 3))
+    tmp8 = np.copy(tmp4)
+    del tmp4
+    tmp5 = einsum(v.oovv, (0, 1, 2, 3), t2, (4, 1, 5, 3), (4, 0, 5, 2))
+    tmp8 += tmp5
+    del tmp5
+    tmp6 = np.copy(np.transpose(t2, (0, 1, 3, 2))) * -1
+    tmp6 += t2 * 2
+    tmp7 = einsum(tmp6, (0, 1, 2, 3), v.ovov, (4, 5, 0, 2), (4, 1, 5, 3))
+    del tmp6
+    tmp8 += np.transpose(tmp7, (1, 0, 3, 2)) * -1
+    del tmp7
+    t2new += tmp8 * -1
+    t2new += np.transpose(tmp8, (1, 0, 3, 2)) * -1
+    del tmp8
+    tmp9 = einsum(t2, (0, 1, 2, 3), f.oo, (4, 1), (4, 0, 2, 3))
+    tmp11 = np.copy(tmp9)
+    del tmp9
+    tmp10 = einsum(v.oovv, (0, 1, 2, 3), t2, (4, 1, 3, 5), (4, 0, 5, 2))
+    tmp11 += tmp10
+    del tmp10
+    t2new += np.transpose(tmp11, (0, 1, 3, 2)) * -1
+    t2new += np.transpose(tmp11, (1, 0, 2, 3)) * -1
+    del tmp11
+    tmp12 = einsum(f.vv, (0, 1), t2, (2, 3, 4, 1), (2, 3, 0, 4))
+    tmp14 = np.copy(tmp12)
+    del tmp12
+    tmp13 = einsum(t1, (0, 1), v.ovvv, (2, 3, 4, 1), (0, 2, 3, 4))
+    tmp14 += tmp13
+    del tmp13
+    t2new += np.transpose(tmp14, (0, 1, 3, 2))
+    t2new += np.transpose(tmp14, (1, 0, 2, 3))
+    del tmp14
+
+    return {"t1new": t1new, "t2new": t2new}
+
diff --git a/ebcc/codegen/UDFLCCD.py b/ebcc/codegen/UDFLCCD.py
new file mode 100644
index 00000000..6c19fbad
--- /dev/null
+++ b/ebcc/codegen/UDFLCCD.py
@@ -0,0 +1,161 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-18T00:18:06.986031
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    tmp0 = einsum(t2.bbbb, (0, 1, 2, 3), v.bb.xov, (4, 1, 3), (0, 2, 4))
+    e_cc = einsum(v.bb.xov, (0, 1, 2), tmp0, (1, 2, 0), ())
+    del tmp0
+    tmp1 = einsum(t2.aaaa, (0, 1, 2, 3), v.aa.xov, (4, 1, 3), (0, 2, 4))
+    tmp1 += einsum(t2.abab, (0, 1, 2, 3), v.bb.xov, (4, 1, 3), (0, 2, 4))
+    e_cc += einsum(v.aa.xov, (0, 1, 2), tmp1, (1, 2, 0), ())
+    del tmp1
+
+    return e_cc
+
+def update_amps(f=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t2: 
+        v: 
+
+    Returns:
+        t2new: 
+    """
+
+    tmp0 = einsum(f.aa.oo, (0, 1), t2.aaaa, (2, 1, 3, 4), (0, 2, 3, 4))
+    t2new = Namespace()
+    t2new.aaaa = np.copy(np.transpose(tmp0, (1, 0, 3, 2))) * 2
+    t2new.aaaa += np.transpose(tmp0, (0, 1, 3, 2)) * -2
+    del tmp0
+    tmp1 = einsum(v.aa.xoo, (0, 1, 2), v.aa.xoo, (0, 3, 4), (1, 3, 4, 2))
+    t2new.aaaa += einsum(tmp1, (0, 1, 2, 3), t2.aaaa, (1, 3, 4, 5), (2, 0, 4, 5)) * 2
+    del tmp1
+    tmp2 = einsum(t2.aaaa, (0, 1, 2, 3), f.aa.vv, (4, 3), (0, 1, 4, 2))
+    t2new.aaaa += np.transpose(tmp2, (1, 0, 2, 3)) * 2
+    t2new.aaaa += np.transpose(tmp2, (1, 0, 3, 2)) * -2
+    del tmp2
+    tmp3 = einsum(v.aa.xoo, (0, 1, 2), v.aa.xvv, (0, 3, 4), (1, 2, 3, 4))
+    tmp4 = einsum(tmp3, (0, 1, 2, 3), t2.aaaa, (4, 1, 5, 3), (4, 0, 5, 2))
+    t2new.aaaa += np.transpose(tmp4, (1, 0, 3, 2)) * -2
+    tmp7 = np.copy(np.transpose(v.aa.xov, (1, 2, 0))) * 0.5
+    tmp5 = einsum(t2.aaaa, (0, 1, 2, 3), v.aa.xov, (4, 1, 3), (0, 2, 4))
+    tmp7 += tmp5
+    tmp6 = einsum(t2.abab, (0, 1, 2, 3), v.bb.xov, (4, 1, 3), (0, 2, 4))
+    tmp7 += tmp6 * 0.5
+    t2new.aaaa += einsum(tmp7, (0, 1, 2), v.aa.xov, (2, 3, 4), (3, 0, 4, 1)) * 2
+    del tmp7
+    t2new.aaaa += np.transpose(tmp4, (1, 0, 2, 3)) * 2
+    tmp8 = np.copy(np.transpose(v.aa.xov, (1, 2, 0)))
+    tmp8 += tmp5 * 2
+    tmp8 += tmp6
+    t2new.aaaa += einsum(v.aa.xov, (0, 1, 2), tmp8, (3, 4, 0), (1, 3, 4, 2)) * -1
+    del tmp8
+    t2new.aaaa += np.transpose(tmp4, (0, 1, 3, 2)) * 2
+    tmp9 = np.copy(tmp5)
+    tmp9 += tmp6 * 0.5
+    t2new.aaaa += einsum(v.aa.xov, (0, 1, 2), tmp9, (3, 4, 0), (3, 1, 2, 4)) * -2
+    del tmp9
+    t2new.aaaa += tmp4 * -2
+    del tmp4
+    tmp10 = np.copy(tmp5) * 2
+    del tmp5
+    tmp10 += tmp6
+    del tmp6
+    t2new.aaaa += einsum(tmp10, (0, 1, 2), v.aa.xov, (2, 3, 4), (0, 3, 1, 4))
+    tmp11 = einsum(v.aa.xvv, (0, 1, 2), v.aa.xvv, (0, 3, 4), (3, 1, 2, 4))
+    t2new.aaaa += einsum(t2.aaaa, (0, 1, 2, 3), tmp11, (4, 2, 5, 3), (0, 1, 4, 5)) * -2
+    del tmp11
+    tmp20 = einsum(t2.bbbb, (0, 1, 2, 3), f.bb.oo, (4, 1), (4, 0, 2, 3))
+    t2new.bbbb = np.copy(np.transpose(tmp20, (1, 0, 3, 2))) * 2
+    tmp21 = einsum(v.bb.xoo, (0, 1, 2), v.bb.xoo, (0, 3, 4), (1, 3, 4, 2))
+    t2new.bbbb += einsum(tmp21, (0, 1, 2, 3), t2.bbbb, (1, 3, 4, 5), (2, 0, 4, 5)) * 2
+    del tmp21
+    t2new.bbbb += np.transpose(tmp20, (0, 1, 3, 2)) * -2
+    del tmp20
+    tmp22 = einsum(f.bb.vv, (0, 1), t2.bbbb, (2, 3, 4, 1), (2, 3, 0, 4))
+    t2new.bbbb += np.transpose(tmp22, (1, 0, 2, 3)) * 2
+    t2new.bbbb += np.transpose(tmp22, (1, 0, 3, 2)) * -2
+    del tmp22
+    tmp17 = einsum(v.bb.xoo, (0, 1, 2), v.bb.xvv, (0, 3, 4), (1, 2, 3, 4))
+    tmp23 = einsum(t2.bbbb, (0, 1, 2, 3), tmp17, (4, 1, 5, 3), (0, 4, 2, 5))
+    t2new.bbbb += np.transpose(tmp23, (1, 0, 3, 2)) * -2
+    tmp24 = np.copy(np.transpose(v.bb.xov, (1, 2, 0))) * 0.5
+    tmp13 = einsum(v.aa.xov, (0, 1, 2), t2.abab, (1, 3, 2, 4), (3, 4, 0))
+    tmp24 += tmp13 * 0.5
+    tmp14 = einsum(t2.bbbb, (0, 1, 2, 3), v.bb.xov, (4, 1, 3), (0, 2, 4))
+    tmp24 += tmp14
+    t2new.bbbb += einsum(v.bb.xov, (0, 1, 2), tmp24, (3, 4, 0), (1, 3, 2, 4)) * 2
+    del tmp24
+    t2new.bbbb += np.transpose(tmp23, (1, 0, 2, 3)) * 2
+    tmp15 = np.copy(np.transpose(v.bb.xov, (1, 2, 0)))
+    tmp15 += tmp13
+    tmp15 += tmp14 * 2
+    t2new.bbbb += einsum(tmp15, (0, 1, 2), v.bb.xov, (2, 3, 4), (3, 0, 1, 4)) * -1
+    t2new.bbbb += np.transpose(tmp23, (0, 1, 3, 2)) * 2
+    tmp25 = np.copy(tmp13) * 0.5
+    tmp25 += tmp14
+    t2new.bbbb += einsum(v.bb.xov, (0, 1, 2), tmp25, (3, 4, 0), (3, 1, 2, 4)) * -2
+    del tmp25
+    t2new.bbbb += tmp23 * -2
+    del tmp23
+    tmp26 = np.copy(tmp13)
+    del tmp13
+    tmp26 += tmp14 * 2
+    del tmp14
+    t2new.bbbb += einsum(v.bb.xov, (0, 1, 2), tmp26, (3, 4, 0), (3, 1, 4, 2))
+    del tmp26
+    tmp27 = einsum(v.bb.xvv, (0, 1, 2), v.bb.xvv, (0, 3, 4), (3, 1, 2, 4))
+    t2new.bbbb += einsum(t2.bbbb, (0, 1, 2, 3), tmp27, (4, 2, 5, 3), (0, 1, 4, 5)) * -2
+    del tmp27
+    t2new.abab = einsum(f.bb.oo, (0, 1), t2.abab, (2, 1, 3, 4), (2, 0, 3, 4)) * -1
+    t2new.abab += einsum(f.aa.oo, (0, 1), t2.abab, (1, 2, 3, 4), (0, 2, 3, 4)) * -1
+    tmp12 = einsum(v.bb.xoo, (0, 1, 2), v.aa.xoo, (0, 3, 4), (3, 4, 1, 2))
+    t2new.abab += einsum(tmp12, (0, 1, 2, 3), t2.abab, (1, 3, 4, 5), (0, 2, 4, 5))
+    del tmp12
+    t2new.abab += einsum(f.aa.vv, (0, 1), t2.abab, (2, 3, 1, 4), (2, 3, 0, 4))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), f.bb.vv, (4, 3), (0, 1, 2, 4))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp3, (4, 0, 5, 2), (4, 1, 5, 3)) * -1
+    del tmp3
+    t2new.abab += einsum(v.aa.xov, (0, 1, 2), tmp15, (3, 4, 0), (1, 3, 2, 4))
+    del tmp15
+    tmp16 = einsum(v.bb.xoo, (0, 1, 2), v.aa.xvv, (0, 3, 4), (1, 2, 3, 4))
+    t2new.abab += einsum(tmp16, (0, 1, 2, 3), t2.abab, (4, 1, 3, 5), (4, 0, 2, 5)) * -1
+    del tmp16
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp17, (4, 1, 5, 3), (0, 4, 2, 5)) * -1
+    del tmp17
+    t2new.abab += einsum(v.bb.xov, (0, 1, 2), tmp10, (3, 4, 0), (3, 1, 4, 2))
+    del tmp10
+    tmp18 = einsum(v.aa.xoo, (0, 1, 2), v.bb.xvv, (0, 3, 4), (1, 2, 3, 4))
+    t2new.abab += einsum(tmp18, (0, 1, 2, 3), t2.abab, (1, 4, 5, 3), (0, 4, 5, 2)) * -1
+    del tmp18
+    tmp19 = einsum(v.bb.xvv, (0, 1, 2), v.aa.xvv, (0, 3, 4), (3, 4, 1, 2))
+    t2new.abab += einsum(tmp19, (0, 1, 2, 3), t2.abab, (4, 5, 1, 3), (4, 5, 0, 2))
+    del tmp19
+
+    return {"t2new": t2new}
+
diff --git a/ebcc/codegen/UDFLCCSD.py b/ebcc/codegen/UDFLCCSD.py
new file mode 100644
index 00000000..cd0c6426
--- /dev/null
+++ b/ebcc/codegen/UDFLCCSD.py
@@ -0,0 +1,239 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-18T00:18:21.497352
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(t1.aa, (0, 1), f.aa.ov, (0, 1), ())
+    e_cc += einsum(f.bb.ov, (0, 1), t1.bb, (0, 1), ())
+    tmp0 = einsum(v.aa.xov, (0, 1, 2), t2.aaaa, (3, 1, 4, 2), (3, 4, 0))
+    e_cc += einsum(v.aa.xov, (0, 1, 2), tmp0, (1, 2, 0), ())
+    del tmp0
+    tmp1 = einsum(t2.abab, (0, 1, 2, 3), v.aa.xov, (4, 0, 2), (1, 3, 4))
+    tmp1 += einsum(t2.bbbb, (0, 1, 2, 3), v.bb.xov, (4, 1, 3), (0, 2, 4))
+    e_cc += einsum(tmp1, (0, 1, 2), v.bb.xov, (2, 0, 1), ())
+    del tmp1
+
+    return e_cc
+
+def update_amps(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        t1new: 
+        t2new: 
+    """
+
+    tmp0 = einsum(v.aa.xov, (0, 1, 2), v.aa.xoo, (0, 3, 4), (3, 4, 1, 2))
+    t1new = Namespace()
+    t1new.aa = einsum(t2.aaaa, (0, 1, 2, 3), tmp0, (4, 0, 1, 3), (4, 2)) * -2
+    tmp1 = einsum(v.bb.xov, (0, 1, 2), v.aa.xoo, (0, 3, 4), (3, 4, 1, 2))
+    t1new.aa += einsum(t2.abab, (0, 1, 2, 3), tmp1, (4, 0, 1, 3), (4, 2)) * -1
+    t1new.aa += einsum(t1.aa, (0, 1), f.aa.oo, (2, 0), (2, 1)) * -1
+    t1new.aa += einsum(f.aa.ov, (0, 1), t2.aaaa, (2, 0, 3, 1), (2, 3)) * 2
+    t1new.aa += f.aa.ov
+    tmp2 = einsum(v.aa.xvv, (0, 1, 2), t1.aa, (3, 2), (3, 1, 0))
+    t1new.aa += einsum(tmp2, (0, 1, 2), v.aa.xoo, (2, 3, 0), (3, 1)) * -1
+    t1new.aa += einsum(t2.abab, (0, 1, 2, 3), f.bb.ov, (1, 3), (0, 2))
+    tmp3 = einsum(v.aa.xov, (0, 1, 2), t1.aa, (1, 2), (0,))
+    tmp5 = np.copy(tmp3)
+    del tmp3
+    tmp4 = einsum(t1.bb, (0, 1), v.bb.xov, (2, 0, 1), (2,))
+    tmp5 += tmp4
+    del tmp4
+    t1new.aa += einsum(v.aa.xov, (0, 1, 2), tmp5, (0,), (1, 2))
+    t1new.aa += einsum(f.aa.vv, (0, 1), t1.aa, (2, 1), (2, 0))
+    tmp6 = einsum(v.aa.xov, (0, 1, 2), t2.aaaa, (3, 1, 4, 2), (3, 4, 0))
+    tmp8 = np.copy(tmp6) * 2
+    tmp7 = einsum(t2.abab, (0, 1, 2, 3), v.bb.xov, (4, 1, 3), (0, 2, 4))
+    tmp8 += tmp7
+    t1new.aa += einsum(tmp8, (0, 1, 2), v.aa.xvv, (2, 3, 1), (0, 3))
+    del tmp8
+    t1new.bb = einsum(t2.bbbb, (0, 1, 2, 3), f.bb.ov, (1, 3), (0, 2)) * 2
+    tmp9 = einsum(t1.bb, (0, 1), v.bb.xvv, (2, 3, 1), (0, 3, 2))
+    t1new.bb += einsum(v.bb.xoo, (0, 1, 2), tmp9, (2, 3, 0), (1, 3)) * -1
+    tmp10 = einsum(v.bb.xoo, (0, 1, 2), v.aa.xov, (0, 3, 4), (3, 1, 2, 4))
+    t1new.bb += einsum(tmp10, (0, 1, 2, 3), t2.abab, (0, 2, 3, 4), (1, 4)) * -1
+    del tmp10
+    t1new.bb += f.bb.ov
+    tmp11 = einsum(v.bb.xoo, (0, 1, 2), v.bb.xov, (0, 3, 4), (1, 2, 3, 4))
+    t1new.bb += einsum(t2.bbbb, (0, 1, 2, 3), tmp11, (4, 0, 1, 3), (4, 2)) * -2
+    t1new.bb += einsum(t2.abab, (0, 1, 2, 3), f.aa.ov, (0, 2), (1, 3))
+    t1new.bb += einsum(t1.bb, (0, 1), f.bb.oo, (2, 0), (2, 1)) * -1
+    t1new.bb += einsum(v.bb.xov, (0, 1, 2), tmp5, (0,), (1, 2))
+    del tmp5
+    t1new.bb += einsum(t1.bb, (0, 1), f.bb.vv, (2, 1), (0, 2))
+    tmp12 = einsum(t2.abab, (0, 1, 2, 3), v.aa.xov, (4, 0, 2), (1, 3, 4))
+    tmp14 = np.copy(tmp12) * 0.5
+    tmp13 = einsum(v.bb.xov, (0, 1, 2), t2.bbbb, (3, 1, 4, 2), (3, 4, 0))
+    tmp14 += tmp13
+    t1new.bb += einsum(v.bb.xvv, (0, 1, 2), tmp14, (3, 2, 0), (3, 1)) * 2
+    tmp15 = einsum(t2.aaaa, (0, 1, 2, 3), f.aa.oo, (4, 1), (4, 0, 2, 3))
+    t2new = Namespace()
+    t2new.aaaa = np.copy(np.transpose(tmp15, (1, 0, 3, 2))) * 2
+    tmp16 = einsum(tmp0, (0, 1, 2, 3), t1.aa, (1, 4), (0, 2, 4, 3))
+    del tmp0
+    t2new.aaaa += np.transpose(tmp16, (0, 1, 3, 2))
+    t2new.aaaa += np.transpose(tmp15, (0, 1, 3, 2)) * -2
+    del tmp15
+    tmp17 = einsum(v.aa.xoo, (0, 1, 2), v.aa.xoo, (0, 3, 4), (3, 1, 2, 4))
+    t2new.aaaa += einsum(t2.aaaa, (0, 1, 2, 3), tmp17, (4, 0, 5, 1), (5, 4, 2, 3)) * 2
+    del tmp17
+    tmp18 = einsum(t2.aaaa, (0, 1, 2, 3), f.aa.vv, (4, 3), (0, 1, 4, 2))
+    t2new.aaaa += np.transpose(tmp18, (1, 0, 2, 3)) * 2
+    t2new.aaaa += np.transpose(tmp18, (1, 0, 3, 2)) * -2
+    del tmp18
+    tmp19 = einsum(v.aa.xoo, (0, 1, 2), v.aa.xvv, (0, 3, 4), (1, 2, 3, 4))
+    tmp20 = einsum(tmp19, (0, 1, 2, 3), t2.aaaa, (4, 1, 5, 3), (4, 0, 5, 2))
+    t2new.aaaa += np.transpose(tmp20, (1, 0, 3, 2)) * -2
+    tmp22 = np.copy(np.transpose(v.aa.xov, (1, 2, 0)))
+    tmp21 = einsum(v.aa.xoo, (0, 1, 2), t1.aa, (2, 3), (1, 3, 0))
+    tmp22 += tmp21 * -1
+    del tmp21
+    tmp22 += tmp6 * 2
+    tmp22 += tmp7
+    tmp23 = einsum(tmp22, (0, 1, 2), v.aa.xov, (2, 3, 4), (0, 3, 1, 4))
+    del tmp22
+    t2new.aaaa += np.transpose(tmp23, (1, 0, 3, 2))
+    t2new.aaaa += np.transpose(tmp20, (1, 0, 2, 3)) * 2
+    t2new.aaaa += np.transpose(tmp23, (1, 0, 2, 3)) * -1
+    del tmp23
+    t2new.aaaa += tmp16 * -1
+    del tmp16
+    tmp24 = einsum(v.aa.xov, (0, 1, 2), tmp2, (3, 4, 0), (3, 1, 2, 4))
+    t2new.aaaa += np.transpose(tmp24, (0, 1, 3, 2))
+    t2new.aaaa += tmp24 * -1
+    t2new.aaaa += np.transpose(tmp24, (1, 0, 3, 2)) * -1
+    t2new.aaaa += np.transpose(tmp24, (1, 0, 2, 3))
+    del tmp24
+    t2new.aaaa += np.transpose(tmp20, (0, 1, 3, 2)) * 2
+    tmp25 = np.copy(tmp6)
+    del tmp6
+    tmp25 += tmp7 * 0.5
+    del tmp7
+    tmp26 = einsum(v.aa.xov, (0, 1, 2), tmp25, (3, 4, 0), (3, 1, 4, 2)) * 2
+    t2new.aaaa += np.transpose(tmp26, (0, 1, 3, 2)) * -1
+    t2new.aaaa += tmp20 * -2
+    del tmp20
+    t2new.aaaa += tmp26
+    del tmp26
+    tmp27 = einsum(v.aa.xvv, (0, 1, 2), v.aa.xvv, (0, 3, 4), (3, 1, 2, 4))
+    t2new.aaaa += einsum(tmp27, (0, 1, 2, 3), t2.aaaa, (4, 5, 1, 3), (4, 5, 0, 2)) * -2
+    del tmp27
+    tmp35 = einsum(t2.bbbb, (0, 1, 2, 3), f.bb.oo, (4, 1), (4, 0, 2, 3))
+    t2new.bbbb = np.copy(np.transpose(tmp35, (1, 0, 3, 2))) * 2
+    tmp36 = einsum(v.bb.xoo, (0, 1, 2), v.bb.xoo, (0, 3, 4), (3, 1, 2, 4))
+    t2new.bbbb += einsum(tmp36, (0, 1, 2, 3), t2.bbbb, (1, 3, 4, 5), (2, 0, 4, 5)) * 2
+    del tmp36
+    t2new.bbbb += np.transpose(tmp35, (0, 1, 3, 2)) * -2
+    del tmp35
+    tmp37 = einsum(t1.bb, (0, 1), tmp11, (2, 0, 3, 4), (2, 3, 1, 4))
+    del tmp11
+    t2new.bbbb += np.transpose(tmp37, (0, 1, 3, 2))
+    tmp38 = einsum(t2.bbbb, (0, 1, 2, 3), f.bb.vv, (4, 3), (0, 1, 4, 2))
+    t2new.bbbb += np.transpose(tmp38, (1, 0, 2, 3)) * 2
+    t2new.bbbb += np.transpose(tmp38, (1, 0, 3, 2)) * -2
+    del tmp38
+    tmp32 = einsum(v.bb.xoo, (0, 1, 2), v.bb.xvv, (0, 3, 4), (1, 2, 3, 4))
+    tmp39 = einsum(tmp32, (0, 1, 2, 3), t2.bbbb, (4, 1, 5, 3), (4, 0, 5, 2))
+    t2new.bbbb += np.transpose(tmp39, (1, 0, 3, 2)) * -2
+    tmp40 = np.copy(np.transpose(v.bb.xov, (1, 2, 0)))
+    tmp29 = einsum(v.bb.xoo, (0, 1, 2), t1.bb, (2, 3), (1, 3, 0))
+    tmp40 += tmp29 * -1
+    tmp40 += tmp12
+    tmp40 += tmp13 * 2
+    tmp41 = einsum(tmp40, (0, 1, 2), v.bb.xov, (2, 3, 4), (0, 3, 1, 4))
+    del tmp40
+    t2new.bbbb += np.transpose(tmp41, (1, 0, 3, 2))
+    t2new.bbbb += np.transpose(tmp39, (1, 0, 2, 3)) * 2
+    t2new.bbbb += np.transpose(tmp41, (1, 0, 2, 3)) * -1
+    del tmp41
+    t2new.bbbb += tmp37 * -1
+    del tmp37
+    tmp42 = einsum(v.bb.xov, (0, 1, 2), tmp9, (3, 4, 0), (3, 1, 2, 4))
+    t2new.bbbb += np.transpose(tmp42, (0, 1, 3, 2))
+    t2new.bbbb += tmp42 * -1
+    t2new.bbbb += np.transpose(tmp42, (1, 0, 3, 2)) * -1
+    t2new.bbbb += np.transpose(tmp42, (1, 0, 2, 3))
+    del tmp42
+    t2new.bbbb += np.transpose(tmp39, (0, 1, 3, 2)) * 2
+    tmp43 = np.copy(tmp12)
+    tmp43 += tmp13 * 2
+    t2new.bbbb += einsum(v.bb.xov, (0, 1, 2), tmp43, (3, 4, 0), (3, 1, 2, 4)) * -1
+    del tmp43
+    t2new.bbbb += tmp39 * -2
+    del tmp39
+    t2new.bbbb += einsum(tmp14, (0, 1, 2), v.bb.xov, (2, 3, 4), (0, 3, 1, 4)) * 2
+    del tmp14
+    tmp44 = einsum(v.bb.xvv, (0, 1, 2), v.bb.xvv, (0, 3, 4), (3, 1, 2, 4))
+    t2new.bbbb += einsum(tmp44, (0, 1, 2, 3), t2.bbbb, (4, 5, 1, 3), (4, 5, 0, 2)) * -2
+    del tmp44
+    t2new.abab = einsum(f.bb.oo, (0, 1), t2.abab, (2, 1, 3, 4), (2, 0, 3, 4)) * -1
+    tmp28 = einsum(v.bb.xoo, (0, 1, 2), v.aa.xoo, (0, 3, 4), (3, 4, 1, 2))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp28, (4, 0, 5, 1), (4, 5, 2, 3))
+    del tmp28
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), f.aa.oo, (4, 0), (4, 1, 2, 3)) * -1
+    t2new.abab += einsum(f.aa.vv, (0, 1), t2.abab, (2, 3, 1, 4), (2, 3, 0, 4))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), f.bb.vv, (4, 3), (0, 1, 2, 4))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp19, (4, 0, 5, 2), (4, 1, 5, 3)) * -1
+    del tmp19
+    tmp30 = np.copy(np.transpose(v.bb.xov, (1, 2, 0))) * 0.5
+    tmp30 += tmp29 * -0.5
+    del tmp29
+    tmp30 += tmp12 * 0.5
+    del tmp12
+    tmp30 += tmp13
+    del tmp13
+    t2new.abab += einsum(v.aa.xov, (0, 1, 2), tmp30, (3, 4, 0), (1, 3, 2, 4)) * 2
+    del tmp30
+    t2new.abab += einsum(tmp1, (0, 1, 2, 3), t1.aa, (1, 4), (0, 2, 4, 3)) * -1
+    del tmp1
+    t2new.abab += einsum(tmp2, (0, 1, 2), v.bb.xov, (2, 3, 4), (0, 3, 1, 4))
+    del tmp2
+    t2new.abab += einsum(v.aa.xov, (0, 1, 2), tmp9, (3, 4, 0), (1, 3, 2, 4))
+    del tmp9
+    tmp31 = einsum(v.bb.xoo, (0, 1, 2), v.aa.xvv, (0, 3, 4), (1, 2, 3, 4))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp31, (4, 1, 5, 2), (0, 4, 5, 3)) * -1
+    del tmp31
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp32, (4, 1, 5, 3), (0, 4, 2, 5)) * -1
+    del tmp32
+    t2new.abab += einsum(tmp25, (0, 1, 2), v.bb.xov, (2, 3, 4), (0, 3, 1, 4)) * 2
+    del tmp25
+    tmp33 = einsum(v.bb.xvv, (0, 1, 2), v.aa.xoo, (0, 3, 4), (3, 4, 1, 2))
+    t2new.abab += einsum(tmp33, (0, 1, 2, 3), t2.abab, (1, 4, 5, 3), (0, 4, 5, 2)) * -1
+    del tmp33
+    tmp34 = einsum(v.bb.xvv, (0, 1, 2), v.aa.xvv, (0, 3, 4), (3, 4, 1, 2))
+    t2new.abab += einsum(tmp34, (0, 1, 2, 3), t2.abab, (4, 5, 1, 3), (4, 5, 0, 2))
+    del tmp34
+
+    return {"t1new": t1new, "t2new": t2new}
+
diff --git a/ebcc/codegen/ULCCD.py b/ebcc/codegen/ULCCD.py
new file mode 100644
index 00000000..3ec2881e
--- /dev/null
+++ b/ebcc/codegen/ULCCD.py
@@ -0,0 +1,116 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-17T16:12:07.457835
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(t2.abab, (0, 1, 2, 3), v.aabb.ovov, (0, 2, 1, 3), ())
+    e_cc += einsum(v.aaaa.ovov, (0, 1, 2, 3), t2.aaaa, (0, 2, 3, 1), ()) * -1
+    e_cc += einsum(t2.bbbb, (0, 1, 2, 3), v.bbbb.ovov, (0, 3, 1, 2), ()) * -1
+
+    return e_cc
+
+def update_amps(f=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t2: 
+        v: 
+
+    Returns:
+        t2new: 
+    """
+
+    tmp0 = einsum(t2.aaaa, (0, 1, 2, 3), f.aa.oo, (4, 1), (4, 0, 2, 3))
+    t2new = Namespace()
+    t2new.aaaa = np.copy(np.transpose(tmp0, (1, 0, 3, 2))) * 2
+    t2new.aaaa += np.transpose(tmp0, (0, 1, 3, 2)) * -2
+    del tmp0
+    t2new.aaaa += einsum(t2.aaaa, (0, 1, 2, 3), v.aaaa.oooo, (4, 0, 5, 1), (4, 5, 2, 3)) * 2
+    tmp1 = einsum(t2.aaaa, (0, 1, 2, 3), f.aa.vv, (4, 3), (0, 1, 4, 2))
+    t2new.aaaa += np.transpose(tmp1, (1, 0, 2, 3)) * 2
+    t2new.aaaa += np.transpose(tmp1, (1, 0, 3, 2)) * -2
+    del tmp1
+    t2new.aaaa += np.transpose(v.aaaa.ovov, (0, 2, 1, 3))
+    tmp2 = einsum(t2.abab, (0, 1, 2, 3), v.aabb.ovov, (4, 5, 1, 3), (0, 4, 2, 5))
+    t2new.aaaa += np.transpose(tmp2, (1, 0, 3, 2))
+    tmp3 = np.copy(np.transpose(v.aaaa.ovov, (0, 2, 1, 3)))
+    tmp3 += v.aaaa.oovv * -1
+    tmp4 = einsum(tmp3, (0, 1, 2, 3), t2.aaaa, (4, 0, 5, 2), (4, 1, 5, 3)) * 2
+    t2new.aaaa += np.transpose(tmp4, (1, 0, 3, 2))
+    t2new.aaaa += np.transpose(v.aaaa.ovov, (0, 2, 3, 1)) * -1
+    t2new.aaaa += np.transpose(tmp2, (1, 0, 2, 3)) * -1
+    t2new.aaaa += np.transpose(tmp4, (1, 0, 2, 3)) * -1
+    t2new.aaaa += np.transpose(tmp2, (0, 1, 3, 2)) * -1
+    t2new.aaaa += np.transpose(tmp4, (0, 1, 3, 2)) * -1
+    t2new.aaaa += tmp2
+    del tmp2
+    t2new.aaaa += tmp4
+    del tmp4
+    t2new.aaaa += einsum(v.aaaa.vvvv, (0, 1, 2, 3), t2.aaaa, (4, 5, 1, 3), (4, 5, 0, 2)) * 2
+    t2new.abab = einsum(f.bb.oo, (0, 1), t2.abab, (2, 1, 3, 4), (2, 0, 3, 4)) * -1
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), v.aabb.oooo, (4, 0, 5, 1), (4, 5, 2, 3))
+    t2new.abab += einsum(f.aa.oo, (0, 1), t2.abab, (1, 2, 3, 4), (0, 2, 3, 4)) * -1
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), f.aa.vv, (4, 2), (0, 1, 4, 3))
+    t2new.abab += einsum(f.bb.vv, (0, 1), t2.abab, (2, 3, 4, 1), (2, 3, 4, 0))
+    t2new.abab += einsum(v.aabb.ovov, (0, 1, 2, 3), t2.bbbb, (4, 2, 5, 3), (0, 4, 1, 5)) * 2
+    t2new.abab += np.transpose(v.aabb.ovov, (0, 2, 1, 3))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp3, (0, 4, 2, 5), (4, 1, 5, 3))
+    del tmp3
+    t2new.abab += einsum(v.aabb.vvoo, (0, 1, 2, 3), t2.abab, (4, 3, 1, 5), (4, 2, 0, 5)) * -1
+    t2new.abab += einsum(t2.aaaa, (0, 1, 2, 3), v.aabb.ovov, (1, 3, 4, 5), (0, 4, 2, 5)) * 2
+    tmp5 = np.copy(np.transpose(v.bbbb.ovov, (0, 2, 1, 3)))
+    tmp5 += v.bbbb.oovv * -1
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp5, (1, 4, 3, 5), (0, 4, 2, 5))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), v.aabb.oovv, (4, 0, 5, 3), (4, 1, 2, 5)) * -1
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), v.aabb.vvvv, (4, 2, 5, 3), (0, 1, 4, 5))
+    tmp6 = einsum(f.bb.oo, (0, 1), t2.bbbb, (2, 1, 3, 4), (0, 2, 3, 4))
+    t2new.bbbb = np.copy(np.transpose(tmp6, (1, 0, 3, 2))) * 2
+    t2new.bbbb += np.transpose(tmp6, (0, 1, 3, 2)) * -2
+    del tmp6
+    t2new.bbbb += einsum(t2.bbbb, (0, 1, 2, 3), v.bbbb.oooo, (4, 0, 5, 1), (4, 5, 2, 3)) * 2
+    tmp7 = einsum(t2.bbbb, (0, 1, 2, 3), f.bb.vv, (4, 3), (0, 1, 4, 2))
+    t2new.bbbb += np.transpose(tmp7, (1, 0, 2, 3)) * 2
+    t2new.bbbb += np.transpose(tmp7, (1, 0, 3, 2)) * -2
+    del tmp7
+    tmp8 = einsum(t2.abab, (0, 1, 2, 3), v.aabb.ovov, (0, 2, 4, 5), (1, 4, 3, 5))
+    t2new.bbbb += np.transpose(tmp8, (1, 0, 3, 2))
+    t2new.bbbb += np.transpose(v.bbbb.ovov, (0, 2, 1, 3))
+    tmp9 = einsum(t2.bbbb, (0, 1, 2, 3), tmp5, (1, 4, 3, 5), (0, 4, 2, 5)) * 2
+    del tmp5
+    t2new.bbbb += np.transpose(tmp9, (1, 0, 3, 2))
+    t2new.bbbb += np.transpose(v.bbbb.ovov, (0, 2, 3, 1)) * -1
+    t2new.bbbb += np.transpose(tmp8, (1, 0, 2, 3)) * -1
+    t2new.bbbb += np.transpose(tmp9, (1, 0, 2, 3)) * -1
+    t2new.bbbb += np.transpose(tmp8, (0, 1, 3, 2)) * -1
+    t2new.bbbb += np.transpose(tmp9, (0, 1, 3, 2)) * -1
+    t2new.bbbb += tmp8
+    del tmp8
+    t2new.bbbb += tmp9
+    del tmp9
+    t2new.bbbb += einsum(t2.bbbb, (0, 1, 2, 3), v.bbbb.vvvv, (4, 2, 5, 3), (0, 1, 4, 5)) * 2
+
+    return {"t2new": t2new}
+
diff --git a/ebcc/codegen/ULCCSD.py b/ebcc/codegen/ULCCSD.py
new file mode 100644
index 00000000..44c90b1c
--- /dev/null
+++ b/ebcc/codegen/ULCCSD.py
@@ -0,0 +1,173 @@
+"""Code generated by `albert` version 0.0.0.
+
+ * date: 2024-12-17T16:28:14.750474
+ * python version: 3.10.12 (main, Nov  6 2024, 20:22:13) [GCC 11.4.0]
+ * albert version: 0.0.0
+ * caller: /home/ollie/git/albert/albert/code/einsum.py
+ * node: ollie-desktop
+ * system: Linux
+ * processor: x86_64
+ * release: 6.8.0-49-generic
+"""
+
+from ebcc import numpy as np
+from ebcc.util import pack_2e, einsum, dirsum, Namespace
+
+
+def energy(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        e_cc: 
+    """
+
+    e_cc = einsum(t1.aa, (0, 1), f.aa.ov, (0, 1), ())
+    e_cc += einsum(v.bbbb.ovov, (0, 1, 2, 3), t2.bbbb, (0, 2, 3, 1), ()) * -1
+    e_cc += einsum(v.aaaa.ovov, (0, 1, 2, 3), t2.aaaa, (0, 2, 3, 1), ()) * -1
+    e_cc += einsum(v.aabb.ovov, (0, 1, 2, 3), t2.abab, (0, 2, 1, 3), ())
+    e_cc += einsum(f.bb.ov, (0, 1), t1.bb, (0, 1), ())
+
+    return e_cc
+
+def update_amps(f=None, t1=None, t2=None, v=None, **kwargs):
+    """Code generated by `albert` 0.0.0.
+
+    Args:
+        f: 
+        t1: 
+        t2: 
+        v: 
+
+    Returns:
+        t1new: 
+        t2new: 
+    """
+
+    t1new = Namespace()
+    t1new.aa = einsum(f.aa.oo, (0, 1), t1.aa, (1, 2), (0, 2)) * -1
+    t1new.aa += f.aa.ov
+    t1new.aa += einsum(t2.aaaa, (0, 1, 2, 3), f.aa.ov, (1, 3), (0, 2)) * 2
+    t1new.aa += einsum(v.aaaa.ooov, (0, 1, 2, 3), t2.aaaa, (2, 1, 4, 3), (0, 4)) * 2
+    t1new.aa += einsum(v.aabb.ooov, (0, 1, 2, 3), t2.abab, (1, 2, 4, 3), (0, 4)) * -1
+    t1new.aa += einsum(f.bb.ov, (0, 1), t2.abab, (2, 0, 3, 1), (2, 3))
+    t1new.aa += einsum(t1.bb, (0, 1), v.aabb.ovov, (2, 3, 0, 1), (2, 3))
+    tmp0 = np.copy(np.transpose(v.aaaa.ovov, (0, 2, 1, 3)))
+    tmp0 += v.aaaa.oovv * -1
+    t1new.aa += einsum(tmp0, (0, 1, 2, 3), t1.aa, (0, 2), (1, 3))
+    t1new.aa += einsum(v.aaaa.ovvv, (0, 1, 2, 3), t2.aaaa, (4, 0, 1, 3), (4, 2)) * -2
+    t1new.aa += einsum(f.aa.vv, (0, 1), t1.aa, (2, 1), (2, 0))
+    t1new.aa += einsum(v.aabb.vvov, (0, 1, 2, 3), t2.abab, (4, 2, 1, 3), (4, 0))
+    t1new.bb = einsum(f.aa.ov, (0, 1), t2.abab, (0, 2, 1, 3), (2, 3))
+    t1new.bb += einsum(t1.aa, (0, 1), v.aabb.ovov, (0, 1, 2, 3), (2, 3))
+    t1new.bb += einsum(t1.bb, (0, 1), f.bb.oo, (2, 0), (2, 1)) * -1
+    t1new.bb += f.bb.ov
+    t1new.bb += einsum(v.aabb.ovoo, (0, 1, 2, 3), t2.abab, (0, 3, 1, 4), (2, 4)) * -1
+    t1new.bb += einsum(t2.bbbb, (0, 1, 2, 3), f.bb.ov, (1, 3), (0, 2)) * 2
+    t1new.bb += einsum(t2.bbbb, (0, 1, 2, 3), v.bbbb.ooov, (4, 0, 1, 3), (4, 2)) * -2
+    tmp1 = np.copy(np.transpose(v.bbbb.ovov, (0, 2, 1, 3)))
+    tmp1 += v.bbbb.oovv * -1
+    t1new.bb += einsum(t1.bb, (0, 1), tmp1, (0, 2, 1, 3), (2, 3))
+    t1new.bb += einsum(t2.bbbb, (0, 1, 2, 3), v.bbbb.ovvv, (1, 3, 4, 2), (0, 4)) * 2
+    t1new.bb += einsum(t2.abab, (0, 1, 2, 3), v.aabb.ovvv, (0, 2, 4, 3), (1, 4))
+    t1new.bb += einsum(t1.bb, (0, 1), f.bb.vv, (2, 1), (0, 2))
+    tmp2 = einsum(f.aa.oo, (0, 1), t2.aaaa, (2, 1, 3, 4), (0, 2, 3, 4))
+    t2new = Namespace()
+    t2new.aaaa = np.copy(np.transpose(tmp2, (1, 0, 3, 2))) * 2
+    tmp3 = einsum(v.aaaa.ooov, (0, 1, 2, 3), t1.aa, (1, 4), (0, 2, 4, 3))
+    t2new.aaaa += np.transpose(tmp3, (0, 1, 3, 2))
+    t2new.aaaa += np.transpose(tmp2, (0, 1, 3, 2)) * -2
+    del tmp2
+    t2new.aaaa += einsum(t2.aaaa, (0, 1, 2, 3), v.aaaa.oooo, (4, 1, 5, 0), (4, 5, 2, 3)) * -2
+    tmp4 = einsum(f.aa.vv, (0, 1), t2.aaaa, (2, 3, 4, 1), (2, 3, 0, 4))
+    t2new.aaaa += np.transpose(tmp4, (1, 0, 2, 3)) * 2
+    t2new.aaaa += np.transpose(tmp4, (1, 0, 3, 2)) * -2
+    del tmp4
+    t2new.aaaa += np.transpose(tmp3, (1, 0, 3, 2)) * -1
+    t2new.aaaa += np.transpose(v.aaaa.ovov, (0, 2, 1, 3))
+    tmp5 = einsum(v.aabb.ovov, (0, 1, 2, 3), t2.abab, (4, 2, 5, 3), (4, 0, 5, 1))
+    t2new.aaaa += np.transpose(tmp5, (1, 0, 3, 2))
+    tmp6 = einsum(t2.aaaa, (0, 1, 2, 3), tmp0, (1, 4, 3, 5), (0, 4, 2, 5)) * 2
+    t2new.aaaa += np.transpose(tmp6, (1, 0, 3, 2))
+    t2new.aaaa += np.transpose(tmp3, (1, 0, 2, 3))
+    t2new.aaaa += np.transpose(v.aaaa.ovov, (0, 2, 3, 1)) * -1
+    t2new.aaaa += np.transpose(tmp5, (1, 0, 2, 3)) * -1
+    t2new.aaaa += np.transpose(tmp6, (1, 0, 2, 3)) * -1
+    t2new.aaaa += tmp3 * -1
+    del tmp3
+    tmp7 = einsum(v.aaaa.ovvv, (0, 1, 2, 3), t1.aa, (4, 3), (4, 0, 1, 2))
+    t2new.aaaa += np.transpose(tmp7, (0, 1, 3, 2))
+    t2new.aaaa += tmp7 * -1
+    t2new.aaaa += np.transpose(tmp7, (1, 0, 3, 2)) * -1
+    t2new.aaaa += np.transpose(tmp7, (1, 0, 2, 3))
+    del tmp7
+    t2new.aaaa += np.transpose(tmp5, (0, 1, 3, 2)) * -1
+    t2new.aaaa += np.transpose(tmp6, (0, 1, 3, 2)) * -1
+    t2new.aaaa += tmp5
+    del tmp5
+    t2new.aaaa += tmp6
+    del tmp6
+    t2new.aaaa += einsum(t2.aaaa, (0, 1, 2, 3), v.aaaa.vvvv, (4, 3, 5, 2), (0, 1, 4, 5)) * -2
+    t2new.abab = einsum(f.bb.oo, (0, 1), t2.abab, (2, 1, 3, 4), (2, 0, 3, 4)) * -1
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), f.aa.oo, (4, 0), (4, 1, 2, 3)) * -1
+    t2new.abab += einsum(v.aabb.oooo, (0, 1, 2, 3), t2.abab, (1, 3, 4, 5), (0, 2, 4, 5))
+    t2new.abab += einsum(f.aa.vv, (0, 1), t2.abab, (2, 3, 1, 4), (2, 3, 0, 4))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), f.bb.vv, (4, 3), (0, 1, 2, 4))
+    t2new.abab += einsum(v.aabb.ovoo, (0, 1, 2, 3), t1.bb, (3, 4), (0, 2, 1, 4)) * -1
+    t2new.abab += einsum(t2.bbbb, (0, 1, 2, 3), v.aabb.ovov, (4, 5, 1, 3), (4, 0, 5, 2)) * 2
+    t2new.abab += np.transpose(v.aabb.ovov, (0, 2, 1, 3))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), tmp0, (0, 4, 2, 5), (4, 1, 5, 3))
+    del tmp0
+    t2new.abab += einsum(v.aabb.ooov, (0, 1, 2, 3), t1.aa, (1, 4), (0, 2, 4, 3)) * -1
+    t2new.abab += einsum(v.aabb.vvov, (0, 1, 2, 3), t1.aa, (4, 1), (4, 2, 0, 3))
+    t2new.abab += einsum(t1.bb, (0, 1), v.aabb.ovvv, (2, 3, 4, 1), (2, 0, 3, 4))
+    t2new.abab += einsum(v.aabb.vvoo, (0, 1, 2, 3), t2.abab, (4, 3, 1, 5), (4, 2, 0, 5)) * -1
+    t2new.abab += einsum(t2.aaaa, (0, 1, 2, 3), v.aabb.ovov, (1, 3, 4, 5), (0, 4, 2, 5)) * 2
+    t2new.abab += einsum(tmp1, (0, 1, 2, 3), t2.abab, (4, 0, 5, 2), (4, 1, 5, 3))
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), v.aabb.oovv, (4, 0, 5, 3), (4, 1, 2, 5)) * -1
+    t2new.abab += einsum(t2.abab, (0, 1, 2, 3), v.aabb.vvvv, (4, 2, 5, 3), (0, 1, 4, 5))
+    tmp8 = einsum(t2.bbbb, (0, 1, 2, 3), f.bb.oo, (4, 1), (4, 0, 2, 3))
+    t2new.bbbb = np.copy(np.transpose(tmp8, (1, 0, 3, 2))) * 2
+    tmp9 = einsum(t1.bb, (0, 1), v.bbbb.ooov, (2, 0, 3, 4), (2, 3, 1, 4))
+    t2new.bbbb += np.transpose(tmp9, (0, 1, 3, 2))
+    t2new.bbbb += np.transpose(tmp8, (0, 1, 3, 2)) * -2
+    del tmp8
+    t2new.bbbb += einsum(v.bbbb.oooo, (0, 1, 2, 3), t2.bbbb, (3, 1, 4, 5), (0, 2, 4, 5)) * -2
+    tmp10 = einsum(f.bb.vv, (0, 1), t2.bbbb, (2, 3, 4, 1), (2, 3, 0, 4))
+    t2new.bbbb += np.transpose(tmp10, (1, 0, 2, 3)) * 2
+    t2new.bbbb += np.transpose(tmp10, (1, 0, 3, 2)) * -2
+    del tmp10
+    t2new.bbbb += np.transpose(tmp9, (1, 0, 3, 2)) * -1
+    tmp11 = einsum(v.aabb.ovov, (0, 1, 2, 3), t2.abab, (0, 4, 1, 5), (4, 2, 5, 3))
+    t2new.bbbb += np.transpose(tmp11, (1, 0, 3, 2))
+    t2new.bbbb += np.transpose(v.bbbb.ovov, (0, 2, 1, 3))
+    tmp12 = einsum(tmp1, (0, 1, 2, 3), t2.bbbb, (4, 0, 5, 2), (4, 1, 5, 3)) * 2
+    del tmp1
+    t2new.bbbb += np.transpose(tmp12, (1, 0, 3, 2))
+    t2new.bbbb += np.transpose(tmp9, (1, 0, 2, 3))
+    t2new.bbbb += np.transpose(tmp11, (1, 0, 2, 3)) * -1
+    t2new.bbbb += np.transpose(v.bbbb.ovov, (0, 2, 3, 1)) * -1
+    t2new.bbbb += np.transpose(tmp12, (1, 0, 2, 3)) * -1
+    t2new.bbbb += tmp9 * -1
+    del tmp9
+    tmp13 = einsum(v.bbbb.ovvv, (0, 1, 2, 3), t1.bb, (4, 3), (4, 0, 1, 2))
+    t2new.bbbb += np.transpose(tmp13, (0, 1, 3, 2))
+    t2new.bbbb += tmp13 * -1
+    t2new.bbbb += np.transpose(tmp13, (1, 0, 3, 2)) * -1
+    t2new.bbbb += np.transpose(tmp13, (1, 0, 2, 3))
+    del tmp13
+    t2new.bbbb += np.transpose(tmp11, (0, 1, 3, 2)) * -1
+    t2new.bbbb += np.transpose(tmp12, (0, 1, 3, 2)) * -1
+    t2new.bbbb += tmp11
+    del tmp11
+    t2new.bbbb += tmp12
+    del tmp12
+    t2new.bbbb += einsum(t2.bbbb, (0, 1, 2, 3), v.bbbb.vvvv, (4, 3, 5, 2), (0, 1, 4, 5)) * -2
+
+    return {"t1new": t1new, "t2new": t2new}
+
diff --git a/ebcc/core/ansatz.py b/ebcc/core/ansatz.py
index da8c043f..5da4591b 100644
--- a/ebcc/core/ansatz.py
+++ b/ebcc/core/ansatz.py
@@ -21,6 +21,8 @@
     "CCSD(T)": ("CCSD(T)", "", 0, 0),
     "CC2": ("CC2", "", 0, 0),
     "CC3": ("CC3", "", 0, 0),
+    "LCCD": ("LCCD", "", 0, 0),
+    "LCCSD": ("LCCSD", "", 0, 0),
     "QCISD": ("QCISD", "", 0, 0),
     "DCD": ("DCD", "", 0, 0),
     "DCSD": ("DCSD", "", 0, 0),
diff --git a/tests/test_GLCCD.py b/tests/test_GLCCD.py
new file mode 100644
index 00000000..5f6b9666
--- /dev/null
+++ b/tests/test_GLCCD.py
@@ -0,0 +1,65 @@
+"""Tests for the ULCCD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, GEBCC, NullLogger, Space, util
+
+
+@pytest.mark.regression
+class ULCCD_Tests(unittest.TestCase):
+    """Test ULCCD against RLCCD.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "N 0 0 0; N 0 0 2.2"
+        mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        rlccd = REBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        rlccd.options.e_tol = 1e-10
+        rlccd.kernel()
+
+        rlccd = GEBCC.from_rebcc(rlccd)
+
+        glccd = GEBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        glccd.options.e_tol = 1e-10
+        glccd.kernel()
+
+        cls.mf, cls.rlccd, cls.glccd= mf, rlccd, glccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.rlccd, cls.glccd
+
+    def test_converged(self):
+        self.assertTrue(self.rlccd.converged)
+        self.assertTrue(self.glccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.rlccd.e_tot, self.glccd.e_tot, places=8)
+
+
+if __name__ == "__main__":
+    print("Tests for GLCCD")
+    unittest.main()
diff --git a/tests/test_GLCCSD.py b/tests/test_GLCCSD.py
new file mode 100644
index 00000000..3b79db36
--- /dev/null
+++ b/tests/test_GLCCSD.py
@@ -0,0 +1,65 @@
+"""Tests for the ULCCSD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, GEBCC, NullLogger, Space, util
+
+
+@pytest.mark.regression
+class ULCCSD_Tests(unittest.TestCase):
+    """Test ULCCSD against RLCCSD.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "N 0 0 0; N 0 0 2.2"
+        mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        rlccd = REBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        rlccd.options.e_tol = 1e-10
+        rlccd.kernel()
+
+        rlccd = GEBCC.from_rebcc(rlccd)
+
+        glccd = GEBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        glccd.options.e_tol = 1e-10
+        glccd.kernel()
+
+        cls.mf, cls.rlccd, cls.glccd= mf, rlccd, glccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.rlccd, cls.glccd
+
+    def test_converged(self):
+        self.assertTrue(self.rlccd.converged)
+        self.assertTrue(self.glccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.rlccd.e_tot, self.glccd.e_tot, places=8)
+
+
+if __name__ == "__main__":
+    print("Tests for GLCCSD")
+    unittest.main()
diff --git a/tests/test_RDFLCCD.py b/tests/test_RDFLCCD.py
new file mode 100644
index 00000000..ee04e1eb
--- /dev/null
+++ b/tests/test_RDFLCCD.py
@@ -0,0 +1,57 @@
+"""Tests for the RLCCD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, NullLogger, Space
+
+
+@pytest.mark.reference
+class RLCCD_Tests(unittest.TestCase):
+    """Test RLCCD against reference values.
+
+    Reference values from PyBEST 2.0.0.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "H 0 0 0; F 0 0 1"
+        #mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf = mf.density_fit(auxbasis="aug-cc-pvqz-ri")
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        lccd = REBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        lccd.options.e_tol = 1e-10
+        lccd.kernel()
+
+        cls.mf, cls.lccd = mf, lccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.lccd
+
+    def test_converged(self):
+        self.assertTrue(self.lccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.lccd.e_corr, -0.21336662492171, 5)
+
+
+if __name__ == "__main__":
+    print("Tests for RLCCD")
+    unittest.main()
diff --git a/tests/test_RDFLCCSD.py b/tests/test_RDFLCCSD.py
new file mode 100644
index 00000000..e674bf52
--- /dev/null
+++ b/tests/test_RDFLCCSD.py
@@ -0,0 +1,57 @@
+"""Tests for the RLCCSD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, NullLogger, Space
+
+
+@pytest.mark.reference
+class RLCCSD_Tests(unittest.TestCase):
+    """Test RLCCSD against reference values.
+
+    Reference values from PyBEST 2.0.0.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "H 0 0 0; F 0 0 1"
+        #mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf = mf.density_fit(auxbasis="aug-cc-pvqz-ri")
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        lccd = REBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        lccd.options.e_tol = 1e-10
+        lccd.kernel()
+
+        cls.mf, cls.lccd = mf, lccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.lccd
+
+    def test_converged(self):
+        self.assertTrue(self.lccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.lccd.e_corr, -0.21451573243690, 5)
+
+
+if __name__ == "__main__":
+    print("Tests for RLCCSD")
+    unittest.main()
diff --git a/tests/test_RLCCD.py b/tests/test_RLCCD.py
new file mode 100644
index 00000000..154451a7
--- /dev/null
+++ b/tests/test_RLCCD.py
@@ -0,0 +1,57 @@
+"""Tests for the RLCCD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, NullLogger, Space
+
+
+@pytest.mark.reference
+class RLCCD_Tests(unittest.TestCase):
+    """Test RLCCD against reference values.
+
+    Reference values from PyBEST 2.0.0.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "H 0 0 0; F 0 0 1"
+        #mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        lccd = REBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        lccd.options.e_tol = 1e-10
+        lccd.kernel()
+
+        cls.mf, cls.lccd = mf, lccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.lccd
+
+    def test_converged(self):
+        self.assertTrue(self.lccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.lccd.e_corr, -0.21336662492171, 6)
+        self.assertAlmostEqual(self.lccd.e_tot, -100.2232366580237, 6)
+
+
+if __name__ == "__main__":
+    print("Tests for RLCCD")
+    unittest.main()
diff --git a/tests/test_RLCCSD.py b/tests/test_RLCCSD.py
new file mode 100644
index 00000000..0aa29fcd
--- /dev/null
+++ b/tests/test_RLCCSD.py
@@ -0,0 +1,57 @@
+"""Tests for the RLCCSD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, NullLogger, Space
+
+
+@pytest.mark.reference
+class RLCCSD_Tests(unittest.TestCase):
+    """Test RLCCSD against reference values.
+
+    Reference values from PyBEST 2.0.0.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "H 0 0 0; F 0 0 1"
+        #mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        lccd = REBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        lccd.options.e_tol = 1e-10
+        lccd.kernel()
+
+        cls.mf, cls.lccd = mf, lccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.lccd
+
+    def test_converged(self):
+        self.assertTrue(self.lccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.lccd.e_corr, -0.21451573243690, 6)
+        self.assertAlmostEqual(self.lccd.e_tot, -100.2243857655389, 6)
+
+
+if __name__ == "__main__":
+    print("Tests for RLCCSD")
+    unittest.main()
diff --git a/tests/test_UDFLCCD.py b/tests/test_UDFLCCD.py
new file mode 100644
index 00000000..5b662b8e
--- /dev/null
+++ b/tests/test_UDFLCCD.py
@@ -0,0 +1,73 @@
+"""Tests for the ULCCD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, UEBCC, NullLogger, Space
+
+
+@pytest.mark.regression
+class ULCCD_Tests(unittest.TestCase):
+    """Test ULCCD against RLCCD.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "N 0 0 0; N 0 0 2.2"
+        mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf = mf.density_fit()
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        rlccd = REBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        rlccd.options.e_tol = 1e-10
+        rlccd.kernel()
+
+        rlccd = UEBCC.from_rebcc(rlccd)
+
+        ulccd = UEBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        ulccd.options.e_tol = 1e-10
+        ulccd.kernel()
+
+        cls.mf, cls.rlccd, cls.ulccd= mf, rlccd, ulccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.rlccd, cls.ulccd
+
+    def test_converged(self):
+        self.assertTrue(self.rlccd.converged)
+        self.assertTrue(self.ulccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.rlccd.e_tot, self.ulccd.e_tot, places=8)
+
+    def test_t2_amplitudes(self):
+        a = self.rlccd.t2
+        b = self.ulccd.t2
+        np.testing.assert_allclose(a.aaaa, b.aaaa, atol=1e-7)
+        np.testing.assert_allclose(a.bbbb, b.bbbb, atol=1e-7)
+        np.testing.assert_allclose(a.abab, b.abab, atol=1e-7)
+
+
+if __name__ == "__main__":
+    print("Tests for ULCCD")
+    unittest.main()
diff --git a/tests/test_UDFLCCSD.py b/tests/test_UDFLCCSD.py
new file mode 100644
index 00000000..1726523d
--- /dev/null
+++ b/tests/test_UDFLCCSD.py
@@ -0,0 +1,73 @@
+"""Tests for the ULCCSD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, UEBCC, NullLogger, Space
+
+
+@pytest.mark.regression
+class ULCCSD_Tests(unittest.TestCase):
+    """Test ULCCSD against RLCCSD.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "N 0 0 0; N 0 0 2.2"
+        mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf = mf.density_fit()
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        rlccd = REBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        rlccd.options.e_tol = 1e-10
+        rlccd.kernel()
+
+        rlccd = UEBCC.from_rebcc(rlccd)
+
+        ulccd = UEBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        ulccd.options.e_tol = 1e-10
+        ulccd.kernel()
+
+        cls.mf, cls.rlccd, cls.ulccd= mf, rlccd, ulccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.rlccd, cls.ulccd
+
+    def test_converged(self):
+        self.assertTrue(self.rlccd.converged)
+        self.assertTrue(self.ulccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.rlccd.e_tot, self.ulccd.e_tot, places=8)
+
+    def test_t2_amplitudes(self):
+        a = self.rlccd.t2
+        b = self.ulccd.t2
+        np.testing.assert_allclose(a.aaaa, b.aaaa, atol=1e-7)
+        np.testing.assert_allclose(a.bbbb, b.bbbb, atol=1e-7)
+        np.testing.assert_allclose(a.abab, b.abab, atol=1e-7)
+
+
+if __name__ == "__main__":
+    print("Tests for ULCCSD")
+    unittest.main()
diff --git a/tests/test_ULCCD.py b/tests/test_ULCCD.py
new file mode 100644
index 00000000..76345e52
--- /dev/null
+++ b/tests/test_ULCCD.py
@@ -0,0 +1,72 @@
+"""Tests for the ULCCD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, UEBCC, NullLogger, Space
+
+
+@pytest.mark.regression
+class ULCCD_Tests(unittest.TestCase):
+    """Test ULCCD against RLCCD.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "N 0 0 0; N 0 0 2.2"
+        mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        rlccd = REBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        rlccd.options.e_tol = 1e-10
+        rlccd.kernel()
+
+        rlccd = UEBCC.from_rebcc(rlccd)
+
+        ulccd = UEBCC(
+                mf,
+                ansatz="LCCD",
+                log=NullLogger(),
+        )
+        ulccd.options.e_tol = 1e-10
+        ulccd.kernel()
+
+        cls.mf, cls.rlccd, cls.ulccd= mf, rlccd, ulccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.rlccd, cls.ulccd
+
+    def test_converged(self):
+        self.assertTrue(self.rlccd.converged)
+        self.assertTrue(self.ulccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.rlccd.e_tot, self.ulccd.e_tot, places=8)
+
+    def test_t2_amplitudes(self):
+        a = self.rlccd.t2
+        b = self.ulccd.t2
+        np.testing.assert_allclose(a.aaaa, b.aaaa, atol=1e-7)
+        np.testing.assert_allclose(a.bbbb, b.bbbb, atol=1e-7)
+        np.testing.assert_allclose(a.abab, b.abab, atol=1e-7)
+
+
+if __name__ == "__main__":
+    print("Tests for ULCCD")
+    unittest.main()
diff --git a/tests/test_ULCCSD.py b/tests/test_ULCCSD.py
new file mode 100644
index 00000000..bc14ace5
--- /dev/null
+++ b/tests/test_ULCCSD.py
@@ -0,0 +1,72 @@
+"""Tests for the ULCCSD model.
+"""
+
+import unittest
+
+import numpy as np
+import pytest
+from pyscf import gto, scf
+
+from ebcc import REBCC, UEBCC, NullLogger, Space
+
+
+@pytest.mark.regression
+class ULCCSD_Tests(unittest.TestCase):
+    """Test ULCCSD against RLCCSD.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        mol = gto.Mole()
+        mol.atom = "N 0 0 0; N 0 0 2.2"
+        mol.unit = "B"
+        mol.basis = "cc-pvdz"
+        mol.verbose = 0
+        mol.build()
+
+        mf = scf.RHF(mol)
+        mf.conv_tol = 1e-12
+        mf.kernel()
+
+        rlccd = REBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        rlccd.options.e_tol = 1e-10
+        rlccd.kernel()
+
+        rlccd = UEBCC.from_rebcc(rlccd)
+
+        ulccd = UEBCC(
+                mf,
+                ansatz="LCCSD",
+                log=NullLogger(),
+        )
+        ulccd.options.e_tol = 1e-10
+        ulccd.kernel()
+
+        cls.mf, cls.rlccd, cls.ulccd= mf, rlccd, ulccd
+
+    @classmethod
+    def tearDownClass(cls):
+        del cls.mf, cls.rlccd, cls.ulccd
+
+    def test_converged(self):
+        self.assertTrue(self.rlccd.converged)
+        self.assertTrue(self.ulccd.converged)
+
+    def test_energy(self):
+        self.assertAlmostEqual(self.rlccd.e_tot, self.ulccd.e_tot, places=8)
+
+    def test_t2_amplitudes(self):
+        a = self.rlccd.t2
+        b = self.ulccd.t2
+        np.testing.assert_allclose(a.aaaa, b.aaaa, atol=1e-7)
+        np.testing.assert_allclose(a.bbbb, b.bbbb, atol=1e-7)
+        np.testing.assert_allclose(a.abab, b.abab, atol=1e-7)
+
+
+if __name__ == "__main__":
+    print("Tests for ULCCSD")
+    unittest.main()