Repository: quantumlib/OpenFermion
Pull request #824 (2023-06-20 21:10:28)

[PBC Resources Estimates 4/4] Add costing functions
Final part of #821. Builds on previous 3 PRs to add costing estimates for the various factorizations. Also adds documentation and a notebook.

fdmalone: @mpharrigan PTAL
fdmalone: Thanks for slogging through this @mpharrigan!

Patch:
diff --git a/src/openfermion/resource_estimates/pbc/README.md b/src/openfermion/resource_estimates/pbc/README.md
new file mode 100644
index 0000000..8d0a657
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/README.md
@@ -0,0 +1,86 @@
+# Resource Estimation for Periodic Systems 
+
+Module `openfermion.resource_estimates.pbc` facilitates fault-tolerant (FT) resource estimates for second-quantized symmetry-adapted Hamiltonians with periodic boundary conditions (PBC).
+
+The module provides symmetry adapted sparse, single, double and tensor hypercontraction representations of the Hamiltonians. 
+
+For the methods listed above, there are sub-routines which:
+* factorize the two-electron integrals if appropriate
+* compute the associated lambda values, `compute_lambda()`
+* estimate the number of logical qubits and Toffoli gates required to simulate with this factorization, `compute_cost()`
+
+### Details
+
+Given a pyscf scf calculation of a periodic system with k-points:
+
+```python
+from pyscf.pbc import gto, scf
+
+cell = gto.Cell()
+cell.atom = '''
+C 0.000000000000   0.000000000000   0.000000000000
+C 1.685068664391   1.685068664391   1.685068664391
+'''
+cell.basis = 'gth-szv'
+cell.pseudo = 'gth-hf-rev'
+cell.a = '''
+0.000000000, 3.370137329, 3.370137329
+3.370137329, 0.000000000, 3.370137329
+3.370137329, 3.370137329, 0.000000000'''
+cell.unit = 'B'
+cell.verbose = 0
+cell.build()
+
+kmesh = [1, 1, 3]
+kpts = cell.make_kpts(kmesh)
+nkpts = len(kpts)
+mf = scf.KRHF(cell, kpts).rs_density_fit()
+mf.kernel()
+```
+
+then resource estimates for the SF, DF, and THC factorization schemes can be generated a range of cutoffs. For example:
+
+```python
+from openfermion.resource_estimates.pbc import sf
+
+costs = sf.generate_costing_table(mf, name='carbon_diamond', naux_cutoffs=[20,25,30,35,40,45,50])
+print(costs.to_string(index=False))
+```
+will generate a `pandas.DataFrame` of resource estimates for the single factorization Hamiltonian (`sf`) and MP2 correlation energies for the range of auxiliary dimension (`naux_cutoffs`).
+
+
+Note that the automated costing computes the MP2 correlation energy error as a reference point for monitoring the convergence of the factorization with respect to sparsity or the size of auxiliary dimension. MP2 may be a poor model chemistry depending on the system, changing the option `energy_method = "CCSD"` will use CCSD instead, but this may become too expensive as the system size grows.
+
+The philosophy is that all costing methods are captured in the namespace related to the type of factorization. So if one wanted to repeat the costing for DF or THC factorizations, one could do 
+
+```python
+from openfermion.resource_estimates.pbc import df, thc
+
+# We need to specify eigenvalue threshold for second factorization.
+df_cutoffs = [1e-2, 5e-3, 1e-3, 5e-4, 1e-4, 5e-5, 1e-5]
+df_table = df.generate_costing_table(
+    mf,
+    name='carbon-diamond',
+    cutoffs=df_cutoffs)
+
+# Specify THC rank parameter we wish to scan over.
+# Note THC dimension M = thc_rank_param * N, N = number of spin orbitals in the
+# unit cell. 
+thc_rank_params = [2, 4, 6]
+# if you want to save each THC result to a file, you can set 'save_thc' to True
+thc_table = thc.generate_costing_table(mf, name='carbon-diamond', thc_rank_params=thc_rank_params) 
+```
+
+More fine-grained control is given by subroutines that compute the factorization, the lambda values, and the cost estimates.
+Further details are provided in a [tutorial](./notebooks/resource_estimates.ipynb). Note the THC factorization is more involved and we refer the reader to [thc-tutorial](../../notebooks/isdf.ipynb) for further details.
+
+Similar to the case of molecular resource estimation, we do not wish to burden all OpenFermion users with these dependencies, and testing with GitHub workflows is disabled. Currently we only check if pyscf is available. If it is then pytest will pick up the pbc module and run the tests. Note the tests can be quite slow due to the cost associated with building the integrals.
+
+## Requirements
+Requirements can be found in [resource_estimates.txt](../../../../dev_tools/requirements/deps/resource_estimates.txt)
+```
+pyscf
+jax
+jaxlib
+ase
+```
\ No newline at end of file
diff --git a/src/openfermion/resource_estimates/pbc/df/__init__.py b/src/openfermion/resource_estimates/pbc/df/__init__.py
index 53f003f..70a83ae 100644
--- a/src/openfermion/resource_estimates/pbc/df/__init__.py
+++ b/src/openfermion/resource_estimates/pbc/df/__init__.py
@@ -9,3 +9,11 @@
 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from .compute_lambda_df import compute_lambda
+    from .generate_costing_table_df import generate_costing_table
+    from .df_integrals import DFABKpointIntegrals
+    from .compute_df_resources import compute_cost
diff --git a/src/openfermion/resource_estimates/pbc/df/compute_df_resources.py b/src/openfermion/resource_estimates/pbc/df/compute_df_resources.py
new file mode 100644
index 0000000..43379dc
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/df/compute_df_resources.py
@@ -0,0 +1,342 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Determine costs for DF decomposition"""
+from typing import Tuple, Union, List
+import numpy as np
+from numpy.lib.scimath import arccos, arcsin
+from sympy import factorint
+
+from openfermion.resource_estimates.utils import QI
+
+from openfermion.resource_estimates.pbc.resources import (
+    ResourceEstimates,
+    QR3,
+)
+
+
+def compute_beta_for_resources(num_spin_orbs: int, num_kpts: int,
+                               de_for_qpe: float):
+    """Compute beta (number of bits for controlled rotations).
+
+    Uses expression from https://arxiv.org/pdf/2007.14460.pdf.
+
+    Args:
+        num_spin_orbs: Number of spin orbitals (per k-point)
+        num_kpts: Number of k-points.
+        de_for_qpe: epsilon for phase estimation.
+    """
+    return np.ceil(5.652 + np.log2(num_spin_orbs * num_kpts / de_for_qpe))
+
+
+def compute_cost(
+        num_spin_orbs: int,
+        lambda_tot: float,
+        num_aux: int,
+        num_eig: int,
+        kmesh: List[int],
+        dE_for_qpe: float = 0.0016,
+        chi: int = 10,
+        beta: Union[int, None] = None,
+) -> ResourceEstimates:
+    """Determine fault-tolerant costs using double factorized Hamiltonian.
+
+    Light wrapper around _compute_cost to automate choice of steps parameter.
+
+    Arguments:
+        num_spin_orbs: the number of spin-orbitals
+        lambda_tot: the lambda-value for the Hamiltonian
+        num_aux: number auxiliary vectors for the first factorization.
+        num_eig: number of retained eigenvalues for the second factorization.
+        dE_for_qpe: allowable error in phase estimation
+        chi: the number of bits for the representation of the coefficients
+        beta: the number of bits for controlled rotations
+    Returns:
+        resources: double factorized resource estimates
+    """
+    # run once to determine steps parameter
+    if beta is None:
+        num_kpts = int(np.prod(kmesh))
+        beta = compute_beta_for_resources(num_spin_orbs, num_kpts, dE_for_qpe)
+    init_cost = _compute_cost(
+        n=num_spin_orbs,
+        lam=lambda_tot,
+        dE=dE_for_qpe,
+        L=num_aux,
+        Lxi=num_eig,
+        chi=chi,
+        beta=beta,
+        Nkx=kmesh[0],
+        Nky=kmesh[1],
+        Nkz=kmesh[2],
+        steps=20_000,
+    )
+    steps = init_cost[0]
+    final_cost = _compute_cost(
+        n=num_spin_orbs,
+        lam=lambda_tot,
+        dE=dE_for_qpe,
+        L=num_aux,
+        Lxi=num_eig,
+        chi=chi,
+        beta=beta,
+        Nkx=kmesh[0],
+        Nky=kmesh[1],
+        Nkz=kmesh[2],
+        steps=steps,
+    )
+    estimates = ResourceEstimates(
+        toffolis_per_step=final_cost[0],
+        total_toffolis=final_cost[1],
+        logical_qubits=final_cost[2],
+    )
+    return estimates
+
+
+def _compute_cost(
+        n: int,
+        lam: float,
+        dE: float,
+        L: int,
+        Lxi: int,
+        chi: int,
+        beta: int,
+        Nkx: int,
+        Nky: int,
+        Nkz: int,
+        steps: int,
+        verbose: bool = False,
+) -> Tuple[int, int, int]:
+    """Determine fault-tolerant costs using DF decomposition.
+
+    Arguments:
+        n: the number of spin-orbitals.
+        lam: the lambda-value for the Hamiltonian
+        dE: allowable error in phase estimation
+        L: the rank of the first decomposition
+        Lxi: the total number of eigenvectors
+        chi: the number of bits for the representation of the coefficients
+        beta: the number of bits for the rotations
+        Nkx: Number of k-points in x-direction
+        Nky: Number of k-points in y-direction
+        Nkz:  Number of k-points in z-direction
+        steps: an approximate number of steps to choose the precision of
+            single qubit rotations in preparation of the equal superposition
+            state.
+        verbose: do additional printing of intermediates?
+
+    Returns:
+        step_cost: Toffolis per step
+        total_cost: Total number of Toffolis
+        ancilla_cost: Total ancilla cost
+    """
+    nNk = (max(np.ceil(np.log2(Nkx)), 1) + max(np.ceil(np.log2(Nky)), 1) +
+           max(np.ceil(np.log2(Nkz)), 1))
+    Nk = Nkx * Nky * Nkz
+
+    # The number of bits used for the second register.
+    nxi = np.ceil(np.log2(Nk * n // 2))
+
+    # The number of bits for the contiguous register.
+    nLxi = np.ceil(np.log2(Lxi + Nk * n // 2))
+
+    # The number of bits used for the first register.
+    nL = np.ceil(np.log2(L + 1))
+
+    # The power of 2 that is a factor of L + 1
+    if L % 2 == 1:
+        factors = factorint(L + 1)
+        eta = factors[min(list(sorted(factors.keys())))]
+    else:
+        eta = 0
+
+    oh = [0] * 20
+    for p in range(20):
+        v = np.round(
+            np.power(2, p + 1) / (2 * np.pi) *
+            arccos(np.power(2, nL) / np.sqrt((L + 1) / 2**eta) / 2))
+        oh[p] = np.real(steps * (1 / (np.sin(3 * arcsin(
+            np.cos(v * 2 * np.pi / np.power(2, p + 1)) * np.sqrt(
+                (L + 1) / 2**eta) / np.power(2, nL)))**2) - 1) + 4 * (p + 1))
+
+    # Bits of precision for rotation
+    br = int(np.argmin(oh) + 1)
+
+    # The following costs are from the list starting on page 50.
+
+    # The cost for preparing an equal superposition for preparing the first
+    # register in step 1 (a). We double this cost to account for the inverse.
+    cost1a = 2 * (3 * nL + 2 * br - 3 * eta - 9)
+
+    # The output size for the QROM for the first state preparation in Eq. (C27)
+    bp1 = nL + chi
+
+    # The cost of the QROM for the first state preparation in step 1 (b) and
+    # its inverse.
+    cost1b = QR3(L + 1, bp1)[1] + QI(L + 1)[1]
+
+    # The cost for the inequality test, controlled swap and their inverse in
+    # steps 1 (c) and (d)
+    cost1cd = 2 * (chi + nL)
+
+    # The total cost for preparing the first register in step 1.
+    cost1 = cost1a + cost1b + cost1cd
+
+    # The output size for the QROM for the data to prepare the equal
+    # superposition on the second register, as given in Eq. (C29).
+    bo = nxi + nLxi + br + 1
+    bo = bo + nNk  # account for k-point sampling outputting values of Q.
+
+    # This is step 2. This is the cost of outputting the data to prepare the
+    # equal superposition on the second register. We will assume it is not
+    # uncomputed, because we want to keep the offset for applying the QROM for
+    # outputting the rotations.
+    cost2 = QR3(L + 1, bo)[1] + QI(L + 1)[1]
+
+    # The number of bits for rotating the ancilla for the second preparation.
+    # We are just entering this manually because it is a typical value.
+    br = 7
+
+    # The cost of preparing an equal superposition over the second register in
+    # a controlled way. We pay this cost 4 times.
+    cost3a = 4 * (7 * nxi + 2 * br - 6)
+
+    # The cost of the offset to apply the QROM for state preparation on the
+    # second register.
+    cost3b = 4 * (nLxi - 1)
+
+    bp2 = nxi + chi + 2
+
+    # The cost of the QROMs and inverse QROMs for the state preparation, where
+    # in the first one we need + n/2 to account for the one-electron terms.
+    cost3c = (QR3(Lxi + Nk * n // 2, bp2)[1] + QI(Lxi + Nk * n // 2)[1] +
+              QR3(Lxi, bp2)[1] + QI(Lxi)[1])
+
+    # The inequality test and state preparations.
+    cost3d = 4 * (nxi + chi)
+
+    # The total costs for state preparations on register 2.
+    cost3 = cost3a + cost3b + cost3c + cost3d
+
+    # The cost of adding offsets in steps 4 (a) and (h).
+    cost4ah = 4 * (nLxi - 1)
+
+    # The costs of the QROMs and their inverses in steps 4 (b) and (g).
+    cost4bg = (QR3(Lxi + Nk * n // 2, 2 * n * beta + nNk)[1] +
+               QI(Lxi + Nk * n // 2)[1] + QR3(Lxi, 4 * n)[1] + QI(Lxi)[1])
+
+    # The cost of the controlled swaps based on the spin qubit in steps 4c and f
+    cost4cf = 2 * n * Nk
+
+    # The controlled rotations in steps 4 (d) and (f).
+    cost4df = 16 * n * (beta - 2)
+
+    # The controlled Z operations in the middle for step 4 (e).
+    cost4e = 3
+    cost4e = (
+        cost4e + 4 * n * Nk + 12 * nNk
+    )  # extra cost of swapping into working registers and computing k-Q.
+
+    # (*Adjustment for modular arithmetic costs.*)
+    if Nkx == 2**np.ceil(np.log2(Nkx)):
+        cost4e = cost4e - 2 * np.ceil(np.log2(Nkx))
+    if Nky == 2**np.ceil(np.log2(Nky)):
+        cost4e = cost4e - 2 * np.ceil(np.log2(Nky))
+    if Nkz == 2**np.ceil(np.log2(Nkz)):
+        cost4e = cost4e - 2 * np.ceil(np.log2(Nkz))
+
+    # This is the cost of the controlled rotations for step 4.
+    cost4 = cost4ah + cost4bg + cost4cf + cost4df + cost4e
+
+    # This is the cost of the reflection on the second register from step 6.
+    cost6 = nxi + chi + 2
+
+    # The cost of the final reflection req'd to construct the step of the
+    # quantum walk from step 9.
+    cost9 = nL + nxi + chi + 1
+
+    # The extra two qubits for unary iteration and making the rflxn controlled.
+    cost10 = 2
+
+    # The Toffoli cost for a single step
+    cost = cost1 + cost2 + cost3 + cost4 + cost6 + cost9 + cost10
+
+    # The number of steps needed
+    iters = np.ceil(np.pi * lam / (2 * dE))
+
+    # Now the number of qubits from the list on page 54.
+
+    k1 = np.power(2, QR3(Lxi + Nk * n // 2, 2 * n * beta + nNk)[0])
+
+    # The control register for phase estimation and iteration on it.
+    ac1 = np.ceil(np.log2(iters + 1)) * 2 - 1
+
+    # The system qubits
+    ac2 = n * Nk
+
+    # The first register prepared, a rotated qubit and a flag qubit.
+    ac3 = nL + 2
+
+    # The output of the QROM, the equal superposition state and a flag qubit.
+    ac4 = nL + chi * 2 + 1
+
+    # The data used for preparing the equal superposition state on the 2nd reg
+    ac5 = bo
+
+    # The second register, a rotated qubit and a flag qubit.
+    ac6 = nxi + 2
+
+    # The second preparation QROM output.
+    ac8 = bp2
+
+    # The equal superposition state and the result of the inequality test.
+    ac9 = chi + 1
+
+    # The angles for rotations.
+    ac10 = k1 * (2 * n * beta + nNk)
+
+    # The phase gradient state.
+    ac11 = beta
+
+    # A control qubit for the spin.
+    ac12 = 1
+
+    # A T state.
+    ac13 = 1
+
+    # extra for kpoint
+    ac14 = nNk + 2 * n
+
+    if verbose:
+        print("[*] Top of routine")
+        print("  [+] nxi = ", nxi)
+        print("  [+] nLxi = ", nLxi)
+        print("  [+] nL = ", nL)
+        print("  [+] eta = ", eta)
+        print("  [+] cost3 = ", cost3)
+        print("  [+] cost4 = ", cost4)
+        print("  [+] cost = ", cost)
+        print("  [+] iters = ", iters)
+
+    ancilla_cost = (ac1 + ac2 + ac3 + ac4 + ac5 + ac6 + ac8 + ac9 + ac10 +
+                    ac11 + ac12 + ac13 + ac14)
+
+    # Sanity checks before returning as int
+    assert cost.is_integer()
+    assert iters.is_integer()
+    assert ancilla_cost.is_integer()
+
+    step_cost = int(cost)
+    total_cost = int(cost * iters)
+    ancilla_cost = int(ancilla_cost)
+
+    return step_cost, total_cost, ancilla_cost
diff --git a/src/openfermion/resource_estimates/pbc/df/generate_costing_table_df.py b/src/openfermion/resource_estimates/pbc/df/generate_costing_table_df.py
new file mode 100644
index 0000000..f132eb3
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/df/generate_costing_table_df.py
@@ -0,0 +1,148 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+from dataclasses import dataclass
+from functools import reduce
+import pandas as pd
+
+import numpy as np
+
+from pyscf.pbc import scf
+from pyscf.pbc.tools.k2gamma import kpts_to_kmesh
+
+from openfermion.resource_estimates.pbc.resources import PBCResources
+from openfermion.resource_estimates.pbc.df.df_integrals import (
+    DFABKpointIntegrals,)
+from openfermion.resource_estimates.pbc.hamiltonian import (
+    build_hamiltonian,)
+from openfermion.resource_estimates.pbc.hamiltonian.cc_extensions import (
+    build_approximate_eris,
+    build_cc_inst,
+    build_approximate_eris_rohf,
+)
+from openfermion.resource_estimates.pbc.df.compute_lambda_df import (
+    compute_lambda,)
+from openfermion.resource_estimates.pbc.df.compute_df_resources import (
+    compute_cost, compute_beta_for_resources)
+
+
+@dataclass
+class DFResources(PBCResources):
+    """Resource estimates for DF factorization.
+
+    Inherits from PBCReousrces data class and adds the additional attributes:
+
+    Attributes:
+        num_aux: number of auxiliary vectors for first factorization. By default
+            this is not truncated.
+        beta: The number of bits for the controlled rotations.
+    """
+    num_aux: int = -1
+    beta: int = 20
+
+
+def generate_costing_table(
+        pyscf_mf: scf.HF,
+        cutoffs: np.ndarray,
+        name: str = "pbc",
+        chi: int = 10,
+        beta: int = 20,
+        dE_for_qpe: float = 0.0016,
+        energy_method: str = "MP2",
+) -> pd.DataFrame:
+    """Generate resource estimate costing table given a set of cutoffs for
+        double-factorized Hamiltonian.
+
+    Arguments:
+        pyscf_mf: k-point pyscf mean-field object
+        cutoffs: Array of (integer) auxiliary index cutoff values
+        name: Optional descriptive name for simulation.
+        chi: the number of bits for the representation of the coefficients
+        beta: the number of bits for controlled rotation. If None we estimate
+            the value given N and dE_for_qpe. See compute_beta_for_resources for
+            details.
+        dE_for_qpe: Phase estimation epsilon.
+        energy_method: Which model chemistry to use to estimate convergence of
+            factorization with respect to threshold. values are MP2 or CCSD.
+
+    Returns
+        resources: Table of resource estimates.
+
+    Raises:
+        ValueError if energy_method is unknown.
+    """
+    kmesh = kpts_to_kmesh(pyscf_mf.cell, pyscf_mf.kpts)
+    cc_inst = build_cc_inst(pyscf_mf)
+    exact_eris = cc_inst.ao2mo()
+    if energy_method.lower() == "mp2":
+        energy_function = lambda x: cc_inst.init_amps(x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    elif energy_method.lower() == "ccsd":
+        energy_function = lambda x: cc_inst.kernel(eris=x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    else:
+        raise ValueError(f"Unknown value for energy_method: {energy_method}")
+
+    hcore, chol = build_hamiltonian(pyscf_mf)
+    num_spin_orbs = 2 * hcore[0].shape[-1]
+    num_kpts = np.prod(kmesh)
+
+    num_aux = chol[0, 0].shape[0]
+    # This is constant as we don't truncate first factorization for DF
+    num_aux_df = 2 * num_aux * num_kpts
+    df_resource_obj = DFResources(
+        system_name=name,
+        num_spin_orbitals=num_spin_orbs,
+        num_kpts=num_kpts,
+        dE=dE_for_qpe,
+        chi=chi,
+        beta=beta,
+        energy_method=energy_method,
+        exact_energy=np.real(reference_energy),
+        num_aux=num_aux_df,
+    )
+    # Save some space and overwrite eris object from exact CC
+    approx_eris = exact_eris
+    if beta is None:
+        num_kpts = np.prod(kmesh)
+        beta = compute_beta_for_resources(num_spin_orbs, num_kpts, dE_for_qpe)
+    for cutoff in cutoffs:
+        df_helper = DFABKpointIntegrals(cholesky_factor=chol, kmf=pyscf_mf)
+        df_helper.double_factorize(thresh=cutoff)
+        df_lambda = compute_lambda(hcore, df_helper)
+        if pyscf_mf.cell.spin == 0:
+            approx_eris = build_approximate_eris(cc_inst,
+                                                 df_helper,
+                                                 eris=approx_eris)
+        else:
+            approx_eris = build_approximate_eris_rohf(cc_inst,
+                                                      df_helper,
+                                                      eris=approx_eris)
+        approx_energy, _, _ = energy_function(approx_eris)
+        df_res_cost = compute_cost(
+            num_spin_orbs,
+            df_lambda.lambda_total,
+            num_aux_df,
+            df_lambda.num_eig,
+            list(kmesh),
+            chi=chi,
+            beta=beta,
+            dE_for_qpe=dE_for_qpe,
+        )
+        df_resource_obj.add_resources(
+            ham_properties=df_lambda,
+            resource_estimates=df_res_cost,
+            cutoff=cutoff,
+            approx_energy=np.real(approx_energy),
+        )
+
+    return df_resource_obj.to_dataframe()
diff --git a/src/openfermion/resource_estimates/pbc/notebooks/isdf.ipynb b/src/openfermion/resource_estimates/pbc/notebooks/isdf.ipynb
index 7c194c9..014c021 100644
--- a/src/openfermion/resource_estimates/pbc/notebooks/isdf.ipynb
+++ b/src/openfermion/resource_estimates/pbc/notebooks/isdf.ipynb
@@ -112,7 +112,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from openfermion.resource_estimates.pbc.utils import build_momentum_transfer_mapping \n",
+    "from openfermion.resource_estimates.pbc.hamiltonian import build_momentum_transfer_mapping \n",
     "\n",
     "momentum_map = build_momentum_transfer_mapping(cell, kpts)\n",
     "num_spatial_orbs = mf.mo_coeff[0].shape[-1]\n",
@@ -154,7 +154,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from openfermion.resource_estimates.pbc.utils import compute_emp2_approx, build_cc_inst\n",
+    "from openfermion.resource_estimates.pbc.hamiltonian.cc_extensions import compute_emp2_approx, build_cc_inst\n",
     "cc_inst = build_cc_inst(mf)\n",
     "eri_exact = cc_inst.ao2mo()\n",
     "emp2_exact, _, _ = cc_inst.init_amps(eri_exact)\n",
@@ -165,7 +165,7 @@
     "for cthc in thc_ranks:\n",
     "    num_thc = cthc * num_spatial_orbs\n",
     "    kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
-    "    kthc_eri_helper = kthc.KPTHCHelperDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
+    "    kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
     "    eri_thc = kthc_eri_helper.get_eri(eri_kindices)\n",
     "    eri_exact = kthc_eri_helper.get_eri_exact(eri_kindices)\n",
     "    # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
@@ -293,7 +293,7 @@
     "        print(f\"Running mesh = {mesh}, cthc = {cthc}\")\n",
     "        num_thc = cthc * num_spatial_orbs\n",
     "        kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
-    "        kthc_eri_helper = kthc.KPTHCHelperDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
+    "        kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
     "        emp2_approx = compute_emp2_approx(mf, kthc_eri_helper)\n",
     "        # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
     "        results[mesh].append(emp2_exact - emp2_approx)\n",
@@ -354,7 +354,7 @@
     "        print(f\"Running basis = {basis}, cthc = {cthc}\")\n",
     "        num_thc = cthc * num_spatial_orbs\n",
     "        kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
-    "        kthc_eri_helper = kthc.KPTHCHelperDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
+    "        kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
     "        emp2_approx = compute_emp2_approx(mf, kthc_eri_helper)\n",
     "        # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
     "        basis_results[basis].append(emp2_exact - emp2_approx)\n",
@@ -430,7 +430,7 @@
     "for cthc in thc_ranks:\n",
     "    num_thc = cthc * num_spatial_orbs\n",
     "    kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
-    "    kthc_eri_helper = kthc.KPTHCHelperDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
+    "    kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
     "    # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
     "    thc_lambda = kthc.compute_lambda(hcore, kthc_eri_helper)\n",
     "    emp2_approx = compute_emp2_approx(mf, kthc_eri_helper) - mf.e_tot\n",
@@ -469,7 +469,7 @@
     "# Here we use the helper function kpoint_thc_via_isdf which will first find\n",
     "# the ISDF factors and feed these into the BFGS and AdaGrad solvers.\n",
     "kpt_thc, _ = kthc.kpoint_thc_via_isdf(mf, Luv, num_thc, verbose=False)\n",
-    "kthc_eri_helper = kthc.KPTHCHelperDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
+    "kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
     "# Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
     "thc_lambda = kthc.compute_lambda(hcore, kthc_eri_helper)\n",
     "emp2_approx = utils.compute_emp2_approx(mf, kthc_eri_helper) - mf.e_tot # only compare correlation energy.\n",
@@ -494,7 +494,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.2"
+   "version": "3.11.4"
   },
   "orig_nbformat": 4,
   "vscode": {
diff --git a/src/openfermion/resource_estimates/pbc/notebooks/resource_estimates.ipynb b/src/openfermion/resource_estimates/pbc/notebooks/resource_estimates.ipynb
new file mode 100644
index 0000000..20b0fce
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/notebooks/resource_estimates.ipynb
@@ -0,0 +1,546 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Estimating Fault Tolerant Resources for Periodic Systems\n",
+    "\n",
+    "The resource estimation code provided in this module relies on `pyscf` to compute the required symmetry adapted molecular orbitals, one- and two-electron integrals, and correlated wavefunction calculations to determine the accuracy of different truncation schemes. So, to start the tutorial let's run a periodic restricted Hartree-Fock (RHF) simulation of diamond. \n",
+    "\n",
+    "We assume a knowledge of electronic structure theory for solids, and in particular an understanding that any results need to be converged to the infinite $k$-point (thermodynamic) and complete basis set limit. For the purposes of this tutorial we will focus on results far away from this limit and simulate simple systems in minimal basis sets use very small $k$-point meshes. The module also assumes the use of density fitted integrals and we use range-separated density fitting throughout. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "plt.rcParams['figure.dpi'] = 130 \n",
+    "\n",
+    "from ase.build import bulk\n",
+    "import numpy as np\n",
+    "\n",
+    "from pyscf.pbc import gto, scf\n",
+    "from pyscf.pbc.tools import pyscf_ase\n",
+    "\n",
+    "\n",
+    "# Build a 2 atom unit cell for carbon in the diamond structure near it's\n",
+    "# equilibrium lattice constant.\n",
+    "ase_atom = bulk(\"C\", \"diamond\", a=3.5)\n",
+    "cell = gto.Cell()\n",
+    "cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)\n",
+    "cell.a = ase_atom.cell[:].copy()\n",
+    "# Using a minimal basis set for expediency.\n",
+    "cell.basis = \"gth-szv\"\n",
+    "cell.pseudo = \"gth-hf-rev\"\n",
+    "cell.verbose = 0\n",
+    "cell.build()\n",
+    "\n",
+    "# We are using a very small k-point mesh for speed purposes too.\n",
+    "kmesh = [1, 1, 3]\n",
+    "kpts = cell.make_kpts(kmesh)\n",
+    "num_kpts = len(kpts)\n",
+    "mf = scf.KRHF(cell, kpts).rs_density_fit()\n",
+    "mf.kernel()\n",
+    "print(\"SCF energy: \", mf.e_tot)\n",
+    "\n",
+    "# converged SCF energy: -10.39193609748544"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Armed with our SCF solution we can now generate the one and two-electron integrals required to compute $\\lambda$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from openfermion.resource_estimates.pbc.hamiltonian import build_hamiltonian\n",
+    "# Get molecular orbital \"Cholesky\" integrals from RSGDF object, these are just\n",
+    "# 3-centre integrals (X|pq).\n",
+    "hcore_mo, chol = build_hamiltonian(mf)\n",
+    "print(\"(nkpts, nmo, nmo) = {}\".format(hcore_mo.shape))\n",
+    "print(\"(nkpts, nkpts) = {}\".format(chol.shape))\n",
+    "print(\"(naux, nmo, nmo) = {}\".format(chol[0,0].shape))\n",
+    "num_mo = hcore_mo.shape[-1]\n",
+    "num_aux = chol[0,0].shape[0]"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Hamiltonian Representation\n",
+    "The resource estimation module provides four representations for the Hamiltonian: sparse, single-factorization (SF), double-factorization (SF), and tensor hypercontraction (THC). Each of these approaches introduces a different parameter which controls the accuracy of the factorization, and results need to be monitored with respect to these truncation parameters.\n",
+    "\n",
+    "## Sparse Hamiltonian\n",
+    "\n",
+    "The sparse Hamiltonian takes the usual form for a second-quantized $k$-point dependent Hamiltonian\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "H = \\sum_{pq\\mathrm{k}} h_{pq}(\\mathrm{k}) a_{p\\mathbf{k}}^{\\dagger} a_{q\\mathbf{k}} + \\frac{1}{2} \\sum_{\\mathbf{k}_p\\mathbf{k}_q\\mathbf{k}_r\\mathbf{k}_s}\\sum_{pqrs} (p\\mathbf{k}_pq\\mathbf{k}_q|r\\mathbf{k}_rs\\mathbf{k}_s)  a_{p\\mathbf{k}_p}^{\\dagger} a_{r\\mathbf{k}_r}^{\\dagger} a_{s\\mathbf{k}_s}a_{q\\mathbf{k}_q} \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "Utilizing conservation of crystal momentum $\\mathbf{k}_p + \\mathbf{k}_r - \\mathbf{k}_q -\\mathbf{k}_s = \\mathbf{G}$, where $\\mathbf{G}$ is a reciprocal lattice vector, we can write \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "H = \\sum_{pq\\mathrm{k}} h_{pq}(\\mathrm{k}) a_{p\\mathbf{k}}^{\\dagger} a_{q\\mathbf{k}} + \\frac{1}{2} \\sum_{\\mathbf{Q}\\mathbf{k}\\mathbf{k}'}\\sum_{pqrs} (p\\mathbf{k}q\\mathbf{k}-\\mathbf{Q}|r\\mathbf{k}'-\\mathbf{Q}s\\mathbf{k}')  a_{p\\mathbf{k}}^{\\dagger} a_{r\\mathbf{k}'-\\mathbf{Q}}^{\\dagger} a_{s\\mathbf{k}'}a_{q\\mathbf{k}-\\mathbf{Q}} \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "where $\\mathbf{Q}$ is the momentum transfer vector which we chose to live in our set of $k$-points. Note the subtraction in the above expression is really modulo a $\\mathbf{G}$ vector.\n",
+    "\n",
+    "The Hamiltonian above has $N_k^3 N^4$ terms, where $N_k$ is the number of $k$-points and $N$ is the number of spin orbitals. The sparse representation attempts to approximate the Hamiltonian by zeroing elements of $H$ which are below some threshold. This will yield $\\mathcal{O}(s N_k^3 N^4)$ terms where $s$ is a sparsity factor.\n",
+    "\n",
+    "We provide helper functions that will sparsify the Hamiltonian to aid in resource estimation:\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from openfermion.resource_estimates.pbc import sparse\n",
+    "sparse_ham = sparse.SparseFactorization(chol, mf, threshold=1e-3)\n",
+    "# look at eri block\n",
+    "kpts = [0]*4\n",
+    "eri_approx = sparse_ham.get_eri(kpts)\n",
+    "eri_exact = sparse_ham.get_eri_exact(kpts)\n",
+    "# With a sparsity threshold of 1e-3, the approximate eris and \"exact\" eris (i.e.\n",
+    "# those with no truncation)  should yield different results.\n",
+    "assert not np.allclose(eri_approx, eri_exact)\n",
+    "print(\"Total number of elements N_k^3 N^4 = {:d}\".format(num_kpts**3*num_mo**4))\n",
+    "print(\"number of symmetry unique non zero = {}\".format(sparse_ham.get_total_unique_terms_above_thresh()))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With the Hamiltonian at hand we can compute the 1-norm (called $\\lambda$) of the Hamiltonian which is essential for computing resource estimates:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sparse_lambda = sparse.compute_lambda(hcore_mo, sparse_ham)\n",
+    "print(sparse_lambda)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally we can compute the total resource costs. In particular the following code will compute the number of Toffoli gates required for a single step of phase estimation (`toffolis_per_step`), the total Toffoli count (`total_toffolis`), and the number of logical qubits (`logical_qubits`).\n",
+    "\n",
+    "The total Toffoli count is given by the toffoli per step cost multiplied by the factor \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "\\left\\lceil \\frac{\\pi \\lambda}{2 \\epsilon_{\\mathrm{QPE}}} \\right\\rceil\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "where $\\epsilon_{\\mathrm{QPE}}$ is our variable `dE_for_qpe`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_spin_orbs = 2*num_mo\n",
+    "resources = sparse.compute_cost(num_spin_orbs, sparse_lambda.lambda_total, sparse_lambda.num_sym_unique , kmesh, dE_for_qpe=0.0016)\n",
+    "print(resources)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We see there are the following steps required:\n",
+    "\n",
+    "1. Run an SCF calculation.\n",
+    "2. Generate the one- and two-electron matrix elements.\n",
+    "3. Compute the lambda value of the Hamiltonian. \n",
+    "4. Compute the resource estimates.\n",
+    "\n",
+    "These steps are required for all four factorizations and involve a lot of boilerplate code. As such we provide utility functions which perform these necessary steps, and scan over the value of the corresponding threshold parameter computing either the MP2 or CCSD correlation energy for each threshold value.\n",
+    "\n",
+    "Let's see how this works for the sparse factorization first."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "thresholds = np.logspace(-1, -5, 5)\n",
+    "sparse_costing_table = sparse.generate_costing_table(mf, thresholds=thresholds)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The costing table is a `pandas.DataFrame` which can be convenient for saving to `csv` `(to_csv)`, string `(to_sting())` or outputting results to $\\LaTeX$ (`to_latex()`)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(sparse_costing_table.to_string(index=False))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A crucial step in resource estimation is determining when the Hamiltonian is sufficiently accurate given a certain truncation. One way to check this is to use a correlated wavefunction method and monitor the convergence of the correlation energy with the truncation parameter. Let's look at the convergence of the MP2 error with the sparsity threshold. Note that MP2 is not expected to be a particularly faithful model to monitor convergence as it only requires a subset of the integral blocks. CCSD is a better option but comes with a considerable overhead."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(sparse_costing_table.cutoff, np.abs(sparse_costing_table.approx_energy-sparse_costing_table.exact_energy), marker=\"o\")\n",
+    "plt.xscale(\"log\")\n",
+    "plt.yscale(\"log\")\n",
+    "plt.xlabel(\"sparse threshold\")\n",
+    "plt.ylabel(\"MP2 Energy Error (Ha)\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(sparse_costing_table.cutoff, np.abs(sparse_costing_table.lambda_total), marker=\"o\")\n",
+    "plt.xscale(\"log\")\n",
+    "plt.xlabel(\"sparse threshold\")\n",
+    "plt.ylabel(\"$\\lambda$\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(sparse_costing_table.cutoff, np.abs(sparse_costing_table.toffolis_per_step), marker=\"o\")\n",
+    "plt.xscale(\"log\")\n",
+    "plt.xlabel(\"sparse threshold\")\n",
+    "plt.ylabel(\"Toffolis per step\")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Single Factorization \n",
+    "\n",
+    "The single factorization approach follows from the usual density fitting or Cholesky factorization of the two-electron integrals:\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "(p\\mathbf{k}q\\mathbf{k}-\\mathbf{Q}|r\\mathbf{k}'-\\mathbf{Q}s\\mathbf{k}') = \\sum_n^{M} L_{p\\mathbf{k}q\\mathbf{k}-\\mathbf{Q}}^n L_{s\\mathbf{k}'r\\mathbf{k}'-\\mathbf{Q}}^{n*} \n",
+    "\n",
+    "\n",
+    "$$\n",
+    "where $M$ is the dimension of an auxiliary index $n$. With this factorization we can define \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "\\hat{\\rho}_{n}(\\mathbf{Q}, \\mathbf{k}) = \\left(\\sum_{\\sigma \\in \\{\\uparrow, \\downarrow\\}}\\sum_{pq}^{N/2}L_{p \\mathbf{k} q (\\mathbf{k}-\\mathbf{Q}), n} a_{p\\mathbf{k}\\sigma}^\\dagger a_{q(\\mathbf{k}-\\mathbf{Q})\\sigma}  \\right)  ,\\qquad \\hat{\\rho}^\\dagger_{n}(\\mathbf{Q}, \\mathbf{k}) = \\left( \\sum_{\\sigma \\in \\{\\uparrow, \\downarrow\\}} \\sum_{pq}^{N/2}L^{*}_{p \\mathbf{k} q (\\mathbf{k}-\\mathbf{Q}), n} a_{q(\\mathbf{k}-\\mathbf{Q})\\sigma}^\\dagger a_{p\\mathbf{k}\\sigma}\\right) \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "to write\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "\\hat{H}_2 = \\frac{1}{2}  \\sum_{\\mathbf{Q}}^{N_{k}}\\sum_{n}^{M} \\left(\\hat{A}^2_{n}(\\mathbf{Q}) + \\hat{B}^2_{n}(\\mathbf{Q})\\right) + \\mathrm{one\\ body\\ term},\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "where\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "\\hat{A}_{n}(\\mathbf{Q}) =\\frac{1}{2}(\\hat{\\rho}_{n}(\\mathbf{Q}) + \\hat{\\rho}^\\dagger_{n}(\\mathbf{Q})),\\\\\n",
+    "\\hat{B}_{n}(\\mathbf{Q}) = \\frac{i}{2}(\\hat{\\rho}_{n}(\\mathbf{Q}) - \\hat{\\rho}^\\dagger_{n}(\\mathbf{Q})).\n",
+    "\n",
+    "$$\n",
+    "\n",
+    "Given this representation of the Hamiltonian it remains to check how rapidly the factorization converges with $M$ typically called `num_aux` below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from openfermion.resource_estimates.pbc import sf \n",
+    "cutoffs = np.arange(10, num_aux, 10)\n",
+    "sf_costing_table = sf.generate_costing_table(mf, naux_cutoffs=cutoffs)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(sf_costing_table.to_string(index=None))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(sf_costing_table.cutoff, np.abs(sf_costing_table.approx_energy-sf_costing_table.exact_energy), marker=\"o\")\n",
+    "plt.yscale(\"log\")\n",
+    "plt.xlabel(\"$M$\")\n",
+    "plt.ylabel(\"MP2 Energy Error (Ha)\")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "## Double Factorization \n",
+    "\n",
+    "The double factorization follows on from the single-factorized Hamiltonian, and one performs a further factorization on the $A$ and $B$ operators defined above. \n",
+    "\n",
+    "$$\n",
+    "\n",
+    "\\hat{A}_{n}(\\mathbf{Q}) = \\sum_{\\mathbf{k}} \\left[ U^A_{n}(\\mathbf{Q}, \\mathbf{k}) \\left( \\sum_{\\sigma}\\sum_{p}^{\\Xi_{\\mathbf{Q}, n,\\mathbf{k},A}}f^A_{p}(\\mathbf{Q}, n, \\mathbf{k})n_{p\\mathbf{k}\\sigma} \\right) U^A_{n}(\\mathbf{Q}, \\mathbf{k})^{\\dagger} \\right]\n",
+    "\n",
+    "$$\n",
+    "where $U^A$ is a unitary that diagonalizes $A$, $f_p$ are the corresponding eigenvalues, and $n_{p\\mathbf{k}}$ is the number operator. Note that the sum over $\\mathbf{k}$ is factored out of the expression which differs from the single factorized case and is crucial to observe a $\\sqrt{N_k}$ speedup in the quantum algorithm. We thus need to monitor convergence with respect to $M$ and the number of eigenvalues $\\sum_{\\mathbf{Q},n,\\mathbf{k}}\\Xi_{\\mathbf{Q}, n,\\mathbf{k},A}$, typically called `num_eigs` below. Note that we fix $M$ to be the full rank and just monitor convergence w.r.t the second factorization for simplicity."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from openfermion.resource_estimates.pbc import df\n",
+    "cutoffs = np.logspace(-1, -5, 5)\n",
+    "df_costing_table = df.generate_costing_table(mf, cutoffs=cutoffs)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(df_costing_table.to_string(index=False))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.rcParams['figure.dpi'] = 130 \n",
+    "plt.plot(df_costing_table.cutoff, np.abs(df_costing_table.approx_energy-df_costing_table.exact_energy), marker=\"o\")\n",
+    "plt.yscale(\"log\")\n",
+    "plt.xscale(\"log\")\n",
+    "plt.xlabel(\"eigenvalue cutoff\")\n",
+    "plt.ylabel(\"MP2 Energy Error (Ha)\")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Tensor hypercontraction factorization\n",
+    "\n",
+    "The tensor hypercontraction factorization (THC) is somewhat involved and we refer the reader to the [isdf_notebook](./isdf.ipynb) for further details, however the procedure follows much like before (except it is considerably more expensive to generate the factorization). For expediency we will forego reoptimizing the THC factors in the example below and just look at the ISDF convergence. This may take a few minutes to run."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from openfermion.resource_estimates.pbc import thc \n",
+    "thc_rank_params = [2, 4, 6, 8] \n",
+    "thc_costing_table = thc.generate_costing_table(mf, thc_rank_params=thc_rank_params, reoptimize=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(thc_costing_table.to_string(index=False))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.rcParams['figure.dpi'] = 130 \n",
+    "plt.plot(thc_costing_table.cutoff, np.abs(thc_costing_table.approx_energy-thc_costing_table.exact_energy), marker=\"o\")\n",
+    "plt.yscale(\"log\")\n",
+    "plt.xlabel(\"THC rank parameter\")\n",
+    "plt.ylabel(\"MP2 Energy Error (Ha)\")"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Other Considerations\n",
+    "\n",
+    "\n",
+    "### Spin polarization\n",
+    "The discussion so far assumed we were performing a RHF calculation to generate spin-free integrals and a closed-shell mean field solution. While unrestricted solutions are not supported, the resource estimation code will work for ROHF solutions. This can be helpful if the system of interest has a different number of alpha and beta electrons.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ase.build import bulk\n",
+    "import numpy as np\n",
+    "\n",
+    "from pyscf.pbc import gto, scf\n",
+    "from pyscf.pbc.tools import pyscf_ase\n",
+    "# Build a 2 atom unit cell for carbon in the diamond structure near it's\n",
+    "# equilibrium lattice constant.\n",
+    "ase_atom = bulk(\"C\", \"diamond\", a=3.5)\n",
+    "cell = gto.Cell()\n",
+    "cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)\n",
+    "cell.a = ase_atom.cell[:].copy()\n",
+    "# Using a minimal basis set for expediency.\n",
+    "cell.basis = \"gth-szv\"\n",
+    "cell.pseudo = \"gth-hf-rev\"\n",
+    "cell.verbose = 0\n",
+    "cell.spin = 2  # Force high spin for ROHF\n",
+    "cell.build()\n",
+    "\n",
+    "# We are using a very small k-point mesh for speed purposes too.\n",
+    "kmesh = [1, 1, 3]\n",
+    "kpts = cell.make_kpts(kmesh)\n",
+    "num_kpts = len(kpts)\n",
+    "mf = scf.KRHF(cell, kpts).rs_density_fit()\n",
+    "mf.kernel()\n",
+    "print(\"SCF energy: \", mf.e_tot)\n",
+    "print(mf.cell.nelec)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutoffs = np.logspace(-1, -5, 5)\n",
+    "df_costing_table = df.generate_costing_table(mf, cutoffs=cutoffs)\n",
+    "print(df_costing_table.to_string(index=False))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "### Threshold convergence\n",
+    "\n",
+    "So far we have used MP2 to monitor the convergence of the factorization, however we can instead use CCSD as a better model chemistry."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "thresholds = np.logspace(-1, -5, 5)\n",
+    "sparse_costing_table = sparse.generate_costing_table(mf, thresholds=thresholds, energy_method=\"CCSD\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(sparse_costing_table.to_string(index=False))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "kpoint_eri",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.4"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/src/openfermion/resource_estimates/pbc/resources/__init__.py b/src/openfermion/resource_estimates/pbc/resources/__init__.py
new file mode 100644
index 0000000..1cf408b
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/resources/__init__.py
@@ -0,0 +1,17 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from .data_types import PBCResources, ResourceEstimates
+    from .qrom import QR2, QR3, QI2
diff --git a/src/openfermion/resource_estimates/pbc/resources/data_types.py b/src/openfermion/resource_estimates/pbc/resources/data_types.py
new file mode 100644
index 0000000..64f2f25
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/resources/data_types.py
@@ -0,0 +1,97 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Dataclasses for resource estimation."""
+from dataclasses import dataclass, asdict, field
+import pandas as pd
+
+from openfermion.resource_estimates.pbc.hamiltonian import (
+    HamiltonianProperties,)
+
+
+@dataclass(frozen=True)
+class ResourceEstimates:
+    """Dataclass to hold return values from compute_cost functions.
+
+    Attributes:
+        toffolis_per_step: Toffolis per step
+        total_toffolis: Total number of Toffolis
+        logical_qubits: Total ancilla cost
+    """
+
+    toffolis_per_step: int
+    total_toffolis: int
+    logical_qubits: int
+
+    dict = asdict
+
+
+@dataclass
+class PBCResources:
+    """Helper class to hold resource estimates given a range of cutoffs.
+
+    Attributes:
+        system_name: Descriptive name for calculation.
+        num_spin_orbitals: Number of spin orbitals.
+        num_kpts: Number of k-points.
+        dE: Epsilon for phase estimation.
+        chi: A number of bits controlling the precision of the factorization.
+            What this means is dependent on the factorization.
+        exact_energy: Exact energy (no cutoff).
+        energy_method: Which method to use to copute energy (MP2 or CCSD).
+        cutoff: list of cuttofs.
+        ham_props: list of lambda values for each cutoff.
+        resources: list of resource estimates for each cutoff.
+    """
+
+    system_name: str
+    num_spin_orbitals: int
+    num_kpts: int
+    dE: float
+    chi: int
+    exact_energy: float
+    energy_method: str = "MP2"
+    cutoff: list = field(default_factory=list)
+    approx_energy: list = field(default_factory=list)
+    ham_props: list = field(default_factory=list)
+    resources: list = field(default_factory=list)
+
+    dict = asdict
+
+    def to_dataframe(self) -> pd.DataFrame:
+        """Convert PBCResources instance to pandas DataFrame."""
+        df = pd.DataFrame(self.dict())
+        lambdas = pd.json_normalize(df.pop("ham_props"))
+        resources = pd.json_normalize(df.pop("resources"))
+        df = df.join(pd.DataFrame(lambdas))
+        df = df.join(pd.DataFrame(resources))
+        return df
+
+    def add_resources(
+            self,
+            ham_properties: HamiltonianProperties,
+            resource_estimates: ResourceEstimates,
+            cutoff: float,
+            approx_energy: float,
+    ) -> None:
+        """Add resource estimates to container for given cutoff value.
+
+        Args:
+            ham_properties: lambda values
+            resource_estimates: resource estimates for this value of cutoff
+            cutoff: Current factorization/representation cutoff value.
+            approx_energy: Current approximate energy estimate given cutoff.
+        """
+        self.ham_props.append(ham_properties)
+        self.resources.append(resource_estimates)
+        self.cutoff.append(cutoff)
+        self.approx_energy.append(approx_energy)
diff --git a/src/openfermion/resource_estimates/pbc/resources/qrom.py b/src/openfermion/resource_estimates/pbc/resources/qrom.py
new file mode 100644
index 0000000..92b674a
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/resources/qrom.py
@@ -0,0 +1,75 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Functions to determine optimal values for qro(a)m.
+
+These were adapted from corresponding Mathematica functions which are given in
+the docstrings.
+"""
+import numpy as np
+
+
+def QR3(L, M1):
+    r"""Find optimal cost for QROM based on L value of size M1.
+
+    QR[Ll_, m_] := Ceiling[MinValue[{Ll/2^k + m*(2^k - 1), k >= 0}, k
+    \[Element] Integers]];
+    """
+    k = 0.5 * np.log2(L / M1)
+    value = lambda k: L / np.power(2, k) + M1 * (np.power(2, k) - 1)
+    if k < 0:
+        k_opt = 0
+        val_opt = np.ceil(value(k_opt))
+        assert val_opt.is_integer()
+        return int(k_opt), int(val_opt)
+    k_int = [np.floor(k), np.ceil(k)]  # restrict optimal k to integers
+    k_opt = k_int[np.argmin(value(k_int))]  # obtain optimal k
+    val_opt = np.ceil(value(k_opt))  # obtain ceiling of optimal value given k
+    assert k_opt.is_integer()
+    assert val_opt.is_integer()
+    return int(k_opt), int(val_opt)
+
+
+def QR2(L1, L2, M):
+    """Find optimal cost for QROM over two values of L of size M.
+
+     Table[Ceiling[L1/2^k1]*Ceiling[L2/2^k2] + M*(2^(k1 + k2) - 1), {k1, 1,
+    10}, {k2, 1, 10}]
+    """
+    min_val = np.inf
+    for k1 in range(1, 11):
+        for k2 in range(1, 11):
+            test_val = np.ceil(L1 /
+                               (2**k1)) * np.ceil(L2 /
+                                                  (2**k2)) + M * (2**
+                                                                  (k1 + k2) - 1)
+            if test_val < min_val:
+                min_val = test_val
+    return int(min_val)
+
+
+def QI2(L1: int, Lv2: int):
+    """Find optimal value for inverse QROM over two L values.
+
+    QI2[L1_, L2_] :=
+    Min[Table[
+    Ceiling[L1/2^k1]*Ceiling[L2/2^k2] + 2^(k1 + k2), {k1, 1, 10}, {k2,
+      1, 10}]];
+    """
+    min_val = np.inf
+    for k1 in range(1, 11):
+        for k2 in range(1, 11):
+            test_val = np.ceil(L1 / (2**k1)) * np.ceil(Lv2 /
+                                                       (2**k2)) + 2**(k1 + k2)
+            if test_val < min_val:
+                min_val = test_val
+    return int(min_val)
diff --git a/src/openfermion/resource_estimates/pbc/sf/__init__.py b/src/openfermion/resource_estimates/pbc/sf/__init__.py
index 53f003f..3a6b416 100644
--- a/src/openfermion/resource_estimates/pbc/sf/__init__.py
+++ b/src/openfermion/resource_estimates/pbc/sf/__init__.py
@@ -9,3 +9,11 @@
 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from .compute_lambda_sf import compute_lambda
+    from .generate_costing_table_sf import generate_costing_table
+    from .sf_integrals import SingleFactorization
+    from .compute_sf_resources import compute_cost
diff --git a/src/openfermion/resource_estimates/pbc/sf/compute_sf_resources.py b/src/openfermion/resource_estimates/pbc/sf/compute_sf_resources.py
new file mode 100644
index 0000000..93c3029
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sf/compute_sf_resources.py
@@ -0,0 +1,286 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""
+Estimate the physical, logical, and Toffoli gate requirements for
+single-factorization qubitization.  Single-Factorization uses
+an LCU formed from a symmeterized Cholesky decomposition of the integrals
+"""
+from typing import Tuple
+
+import numpy as np
+from numpy.lib.scimath import arccos, arcsin
+from sympy import factorint
+
+from openfermion.resource_estimates.utils import QI
+from openfermion.resource_estimates.utils import QR2 as QR2_of
+from openfermion.resource_estimates.pbc.resources import (
+    ResourceEstimates,
+    QR3,
+    QR2,
+    QI2,
+)
+
+
+def compute_cost(
+        num_spin_orbs: int,
+        lambda_tot: float,
+        num_aux: int,
+        kmesh: list[int],
+        dE_for_qpe: float = 0.0016,
+        chi: int = 10,
+) -> ResourceEstimates:
+    """Determine fault-tolerant costs using single factorizated Hamiltonian.
+
+    Light wrapper around _compute_cost to automate choice of stps paramter.
+
+    Arguments:
+        num_spin_orbs: the number of spin-orbitals
+        lambda_tot: the lambda-value for the Hamiltonian
+        num_sym_unique: number of symmetry unique terms kept in the sparse
+            Hamiltonian
+        dE_for_qpe: allowable error in phase estimation
+        chi: the number of bits for the representation of the coefficients
+    Returns:
+        resources: sparse resources
+    """
+    # run once to determine stps parameter
+    init_cost = _compute_cost(
+        num_spin_orbs,
+        lambda_tot,
+        num_aux,
+        dE_for_qpe,
+        chi,
+        20_000,
+        kmesh[0],
+        kmesh[1],
+        kmesh[2],
+    )
+    steps = init_cost[0]
+    final_cost = _compute_cost(
+        num_spin_orbs,
+        lambda_tot,
+        num_aux,
+        dE_for_qpe,
+        chi,
+        steps,
+        kmesh[0],
+        kmesh[1],
+        kmesh[2],
+    )
+    estimates = ResourceEstimates(
+        toffolis_per_step=final_cost[0],
+        total_toffolis=final_cost[1],
+        logical_qubits=final_cost[2],
+    )
+    return estimates
+
+
+def _compute_cost(
+        n: int,
+        lam: float,
+        M: int,
+        dE: float,
+        chi: int,
+        stps: int,
+        Nkx: int,
+        Nky: int,
+        Nkz: int,
+        verbose: bool = False,
+) -> Tuple[int, int, int]:
+    """Determine fault-tolerant costs using SF decomposition in quantum chem
+
+    Args:
+        n: - the number of spin-orbitals. When using this for
+                  k-point sampling, this is equivalent to N times N_k.
+        lam: the lambda-value for the Hamiltonian
+        M: The combined number of values of n
+        dE:  allowable error in phase estimation
+        chi: equivalent to aleph_1 and aleph_2 in the document, the
+            number of bits for the representation of the coefficients
+        stps: an approximate number of steps to choose the precision of
+            single qubit rotations in preparation of the equal superposn state
+        Nkx: num k-points in x-direction
+        Nky: num k-points in y-direction
+        Nkz: num k-points in z-direction
+        verbose: do additional printing of intermediates?
+
+    Returns:
+        step_cost: Toffolis per step
+        total_cost: Total number of Toffolis
+        total_qubit_count: Total qubit count
+    """
+    nNk = (max(np.ceil(np.log2(Nkx)), 1) + max(np.ceil(np.log2(Nky)), 1) +
+           max(np.ceil(np.log2(Nkz)), 1))
+    Nk = Nkx * Nky * Nkz
+    L = Nk * n**2 // 2
+
+    # Number of trailing zeros by computing smallest prime factor exponent
+    # Number of trailing zeros
+    factors = factorint(L)
+    eta = factors[min(list(sorted(factors.keys())))]
+    if L % 2 == 1:
+        eta = 0
+
+    # Number of qubits for the first register
+    nL = np.ceil(np.log2(L))
+
+    # Number of qubits for p and q registers
+    nN = np.ceil(np.log2(n // 2))
+
+    nMN = np.ceil(np.log2(M * Nk + 1))
+
+    oh = [0] * 20
+    for p in range(20):
+        # JJG note: arccos arg may be > 1
+        # v = Round[2^p/(2*\[Pi])*ArcCos[2^nL/Sqrt[(L + 1)/2^\[Eta]]/2]];
+        v = np.round(
+            np.power(2, p + 1) / (2 * np.pi) *
+            arccos(np.power(2, nMN) / np.sqrt(M * Nk + 1) / 2))
+        #   oh[[p]] = stps*(1/N[Sin[3*ArcSin[Cos[v*2*\[Pi]/2^p]*Sqrt[(L +
+        #   1)/2^\[Eta]]/2^nL]]^2] - 1) + 4*p];
+        #   stps*(1/N[Sin[3*ArcSin[Cos[v*2*\[Pi]/2^p]*Sqrt[M*Nk + 1]/2^nMN]]^2]
+        #   - 1) + 4*p]
+        oh[p] = np.real(stps * (1 / (np.sin(3 * arcsin(
+            np.cos(v * 2 * np.pi / np.power(2, p + 1)) * np.sqrt(M * Nk + 1) /
+            np.power(2, nMN)))**2) - 1) + 4 * (p + 1))
+
+    # Bits of precision for rotation
+    br1 = int(np.argmin(oh) + 1)
+
+    oh2 = [0] * 20
+    for p in range(20):
+        # JJG note: arccos arg may be > 1
+        # v = Round[2^p/(2*\[Pi])*ArcCos[2^nL/Sqrt[(L + 1)/2^\[Eta]]/2]];
+        v = np.round(
+            np.power(2, p + 1) / (2 * np.pi) *
+            arccos(np.power(2, nL) / np.sqrt(L / np.power(2, eta)) / 2))
+        #   oh[[p]] = stps*(1/N[Sin[3*ArcSin[Cos[v*2*\[Pi]/2^p]*Sqrt[(L +
+        #   1)/2^\[Eta]]/2^nL]]^2] - 1) + 4*p]; oh2[[p]] =
+        #   stps*(1/N[Sin[3*ArcSin[Cos[v*2*\[Pi]/2^p]*Sqrt[L/2^\[Eta]]/2^nL]]^2]
+        #   - 1) + 4*p];
+        oh2[p] = np.real(stps * (1 / (np.sin(3 * arcsin(
+            np.cos(v * 2 * np.pi / np.power(2, p + 1)) *
+            np.sqrt(L / np.power(2, eta)) / np.power(2, nL)))**2) - 1) + 4 *
+                         (p + 1))
+    br2 = int(np.argmin(oh2) + 1)
+
+    # Cost of preparing the equal superposition state on the 1st register in 1a
+    cost1a = 2 * (3 * nMN + 2 * br1 - 9)
+
+    # We have added two qubits for ell <> 0.
+    # bMN = 2*(np.ceil(np.log2(Nk)) + np.ceil(np.log2(M))) + chi + 2
+    bMN = nMN + 2 * np.ceil(np.log2(Nk)) + chi + 2
+
+    # QROM costs for first register preparation in step 1(b)
+    cost1b = QR3(M * Nk + 1, bMN)[-1] + QI(M * Nk + 1)[-1]
+
+    # The inequality test of cost chi in step 1(c) and the controlled swap with
+    # cost nL + 1 in step 1(d) (and their inversions).
+    cost1cd = 2 * (chi + np.ceil(np.log2((Nk)) + np.ceil(np.log2(M)) + 1))
+
+    # Cost of preparing equal superposition of p and q registers in step 2(a).
+    cost2a = 4 * (3 * nL - 3 * eta + 2 * br2 - 9)
+
+    # The output size for the QROM for the second state preparation.
+    bp = 6 * np.ceil(np.log2(Nk) / 3) + 4 * nN + chi + 3
+
+    # Cost of QROMs for state preparation on p and q in step 2 (c).
+    cost2b = (QR2(M * Nk + 1, L, bp) + QI2(M * Nk + 1, L) + QR2(M * Nk, L, bp) +
+              QI2(M * Nk, L))
+
+    # The cost of the inequality test and controlled swap for the
+    # quantum alias sampling in steps 2(c) and (d).
+    cost2cd = 4 * (chi + 3 * np.ceil(np.log2(Nk) / 3) + 2 * nN + 1)
+
+    # Computing k - Q and k' - Q, then below we perform adjustments for non -
+    # modular arithmetic.
+    cost3 = 4 * nNk
+    if Nkx == 2**np.ceil(np.log2(Nkx)):
+        cost3 = cost3 - 2 * np.ceil(np.log2(Nkx))
+    if Nky == 2**np.ceil(np.log2(Nky)):
+        cost3 = cost3 - 2 * np.ceil(np.log2(Nky))
+    if Nkz == 2**np.ceil(np.log2(Nkz)):
+        cost3 = cost3 - 2 * np.ceil(np.log2(Nkz))
+
+    # The SELECT operation in step 4, which needs to be done twice.
+    cost4 = 2 * (3 * n * Nk - 4)
+
+    # The reflection used on the 2nd register in the middle, with cost in step 6
+    cost6 = nL + chi + 5
+
+    # Step 7 involves performing steps 2 to 5 again, which are accounted for
+    # above, but there is one extra Toffoli for checking that l <> 0.
+    cost7 = 4
+
+    # The total reflection for the quantum walk in step 9.
+    # with an extra + 1 for k - point sampling.
+    cost9 = nMN + nL + 2 * chi + 4
+
+    # The two Toffolis for the control for phase estimation and iteration,
+    # given in step 10.
+    cost10 = 2
+
+    # The number of steps.
+    iters = np.ceil(np.pi * lam / (dE * 2))
+
+    # The total Toffoli costs for a step.
+    cost = (cost1a + cost1b + cost1cd + cost2a + cost2b + cost2cd + cost3 +
+            cost4 + cost6 + cost7 + cost9 + cost10)
+
+    # Control for phase estimation and its iteration
+    ac1 = 2 * np.ceil(np.log2(iters)) - 1
+
+    # System qubits
+    ac2 = n * Nk
+
+    # First ell register that rotated and the flag for success
+    ac3 = nMN + 2
+
+    # State preparation on the first register
+    ac4 = bMN + chi + 1
+
+    # For p and q registers, the rotated qubit and the success flag
+    ac5 = nL + 2
+
+    kp = QR2_of(M * Nk + 1, L, bp)[:2]
+
+    ac6 = chi
+
+    # The equal superposition state for the second state preparation
+    ac7 = np.max([br1, br2])
+
+    # The phase gradient state
+    ac8 = 4
+
+    # The QROM on the p & q registers
+    ac9 = (kp[0] * kp[1] * bp + np.ceil(np.log2(
+        (M * Nk + 1) / kp[0])) + np.ceil(np.log2(L / kp[1])))
+
+    if verbose:
+        print("[*] Top of routine")
+        print("  [+] eta = ", eta)
+        print("  [+] nL = ", nL)
+        print("  [+] nN = ", nN)
+
+    total_qubit_count = ac1 + ac2 + ac3 + ac4 + ac5 + ac6 + ac7 + ac8 + ac9
+
+    # Sanity checks before returning as int
+    assert cost.is_integer()
+    assert iters.is_integer()
+    assert total_qubit_count.is_integer()
+
+    step_cost = int(cost)
+    total_cost = int(cost * iters)
+    total_qubit_count = int(total_qubit_count)
+
+    return step_cost, total_cost, total_qubit_count
diff --git a/src/openfermion/resource_estimates/pbc/sf/generate_costing_table_sf.py b/src/openfermion/resource_estimates/pbc/sf/generate_costing_table_sf.py
new file mode 100644
index 0000000..764c354
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sf/generate_costing_table_sf.py
@@ -0,0 +1,114 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import numpy.typing as npt
+import pandas as pd
+
+from pyscf.pbc import scf
+from pyscf.pbc.tools.k2gamma import kpts_to_kmesh
+
+from openfermion.resource_estimates.pbc.hamiltonian import (
+    build_hamiltonian,)
+from openfermion.resource_estimates.pbc.hamiltonian.cc_extensions import (
+    build_approximate_eris,
+    build_cc_inst,
+    build_approximate_eris_rohf,
+)
+from openfermion.resource_estimates.pbc.sf.compute_lambda_sf import (
+    compute_lambda,)
+from openfermion.resource_estimates.pbc.sf.compute_sf_resources import (
+    compute_cost,)
+from openfermion.resource_estimates.pbc.sf.sf_integrals import (
+    SingleFactorization,)
+from openfermion.resource_estimates.pbc.resources import PBCResources
+
+
+def generate_costing_table(
+        pyscf_mf: scf.HF,
+        naux_cutoffs: npt.NDArray[np.int32],
+        name: str = "pbc",
+        chi: int = 10,
+        dE_for_qpe=0.0016,
+        energy_method="MP2",
+) -> pd.DataFrame:
+    """Generate resource estimate costing table given a set of cutoffs for
+        single-factorized Hamiltonian.
+
+    Arguments:
+        pyscf_mf: k-point pyscf mean-field object
+        naux_cutoffs: Array of (integer) auxiliary index cutoff values
+        name: Optional descriptive name for simulation.
+        chi: the number of bits for the representation of the coefficients
+        dE_for_qpe: Phase estimation epsilon.
+
+    Returns
+        resources: Table of resource estimates.
+    """
+    kmesh = kpts_to_kmesh(pyscf_mf.cell, pyscf_mf.kpts)
+    cc_inst = build_cc_inst(pyscf_mf)
+    exact_eris = cc_inst.ao2mo()
+    if energy_method == "MP2":
+        energy_function = lambda x: cc_inst.init_amps(x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    elif energy_method == "CCSD":
+        energy_function = lambda x: cc_inst.kernel(eris=x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    else:
+        raise ValueError(f"Unknown value for energy_method: {energy_method}")
+
+    hcore, chol = build_hamiltonian(pyscf_mf)
+    num_spin_orbs = 2 * hcore[0].shape[-1]
+    num_kpts = np.prod(kmesh)
+
+    sf_resource_obj = PBCResources(
+        system_name=name,
+        num_spin_orbitals=num_spin_orbs,
+        num_kpts=num_kpts,
+        dE=dE_for_qpe,
+        chi=chi,
+        energy_method=energy_method,
+        exact_energy=np.real(reference_energy),
+    )
+    approx_eris = exact_eris
+    for cutoff in naux_cutoffs:
+        sf_helper = SingleFactorization(cholesky_factor=chol,
+                                        kmf=pyscf_mf,
+                                        naux=cutoff)
+        if pyscf_mf.cell.spin == 0:
+            approx_eris = build_approximate_eris(cc_inst,
+                                                 sf_helper,
+                                                 eris=approx_eris)
+        else:
+            approx_eris = build_approximate_eris_rohf(cc_inst,
+                                                      sf_helper,
+                                                      eris=approx_eris)
+        approx_energy, _, _ = energy_function(approx_eris)
+
+        sf_lambda = compute_lambda(hcore, sf_helper)
+        sparse_res_cost = compute_cost(
+            num_spin_orbs,
+            sf_lambda.lambda_total,
+            sf_lambda.num_aux,
+            list(kmesh),
+            chi=chi,
+            dE_for_qpe=dE_for_qpe,
+        )
+
+        sf_resource_obj.add_resources(
+            ham_properties=sf_lambda,
+            resource_estimates=sparse_res_cost,
+            cutoff=sf_lambda.num_aux,
+            approx_energy=np.real(approx_energy),
+        )
+
+    return sf_resource_obj.to_dataframe()
diff --git a/src/openfermion/resource_estimates/pbc/sparse/__init__.py b/src/openfermion/resource_estimates/pbc/sparse/__init__.py
index 53f003f..e232a78 100644
--- a/src/openfermion/resource_estimates/pbc/sparse/__init__.py
+++ b/src/openfermion/resource_estimates/pbc/sparse/__init__.py
@@ -9,3 +9,11 @@
 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from .compute_lambda_sparse import compute_lambda
+    from .generate_costing_table_sparse import generate_costing_table
+    from .sparse_integrals import SparseFactorization
+    from .compute_sparse_resources import compute_cost
diff --git a/src/openfermion/resource_estimates/pbc/sparse/compute_sparse_resources.py b/src/openfermion/resource_estimates/pbc/sparse/compute_sparse_resources.py
new file mode 100644
index 0000000..1883e96
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sparse/compute_sparse_resources.py
@@ -0,0 +1,189 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""
+Compute resource estimates for sparse LCU of k-point Hamiltonian
+"""
+from typing import Tuple
+import numpy as np
+from numpy.lib.scimath import arccos, arcsin  # has analytc continuation to cplx
+from sympy import factorint
+from openfermion.resource_estimates.utils import QI
+
+from openfermion.resource_estimates.pbc.resources import (ResourceEstimates,
+                                                          QR3)
+
+
+def compute_cost(
+        num_spin_orbs: int,
+        lambda_tot: float,
+        num_sym_unique: int,
+        kmesh: list[int],
+        dE_for_qpe: float = 0.0016,
+        chi: int = 10,
+) -> ResourceEstimates:
+    """Determine fault-tolerant costs using sparse representaion of Hamiltonian.
+
+    Light wrapper around _compute_cost to automate choice of stps paramter.
+
+    Arguments:
+        num_spin_orbs: the number of spin-orbitals
+        lambda_tot: the lambda-value for the Hamiltonian
+        num_sym_unique: number of symmetry unique terms kept in the sparse
+            Hamiltonian
+        dE_for_qpe: allowable error in phase estimation
+        chi: the number of bits for the representation of the coefficients
+    Returns:
+        resources: sparse resources
+    """
+    # run once to determine stps parameter
+    init_cost = _compute_cost(
+        num_spin_orbs,
+        lambda_tot,
+        num_sym_unique,
+        dE_for_qpe,
+        chi,
+        20_000,
+        *kmesh,
+    )
+    steps = init_cost[0]
+    final_cost = _compute_cost(
+        num_spin_orbs,
+        lambda_tot,
+        num_sym_unique,
+        dE_for_qpe,
+        chi,
+        steps,
+        *kmesh,
+    )
+    estimates = ResourceEstimates(
+        toffolis_per_step=final_cost[0],
+        total_toffolis=final_cost[1],
+        logical_qubits=final_cost[2],
+    )
+    return estimates
+
+
+def _compute_cost(
+        n: int,
+        lam: float,
+        d: int,
+        dE: float,
+        chi: int,
+        stps: int,
+        Nkx: int,
+        Nky: int,
+        Nkz: int,
+) -> Tuple[int, int, int]:
+    """Determine fault-tolerant costs using sparse representaion of Hamiltonian.
+
+    Arguments:
+        n: the number of spin-orbitals
+        lam: the lambda-value for the Hamiltonian
+        d: number of symmetry unique terms kept in the sparse Hamiltonian
+        dE: allowable error in phase estimation
+        chi: the number of bits for the representation of the coefficients
+        stps: an approximate number of steps to choose the precision
+            of single qubit rotations in preparation of the equal superposition
+            state
+        Nkx: number of k-points in x-direction
+        Nky: number of k-points in y-direction
+        Nkz: number of k-points in z-direction
+    Returns:
+        step_cost: Toffolis per step
+        total_cost: Total number of Toffolis
+        ancilla_cost: Total ancilla cost
+    """
+    if n % 2 != 0:
+        raise ValueError("The number of spin orbitals is always even!")
+
+    factors = factorint(d)
+    eta = factors[min(list(sorted(factors.keys())))]
+    if d % 2 == 1:
+        eta = 0
+
+    nN = np.ceil(np.log2(n // 2))
+    nNk = (max(np.ceil(np.log2(Nkx)), 1) + max(np.ceil(np.log2(Nky)), 1) +
+           max(np.ceil(np.log2(Nkz)), 1))
+    Nk = Nkx * Nky * Nkz
+
+    m = chi + 8 * nN + 6 * nNk + 5  # Eq 26
+
+    oh = [0] * 20
+
+    nM = (np.ceil(np.log2(d)) - eta) / 2
+
+    for p in range(2, 22):
+        # JJG note: arccos arg may be > 1
+        v = np.round(
+            np.power(2, p + 1) / (2 * np.pi) *
+            arccos(np.power(2, nM) / np.sqrt(d / 2**eta) / 2))
+        oh[p - 2] = np.real(stps * (1 / (np.sin(3 * arcsin(
+            np.cos(v * 2 * np.pi / np.power(2, p + 1)) * np.sqrt(d / 2**eta) /
+            np.power(2, nM)))**2) - 1) + 4 * (p + 1))
+
+    # Bits of precision for rotation
+    br = int(np.argmin(oh) + 1) + 2
+
+    # Hand selecting the k expansion factor
+    k1 = QR3(d, m)[0]
+
+    # Equation (A17)
+    cost = (QR3(d, m)[1] + QI(d)[1] + 6 * n * Nk + 8 * nN + 12 * nNk + 2 * chi +
+            7 * np.ceil(np.log2(d)) - 6 * eta + 4 * br - 8)
+
+    # The following are adjustments if we don't need to do explicit arithmetic
+    # to make subtraction modular
+    if Nkx == 2**np.ceil(np.log2(Nkx)):
+        cost = cost - 2 * np.ceil(np.log2(Nkx))
+    if Nky == 2**np.ceil(np.log2(Nky)):
+        cost = cost - 2 * np.ceil(np.log2(Nky))
+    if Nkz == 2**np.ceil(np.log2(Nkz)):
+        cost = cost - 2 * np.ceil(np.log2(Nkz))
+
+    # Number of iterations needed for the phase estimation.
+    iters = np.ceil(np.pi * lam / (dE * 2))
+
+    # The following are the number of qubits from the listing on page 39.
+
+    # Control registers for phase estimation and iteration on them.
+    ac1 = 2 * np.ceil(np.log2(iters)) - 1
+
+    # System qubits
+    ac2 = n * Nk
+
+    # The register used for the QROM
+    ac3 = np.ceil(np.log2(d))
+
+    # The phase gradient state
+    ac6 = br
+
+    # The equal superposition state for coherent alias sampling
+    ac7 = chi
+
+    # The ancillas used for QROM
+    ac8 = np.ceil(np.log2(d / np.power(2, k1))) + m * np.power(2, k1)
+
+    ac9 = 9
+
+    ancilla_cost = ac1 + ac2 + ac3 + ac6 + ac7 + ac8 + ac9
+
+    # Sanity checks before returning as int
+    assert cost.is_integer()
+    assert iters.is_integer()
+    assert ancilla_cost.is_integer()
+
+    step_cost = int(cost)
+    total_cost = int(cost * iters)
+    ancilla_cost = int(ancilla_cost)
+
+    return step_cost, total_cost, ancilla_cost
diff --git a/src/openfermion/resource_estimates/pbc/sparse/generate_costing_table_sparse.py b/src/openfermion/resource_estimates/pbc/sparse/generate_costing_table_sparse.py
new file mode 100644
index 0000000..9f9e8cb
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sparse/generate_costing_table_sparse.py
@@ -0,0 +1,115 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import pandas as pd
+
+import numpy as np
+
+from pyscf.pbc import scf
+from pyscf.pbc.tools.k2gamma import kpts_to_kmesh
+
+from openfermion.resource_estimates.pbc.resources import PBCResources
+from openfermion.resource_estimates.pbc.sparse.sparse_integrals import (
+    SparseFactorization,)
+from openfermion.resource_estimates.pbc.hamiltonian import (
+    build_hamiltonian,)
+from openfermion.resource_estimates.pbc.hamiltonian.cc_extensions import (
+    build_approximate_eris,
+    build_cc_inst,
+    build_approximate_eris_rohf,
+)
+from openfermion.resource_estimates.pbc.sparse.compute_lambda_sparse import (
+    compute_lambda,)
+from openfermion.resource_estimates.pbc.sparse.compute_sparse_resources import (
+    compute_cost,)
+
+
+def generate_costing_table(
+        pyscf_mf: scf.HF,
+        name="pbc",
+        chi: int = 10,
+        thresholds: np.ndarray = np.logspace(-1, -5, 6),
+        dE_for_qpe=0.0016,
+        energy_method="MP2",
+) -> pd.DataFrame:
+    """Generate resource estimate costing table given a set of cutoffs for
+        sparse Hamiltonian.
+
+    Arguments:
+        pyscf_mf: k-point pyscf mean-field object
+        name: Optional descriptive name for simulation.
+        chi: the number of bits for the representation of the coefficients
+        thresholds: Array of sparse thresholds to generate the table for.
+        dE_for_qpe: Phase estimation epsilon.
+
+    Returns
+        resources: Table of resource estimates.
+    """
+    kmesh = kpts_to_kmesh(pyscf_mf.cell, pyscf_mf.kpts)
+    cc_inst = build_cc_inst(pyscf_mf)
+    exact_eris = cc_inst.ao2mo()
+    if energy_method == "MP2":
+        energy_function = lambda x: cc_inst.init_amps(x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    elif energy_method == "CCSD":
+        energy_function = lambda x: cc_inst.kernel(eris=x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    else:
+        raise ValueError(f"Unknown value for energy_method: {energy_method}")
+
+    hcore, chol = build_hamiltonian(pyscf_mf)
+    num_spin_orbs = 2 * hcore[0].shape[-1]
+
+    num_kpts = np.prod(kmesh)
+
+    sparse_resource_obj = PBCResources(
+        system_name=name,
+        num_spin_orbitals=num_spin_orbs,
+        num_kpts=num_kpts,
+        dE=dE_for_qpe,
+        chi=chi,
+        energy_method=energy_method,
+        exact_energy=np.real(reference_energy),
+    )
+    approx_eris = exact_eris
+    for thresh in thresholds:
+        sparse_helper = SparseFactorization(cholesky_factor=chol,
+                                            kmf=pyscf_mf,
+                                            threshold=thresh)
+        if pyscf_mf.cell.spin == 0:
+            approx_eris = build_approximate_eris(cc_inst,
+                                                 sparse_helper,
+                                                 eris=approx_eris)
+        else:
+            approx_eris = build_approximate_eris_rohf(cc_inst,
+                                                      sparse_helper,
+                                                      eris=approx_eris)
+        approx_energy, _, _ = energy_function(approx_eris)
+
+        sparse_data = compute_lambda(hcore, sparse_helper)
+
+        sparse_res_cost = compute_cost(
+            num_spin_orbs,
+            sparse_data.lambda_total,
+            sparse_data.num_sym_unique,
+            list(kmesh),
+            dE_for_qpe=dE_for_qpe,
+            chi=chi,
+        )
+        sparse_resource_obj.add_resources(
+            ham_properties=sparse_data,
+            resource_estimates=sparse_res_cost,
+            cutoff=thresh,
+            approx_energy=np.real(approx_energy),
+        )
+
+    return sparse_resource_obj.to_dataframe()
diff --git a/src/openfermion/resource_estimates/pbc/thc/__init__.py b/src/openfermion/resource_estimates/pbc/thc/__init__.py
index 53f003f..f19c370 100644
--- a/src/openfermion/resource_estimates/pbc/thc/__init__.py
+++ b/src/openfermion/resource_estimates/pbc/thc/__init__.py
@@ -9,3 +9,11 @@
 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from .compute_lambda_thc import compute_lambda
+    from .generate_costing_table_thc import generate_costing_table
+    from .thc_integrals import KPTHCDoubleTranslation
+    from .compute_thc_resources import compute_cost
diff --git a/src/openfermion/resource_estimates/pbc/thc/compute_thc_resources.py b/src/openfermion/resource_estimates/pbc/thc/compute_thc_resources.py
new file mode 100644
index 0000000..134ee19
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/thc/compute_thc_resources.py
@@ -0,0 +1,308 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+# coverage:ignore
+""" Determine costs for THC decomposition in QC """
+from typing import Tuple, Union
+import numpy as np
+from numpy.lib.scimath import arccos, arcsin
+from sympy import factorint
+from openfermion.resource_estimates.utils import QI
+
+from openfermion.resource_estimates.pbc.resources.data_types import (
+    ResourceEstimates,)
+from openfermion.resource_estimates.pbc.resources.qrom import QR3
+
+
+def compute_cost(
+        num_spin_orbs: int,
+        lambda_tot: float,
+        thc_dim: int,
+        kmesh: list[int],
+        dE_for_qpe: float = 0.0016,
+        chi: int = 10,
+        beta: Union[int, None] = None,
+) -> ResourceEstimates:
+    """Determine fault-tolerant costs using THC factorization representation of
+        symmetry adapted integrals.
+
+    Light wrapper around _compute_cost.
+
+    Arguments:
+        num_spin_orbs: the number of spin-orbitals
+        lambda_tot: the lambda-value for the Hamiltonian
+        kmesh: kpoint mesh.
+        thc_dim: THC dimension (M).
+        dE_for_qpe: allowable error in phase estimation
+        chi: the number of bits for the representation of the coefficients
+        beta: the number of bits for controlled rotations
+    Returns:
+        resources: THC factorized resource estimates
+    """
+    # run once to determine stps parameter
+    thc_costs = _compute_cost(
+        n=num_spin_orbs,
+        lam=lambda_tot,
+        dE=dE_for_qpe,
+        chi=chi,
+        beta=beta,
+        M=thc_dim,
+        Nkx=kmesh[0],
+        Nky=kmesh[1],
+        Nkz=kmesh[2],
+        stps=20_000,  # not used
+    )
+    resources = ResourceEstimates(
+        toffolis_per_step=thc_costs[0],
+        total_toffolis=thc_costs[1],
+        logical_qubits=thc_costs[2],
+    )
+    return resources
+
+
+def _compute_cost(
+        n: int,
+        lam: float,
+        dE: float,
+        chi: int,
+        beta: int,
+        M: int,
+        Nkx: int,
+        Nky: int,
+        Nkz: int,
+        stps: int,
+        verbose: bool = False,
+) -> Tuple[int, int, int]:
+    """Determine fault-tolerant costs using THC decomposition in quantum chem
+
+    Arguments:
+        n: the number of spin-orbitals
+        lam: the lambda-value for the Hamiltonian
+        dE: allowable error in phase estimation
+        chi: equivalent to aleph in the document, the number of bits for
+            the representation of the coefficients
+        beta: equivalent to beth in the document, the number of bits for
+            the rotations
+        M: the dimension for the THC decomposition
+        Nkx: This is the number of values of k for the k - point sampling,
+                with each component
+        Nky: This is the number of values of k for the k - point sampling,
+                with each component
+        Nkz: This is the number of values of k for the k - point sampling,
+                with each component
+        stps: an approximate number of steps to choose the precision of
+            single qubit rotations in preparation of the equal superpositn state
+
+    Returns:
+        step_cost: Toffolis per step
+        total_cost: Total number of Toffolis
+        ancilla_cost: Total ancilla cost
+    """
+    nk = (max(np.ceil(np.log2(Nkx)), 1) + max(np.ceil(np.log2(Nky)), 1) +
+          max(np.ceil(np.log2(Nkz)), 1))
+    Nk = Nkx * Nky * Nkz
+
+    # (*Temporarily set as number of even numbers.*)
+    nc = 3 - Nkx % 2 - Nky % 2 - Nkz % 2
+
+    # The number of steps needed
+    iters = np.ceil(np.pi * lam / (2 * dE))
+
+    # This is the number of distinct items of data we need to output
+    # see Eq. (28).*)
+    d = int(32 * (Nk + 2**nc) * M**2 + n * Nk / 2)
+
+    # The number of bits used for the contiguous register
+    nc = np.ceil(np.log2(d))
+
+    nM = np.ceil(np.log2(M))
+
+    # The output size is 2* Log[M] for the alt values, chi for the keep value,
+    # and 2 for the two sign bits.
+    m = 2 * (2 * nM + nk + 8) + chi
+
+    oh = [0] * 20
+    for p in range(20):
+        # arccos arg may be > 1
+        v = np.round(
+            np.power(2, p + 1) / (2 * np.pi) *
+            arccos(np.power(2, nc) / np.sqrt(d) / 2))
+        oh[p] = stps * (1 / (np.sin(3 * arcsin(
+            np.cos(v * 2 * np.pi / np.power(2, p + 1)) * np.sqrt(d) /
+            np.power(2, nc)))**2) - 1) + 4 * (p + 1)
+
+    # Set it to be the number of bits that minimises the cost, usually 7.
+    # Python is 0-index, so need to add the one back in vs mathematica nb
+    br = np.argmin(oh) + 1
+
+    if d % 2 == 0:
+        factors = factorint(int(d))
+        eta = factors[min(list(sorted(factors.keys())))]
+    else:
+        eta = 0
+
+    cp1 = 2 * (3 * nc - 3 * eta + 2 * br - 9)
+
+    # This is the cost of the QROM for the state preparation in step 3 and its
+    cp3 = QR3(d, m)[1] + QI(d)[1]
+
+    # The cost for the inequality test in step 4 and its inverse.
+    cp4 = 2 * chi
+
+    # the cost of inequality test and controlled swap of mu and nu registers
+    cp5 = 2 * (2 * nM + nk + 7)
+
+    # The cost of an inequality test and controlled - swap of the mu and nu
+    # registers
+    cp6 = 4 * nM + 12
+
+    CPCP = cp1 + cp3 + cp4 + cp5 + cp6
+
+    # The cost of preparing the k superposition. The 7 here is the assumed
+    # number # of bits for the ancilla rotation which makes the probability of
+    # failure negligible.
+    cks = 4 * (Nkx + Nky + Nkz + 8 * nk + 6 * 7 - 24)
+
+    # The cost of the arithmetic computing k - Q, and controlling swaps for
+    # the one-body term
+    cka = 12 * nk
+
+    # This is the cost of swapping based on the spin register
+    cs1 = 3 * n * Nk / 2
+
+    # The cost of controlled swaps into working registers based on the k or
+    # k - Q value.
+    cs2 = 4 * n * (Nk - 1)
+
+    # The QROM for the rotation angles the first time.
+    cs2a = QR3(Nk * (M + n / 2), n * beta)[1] + QI(Nk * (M + n / 2))[1]
+
+    # The QROM for the rotation angles the second time.
+    cs2b = QR3(Nk * M, n * beta)[1] + QI(Nk * M)[1]
+
+    # The cost of the rotations.
+    cs3 = 16 * n * (beta - 2)
+
+    # Cost for constructing contiguous register for outputting rotations.
+    cs4 = (12 * Nk + 4 * np.ceil(np.log2(Nk * (M + n / 2))) +
+           4 * np.ceil(np.log2(Nk * M)))
+
+    # The cost of the controlled selection of the X vs Y.
+    cs5 = 2 + 4 * (Nk - 1)
+
+    # Cost of computing K-Q-G
+    cs6 = 2 * nk
+
+    # Cost of computing contiguous register.
+    cs6 = cs6 + np.ceil(np.log2(Nkx)) * np.ceil(np.log2(Nky))
+    cs6 = cs6 + np.ceil(np.log2(Nkx * Nky))
+    cs6 = cs6 + np.ceil(np.log2(Nkx * Nky)) * np.ceil(np.log2(Nky))
+    cs6 = cs6 + np.ceil(np.log2(Nk))
+    cs6 = cs6 + np.ceil(np.log2(Nk)) * np.ceil(np.log2(Nkx))
+    cs6 = cs6 + np.ceil(np.log2(Nkx * Nk))
+    cs6 = cs6 + np.ceil(np.log2(Nkx * Nk)) * np.ceil(np.log2(Nky))
+    cs6 = cs6 + np.ceil(np.log2(Nkx * Nky * Nk))
+    cs6 = cs6 + np.ceil(np.log2(Nkx * Nky * Nk)) * np.ceil(np.log2(Nkz))
+    cs6 = cs6 + np.ceil(np.log2(Nk**2))
+    cs6 = cs6 + np.ceil(np.log2(Nkx**2)) * np.ceil(np.log2(M))
+    cs6 = cs6 + np.ceil(np.log2(Nk**2 * M))
+
+    # the times for is for computing and uncomputing twice
+    cs6 = cs6 * 4
+
+    # The QROM cost once we computed the contiguous register
+    cs6 = cs6 + 2 * (QR3(Nk**2 * M, nk + chi)[1] + QI(Nk**2 * M)[1])
+
+    # The remaining state preparation cost with coherent alias sampling.
+    cs6 = cs6 + 4 * (nk + chi)
+
+    # The costs of generating the symmetry in Q, with a +1 for the one-body
+    # control.
+    cs6 = cs6 + 4 * nk + 1
+
+    # The total select cost.
+    CS = cks + cka + cs1 + cs2 + cs2a + cs2b + cs3 + cs4 + cs5 + cs6
+
+    # The reflection cost.
+    costref = nc + 2 * nk + 3 * chi + 9
+
+    cost = CPCP + CS + costref
+
+    # Qubits for control for phase estimation
+    ac1 = 2 * np.ceil(np.log2(iters + 1)) - 1
+
+    # system qubits
+    ac2 = n * Nk
+
+    # various control qubits
+    ac3 = nc + chi + nk + 10
+
+    # phase gradient state
+    ac4 = beta
+
+    # T state
+    ac5 = 1
+
+    # kp = 2^QRa[d, m]
+    kp = np.power(2, QR3(d, m)[0])
+
+    # First round of QROM.
+    ac12 = m * kp + np.ceil(np.log2(d / kp)) - 1
+
+    # Temporary qubits from QROM.
+    ac6 = m
+
+    # First round of QROM.
+    ac7 = m * (kp - 1) + np.ceil(np.log2(d / kp)) - 1
+
+    # Qubit from inequality test for state preparation.
+    ac8 = 1
+
+    # Qubit from inequality test for mu, nu
+    ac9 = 1
+
+    # QROM for constructing superposition on k.
+    ac10 = 2 * nk + 3 * 7
+
+    # The contiguous register and k-Q
+    ac11 = np.ceil(np.log2(Nk**2 * M) + nk)
+
+    # output for k-state
+    ac12 = nk + chi
+
+    kn = np.power(2, QR3(Nk**2 * M, nk + chi)[0])
+
+    ac13 = (kn - 1) * (nk + chi) + np.ceil(np.log2(Nk**2 * M)) - 1
+
+    ac14 = chi + 1
+
+    # The contiguous register.
+    ac15 = np.ceil(np.log2(Nk * (M + n / 2)))
+
+    kr = np.power(2, QR3(Nk * (M + n / 2), n * beta)[0])
+
+    #
+    ac16 = n * beta * kr + np.ceil(np.log2(Nk * (M + n / 2))) - 1
+
+    # common ancilla costs
+    aca = ac1 + ac2 + ac3 + ac4 + ac5 + ac6
+
+    acc = np.max([ac13, ac14 + ac15 + ac16])
+
+    acc = np.max([ac7, ac8 + ac9 + ac10 + ac11 + ac12 + acc])
+
+    step_cost = int(cost)
+    total_cost = int(cost * iters)
+    ancilla_cost = int(aca + acc)
+
+    # step-cost, Toffoli count, logical qubits
+    return step_cost, total_cost, ancilla_cost
diff --git a/src/openfermion/resource_estimates/pbc/thc/generate_costing_table_thc.py b/src/openfermion/resource_estimates/pbc/thc/generate_costing_table_thc.py
new file mode 100644
index 0000000..6420a46
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/thc/generate_costing_table_thc.py
@@ -0,0 +1,162 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+from dataclasses import dataclass
+from typing import Union
+
+import numpy as np
+import numpy.typing as npt
+import pandas as pd
+
+from pyscf.pbc import scf
+from pyscf.pbc.tools.k2gamma import kpts_to_kmesh
+
+from openfermion.resource_estimates.pbc.thc.factorizations.thc_jax import (
+    kpoint_thc_via_isdf,)
+from openfermion.resource_estimates.pbc.resources.data_types import PBCResources
+from openfermion.resource_estimates.pbc.thc.thc_integrals import (
+    KPTHCDoubleTranslation,)
+from openfermion.resource_estimates.pbc.hamiltonian.cc_extensions import (
+    build_approximate_eris,
+    build_cc_inst,
+    build_approximate_eris_rohf,
+)
+from openfermion.resource_estimates.pbc.hamiltonian import (
+    build_hamiltonian,)
+from openfermion.resource_estimates.pbc.thc.compute_lambda_thc import (
+    compute_lambda,)
+from openfermion.resource_estimates.pbc.thc.compute_thc_resources import (
+    compute_cost,)
+
+
+@dataclass
+class THCResources(PBCResources):
+    beta: int = 20
+
+
+def generate_costing_table(
+        pyscf_mf: scf.HF,
+        thc_rank_params: Union[list, npt.NDArray],
+        name="pbc",
+        chi: int = 10,
+        beta: int = 20,
+        dE_for_qpe: float = 0.0016,
+        reoptimize: bool = True,
+        bfgs_maxiter: int = 3000,
+        adagrad_maxiter: int = 3000,
+        fft_df_mesh: Union[None, list] = None,
+        energy_method: str = "MP2",
+) -> pd.DataFrame:
+    """Generate resource estimate costing table for THC Hamiltonian.
+
+    Arguments:
+        pyscf_mf: k-point pyscf mean-field object
+        thc_rank_params: Array of (integer) auxiliary index cutoff values
+        name: Optional descriptive name for simulation.
+        chi: the number of bits for the representation of the coefficients
+        beta: the number of bits for rotations.
+        dE_for_qpe: Phase estimation epsilon.
+        reoptimize: Whether or not to perform regularized reoptimization of THC
+            factors.
+        bfgs_maxiter: Max number of BFGS steps.
+        adagrad_maxiter: Max number of AdaGrad steps.
+        fft_df_mesh: FFTDF mesh for ISDF.
+        energy_method: Method to determine energy with (CCSD or MP2.)
+    Returns
+        resources: Table of resource estimates.
+    """
+    kmesh = kpts_to_kmesh(pyscf_mf.cell, pyscf_mf.kpts)
+    cc_inst = build_cc_inst(pyscf_mf)
+    exact_eris = cc_inst.ao2mo()
+    if energy_method.lower() == "mp2":
+        energy_function = lambda x: cc_inst.init_amps(x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    elif energy_method.lower() == "ccsd":
+        energy_function = lambda x: cc_inst.kernel(eris=x)
+        reference_energy, _, _ = energy_function(exact_eris)
+    else:
+        raise ValueError(f"Unknown value for energy_method: {energy_method}")
+
+    hcore, chol = build_hamiltonian(pyscf_mf)
+    num_spin_orbs = 2 * hcore[0].shape[-1]
+    num_kpts = np.prod(kmesh)
+    thc_resource_obj = THCResources(
+        system_name=name,
+        num_spin_orbitals=num_spin_orbs,
+        num_kpts=num_kpts,
+        dE=dE_for_qpe,
+        chi=chi,
+        beta=beta,
+        energy_method=energy_method,
+        exact_energy=np.real(reference_energy),
+    )
+
+    # For the ISDF guess we need an FFTDF MF object (really just need the grids
+    # so a bit of a hack) Since we're fitting to RSGDF it isn't very important
+    # what the value of the FFT mesh is, and there is some tradeoff (grid
+    # density vs comp time) which is not critical (see isdf notebook).
+    # Subsequent optimization attempts to fit to the RSGDF integrals which
+    # should hopefully be somewhat close to the FFTDF ones for ISDF to be a good
+    # starting point.
+    mf_fftdf = scf.KRHF(pyscf_mf.cell, pyscf_mf.kpts)
+    mf_fftdf.max_memory = 180000
+    mf_fftdf.kpts = pyscf_mf.kpts
+    mf_fftdf.e_tot = pyscf_mf.e_tot
+    mf_fftdf.mo_coeff = pyscf_mf.mo_coeff
+    mf_fftdf.mo_energy = pyscf_mf.mo_energy
+    mf_fftdf.mo_occ = pyscf_mf.mo_occ
+    if fft_df_mesh is not None:
+        mf_fftdf.with_df.mesh = fft_df_mesh
+    approx_eris = exact_eris
+    for thc_rank in thc_rank_params:
+        num_thc = thc_rank * num_spin_orbs // 2
+        kpt_thc, _ = kpoint_thc_via_isdf(
+            mf_fftdf,
+            chol,
+            num_thc,
+            perform_adagrad_opt=reoptimize,
+            perform_bfgs_opt=reoptimize,
+            bfgs_maxiter=bfgs_maxiter,
+            adagrad_maxiter=adagrad_maxiter,
+        )
+        thc_helper = KPTHCDoubleTranslation(kpt_thc.chi,
+                                            kpt_thc.zeta,
+                                            pyscf_mf,
+                                            chol=chol)
+        thc_lambda = compute_lambda(hcore, thc_helper)
+        kmesh = kpts_to_kmesh(pyscf_mf.cell, pyscf_mf.kpts)
+        if pyscf_mf.cell.spin == 0:
+            approx_eris = build_approximate_eris(cc_inst,
+                                                 thc_helper,
+                                                 eris=approx_eris)
+        else:
+            approx_eris = build_approximate_eris_rohf(cc_inst,
+                                                      thc_helper,
+                                                      eris=approx_eris)
+        approx_energy, _, _ = energy_function(approx_eris)
+        thc_res_cost = compute_cost(
+            num_spin_orbs,
+            thc_lambda.lambda_total,
+            thc_rank * num_spin_orbs // 2,
+            list(kmesh),
+            chi=chi,
+            beta=beta,
+            dE_for_qpe=dE_for_qpe,
+        )
+        thc_resource_obj.add_resources(
+            ham_properties=thc_lambda,
+            resource_estimates=thc_res_cost,
+            cutoff=thc_rank,
+            approx_energy=np.real(approx_energy),
+        )
+
+    return thc_resource_obj.to_dataframe()


Test patch:
diff --git a/src/openfermion/resource_estimates/pbc/df/compute_df_resources_test.py b/src/openfermion/resource_estimates/pbc/df/compute_df_resources_test.py
new file mode 100644
index 0000000..2672e25
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/df/compute_df_resources_test.py
@@ -0,0 +1,142 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.df.compute_df_resources import (
+        _compute_cost,
+        compute_cost,
+    )
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_costing():
+    nRe = 108
+    lamRe = 294.8
+    dE = 0.001
+    LRe = 360
+    LxiRe = 13031
+    chi = 10
+    betaRe = 16
+
+    # (*The Li et al orbitals.*)
+    nLi = 152
+    lamLi = 1171.2
+    LLi = 394
+    LxiLi = 20115
+    betaLi = 20
+
+    res = _compute_cost(nRe, lamRe, dE, LRe, LxiRe, chi, betaRe, 2, 2, 2,
+                        20_000)
+    res = _compute_cost(nRe, lamRe, dE, LRe, LxiRe, chi, betaRe, 2, 2, 2,
+                        res[0])
+    # 48250, 22343175750, 8174
+    assert np.isclose(res[0], 48250)
+    assert np.isclose(res[1], 22343175750)
+    assert np.isclose(res[2], 8174)
+
+    res = _compute_cost(nRe, lamRe, dE, LRe, LxiRe, chi, betaRe, 3, 5, 1,
+                        20_000)
+    res = _compute_cost(nRe, lamRe, dE, LRe, LxiRe, chi, betaRe, 3, 5, 1,
+                        res[0])
+    # 53146, 24610371366, 8945
+    assert np.isclose(res[0], 53146)
+    assert np.isclose(res[1], 24610371366)
+    assert np.isclose(res[2], 8945)
+
+    res = _compute_cost(nLi, lamLi, dE, LLi, LxiLi, chi, betaLi, 2, 2, 2,
+                        20_000)
+    res = _compute_cost(nLi, lamLi, dE, LLi, LxiLi, chi, betaLi, 2, 2, 2,
+                        res[0])
+    # print(res) # 79212, 145727663004, 13873
+    assert np.isclose(res[0], 79212)
+    assert np.isclose(res[1], 145727663004)
+    assert np.isclose(res[2], 13873)
+    res = _compute_cost(nLi, lamLi, dE, LLi, LxiLi, chi, betaLi, 3, 5, 1,
+                        res[0])
+    # print(res) # 86042, 158292930114, 14952
+    assert np.isclose(res[0], 86042)
+    assert np.isclose(res[1], 158292930114)
+    assert np.isclose(res[2], 14952)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_costing_helper():
+    nRe = 108
+    lamRe = 294.8
+    dE = 0.001
+    LRe = 360
+    LxiRe = 13031
+    chi = 10
+    betaRe = 16
+
+    # (*The Li et al orbitals.*)
+    nLi = 152
+    lamLi = 1171.2
+    LLi = 394
+    LxiLi = 20115
+    betaLi = 20
+
+    res = compute_cost(
+        num_spin_orbs=nRe,
+        lambda_tot=lamRe,
+        num_aux=LRe,
+        num_eig=LxiRe,
+        kmesh=[2, 2, 2],
+        dE_for_qpe=dE,
+        chi=chi,
+        beta=betaRe,
+    )
+    # 48250, 22343175750, 8174
+    assert np.isclose(res.toffolis_per_step, 48250)
+    assert np.isclose(res.total_toffolis, 22343175750)
+    assert np.isclose(res.logical_qubits, 8174)
+
+    res = compute_cost(
+        num_spin_orbs=nRe,
+        lambda_tot=lamRe,
+        num_aux=LRe,
+        num_eig=LxiRe,
+        kmesh=[3, 5, 1],
+        dE_for_qpe=dE,
+        chi=chi,
+        beta=betaRe,
+    )
+    # 53146, 24610371366, 8945
+    assert np.isclose(res.toffolis_per_step, 53146)
+    assert np.isclose(res.total_toffolis, 24610371366)
+    assert np.isclose(res.logical_qubits, 8945)
+
+    res = compute_cost(
+        num_spin_orbs=nLi,
+        lambda_tot=lamLi,
+        num_aux=LLi,
+        num_eig=LxiLi,
+        kmesh=[3, 5, 1],
+        dE_for_qpe=dE,
+        chi=chi,
+        beta=betaLi,
+    )
+    # print(res) # 79212, 145727663004, 13873
+    assert np.isclose(res.toffolis_per_step, 86042)
+    assert np.isclose(res.total_toffolis, 158292930114)
+    assert np.isclose(res.logical_qubits, 14952)
+
+
+if __name__ == "__main__":
+    test_costing()
diff --git a/src/openfermion/resource_estimates/pbc/df/generate_costing_table_df_test.py b/src/openfermion/resource_estimates/pbc/df/generate_costing_table_df_test.py
new file mode 100644
index 0000000..5c09a07
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/df/generate_costing_table_df_test.py
@@ -0,0 +1,42 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.df.\
+        generate_costing_table_df import generate_costing_table
+    from openfermion.resource_estimates.pbc.testing import (
+        make_diamond_113_szv,)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_generate_costing_table_df():
+    mf = make_diamond_113_szv()
+    thresh = np.array([0.1, 1e-2,
+                       1e-14])  # Eigenvalue threshold for second factorization.
+    table = generate_costing_table(mf,
+                                   cutoffs=thresh,
+                                   chi=10,
+                                   beta=22,
+                                   dE_for_qpe=1e-3)
+    assert np.allclose(table.dE, 1e-3)
+    assert np.allclose(table.chi, 10)
+    assert np.allclose(table.beta, 22)
+    assert np.allclose(table.cutoff, thresh)
+    assert np.allclose(table.num_aux, [648] * 3)
+    assert np.isclose(table.approx_energy.values[2],
+                      table.exact_energy.values[0])
diff --git a/src/openfermion/resource_estimates/pbc/resources/data_types_test.py b/src/openfermion/resource_estimates/pbc/resources/data_types_test.py
new file mode 100644
index 0000000..2522456
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/resources/data_types_test.py
@@ -0,0 +1,62 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.resources.data_types import (
+        PBCResources,
+        ResourceEstimates,
+    )
+    from openfermion.resource_estimates.pbc.hamiltonian import (
+        HamiltonianProperties,)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_pbc_resources():
+    nmo = 10
+    np.random.seed(7)
+    pbc_resources = PBCResources(
+        "pbc",
+        num_spin_orbitals=nmo,
+        num_kpts=12,
+        dE=1e-7,
+        chi=13,
+        exact_energy=-13.3,
+    )
+    for cutoff in np.logspace(-1, -3, 5):
+        lv = np.random.random(3)
+        lambdas = HamiltonianProperties(lambda_total=lv[0],
+                                        lambda_one_body=lv[1],
+                                        lambda_two_body=lv[2])
+        resource = ResourceEstimates(
+            toffolis_per_step=np.random.randint(0, 1000),
+            total_toffolis=np.random.randint(0, 1000),
+            logical_qubits=13,
+        )
+        pbc_resources.add_resources(
+            ham_properties=lambdas,
+            resource_estimates=resource,
+            approx_energy=-12,
+            cutoff=cutoff,
+        )
+    df = pbc_resources.to_dataframe()
+    assert np.allclose(
+        df.lambda_total.values,
+        [0.07630829, 0.97798951, 0.26843898, 0.38094113, 0.21338535],
+    )
+    assert (df.toffolis_per_step.values == [919, 366, 895, 787, 949]).all()
+    assert (df.total_toffolis.values == [615, 554, 391, 444, 112]).all()
\ No newline at end of file
diff --git a/src/openfermion/resource_estimates/pbc/resources/qrom_test.py b/src/openfermion/resource_estimates/pbc/resources/qrom_test.py
new file mode 100644
index 0000000..b7b6f71
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/resources/qrom_test.py
@@ -0,0 +1,41 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+
+from openfermion.resource_estimates.pbc.resources.qrom import (QR2, QI2)
+
+
+def test_qr2():
+    L = 728
+    npp = 182
+    bpp = 21
+    test_val = QR2(L + 1, npp, bpp)
+    assert np.isclose(test_val, 3416)
+
+    L = 56
+    npp = 28
+    bpp = 91
+    test_val = QR2(L + 1, npp, bpp)
+    assert np.isclose(test_val, 679)
+
+
+def test_qi2():
+    L1 = 728
+    npp = 182
+    test_val = QI2(L1 + 1, npp)
+    assert np.isclose(test_val, 785)
+
+    L1 = 56
+    npp = 28
+    test_val = QI2(L1 + 1, npp)
+    assert np.isclose(test_val, 88)
diff --git a/src/openfermion/resource_estimates/pbc/sf/compute_sf_resources_test.py b/src/openfermion/resource_estimates/pbc/sf/compute_sf_resources_test.py
new file mode 100644
index 0000000..bccc404
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sf/compute_sf_resources_test.py
@@ -0,0 +1,87 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.sf.compute_sf_resources import (
+        _compute_cost,
+        compute_cost,
+    )
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_estimate():
+    n = 152
+    lam = 3071.8
+    L = 275
+    dE = 0.001
+    chi = 10
+
+    res = _compute_cost(n, lam, L, dE, chi, 20_000, 3, 3, 3)
+    # 1663687, 8027577592851, 438447}
+    assert np.isclose(res[0], 1663687)
+    assert np.isclose(res[1], 8027577592851)
+    assert np.isclose(res[2], 438447)
+    res = _compute_cost(n, lam, L, dE, chi, res[0], 3, 3, 3)
+    assert np.isclose(res[0], 1663707)
+    assert np.isclose(res[1], 8027674096311)
+    assert np.isclose(res[2], 438452)
+
+    res = _compute_cost(n, lam, L, dE, chi, 20_000, 3, 5, 1)
+    # 907828, 4380427154244, 219526
+    assert np.isclose(res[0], 907828)
+    assert np.isclose(res[1], 4380427154244)
+    assert np.isclose(res[2], 219526)
+    res = _compute_cost(n, lam, L, dE, chi, res[0], 3, 5, 1)
+    assert np.isclose(res[0], 907828)
+    assert np.isclose(res[1], 4380427154244)
+    assert np.isclose(res[2], 219526)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_estimate_helper():
+    n = 152
+    lam = 3071.8
+    L = 275
+    dE = 0.001
+    chi = 10
+
+    res = compute_cost(
+        num_spin_orbs=n,
+        lambda_tot=lam,
+        num_aux=L,
+        kmesh=[3, 3, 3],
+        dE_for_qpe=dE,
+        chi=chi,
+    )
+    assert np.isclose(res.toffolis_per_step, 1663707)
+    assert np.isclose(res.total_toffolis, 8027674096311)
+    assert np.isclose(res.logical_qubits, 438452)
+
+    res = compute_cost(
+        num_spin_orbs=n,
+        lambda_tot=lam,
+        num_aux=L,
+        kmesh=[3, 5, 1],
+        dE_for_qpe=dE,
+        chi=chi,
+    )
+    # 1663687, 8027577592851, 438447}
+    assert np.isclose(res.toffolis_per_step, 907828)
+    assert np.isclose(res.total_toffolis, 4380427154244)
+    assert np.isclose(res.logical_qubits, 219526)
diff --git a/src/openfermion/resource_estimates/pbc/sf/generate_costing_table_sf_test.py b/src/openfermion/resource_estimates/pbc/sf/generate_costing_table_sf_test.py
new file mode 100644
index 0000000..c5c6a42
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sf/generate_costing_table_sf_test.py
@@ -0,0 +1,38 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+
+import pytest
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.sf.\
+        generate_costing_table_sf import generate_costing_table
+    from openfermion.resource_estimates.pbc.testing import (
+        make_diamond_113_szv,)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_generate_costing_table_sf():
+    mf = make_diamond_113_szv()
+    thresh = np.array([10, 54, 108])
+    table = generate_costing_table(mf,
+                                   naux_cutoffs=thresh,
+                                   chi=17,
+                                   dE_for_qpe=1e-3)
+    assert np.allclose(table.dE, 1e-3)
+    assert np.allclose(table.chi, 17)
+    assert np.allclose(table.num_aux, [10, 54, 108])
+    assert np.isclose(table.approx_energy.values[2],
+                      table.exact_energy.values[0])
diff --git a/src/openfermion/resource_estimates/pbc/sparse/compute_sparse_resources_test.py b/src/openfermion/resource_estimates/pbc/sparse/compute_sparse_resources_test.py
new file mode 100644
index 0000000..fa14e48
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sparse/compute_sparse_resources_test.py
@@ -0,0 +1,95 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.sparse.\
+        compute_sparse_resources import _compute_cost, compute_cost
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_compute_cost():
+    nRe = 108
+    lam_re = 2135.3
+    dRe = 705831
+    dE = 0.001
+    chi = 10
+
+    nLi = 152
+    lam_Li = 1547.3
+    dLi = 440501
+    dE = 0.001
+    chi = 10
+
+    Nkx = 2
+    Nky = 2
+    Nkz = 2
+
+    res = _compute_cost(nRe, lam_re, dRe, dE, chi, 20_000, Nkx, Nky, Nkz)
+    # res = _compute_cost(nRe, lam_re, dRe, dE, chi, res[0], Nkx, Nky, Nkz)
+    assert res[0] == 22962
+    assert res[1] == 77017349364
+    assert res[2] == 11335
+
+    res = _compute_cost(nRe, lam_re, dRe, dE, chi, 20_000, 3, 5, 1)
+    assert res[0] == 29004
+    assert res[1] == 97282954488
+    assert res[2] == 8060
+
+    res = _compute_cost(nLi, lam_Li, dLi, dE, chi, 20_000, Nkx, Nky, Nkz)
+    res = _compute_cost(nLi, lam_Li, dLi, dE, chi, res[0], Nkx, Nky, Nkz)
+    assert res[0] == 21426
+    assert res[1] == 52075764444
+    assert res[2] == 7015
+
+    res = _compute_cost(nLi, lam_Li, dLi, dE, chi, 20_000, 3, 5, 1)
+    assert res[0] == 28986
+    assert res[1] == 70450299084
+    assert res[2] == 9231
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_compute_cost_helper():
+    dE = 0.001
+    chi = 10
+
+    nLi = 152
+    lam_Li = 1547.3
+    dLi = 440501
+    dE = 0.001
+    chi = 10
+
+    Nkx = 2
+    Nky = 2
+    Nkz = 2
+
+    res = compute_cost(
+        num_spin_orbs=nLi,
+        lambda_tot=lam_Li,
+        num_sym_unique=dLi,
+        kmesh=[Nkx, Nky, Nkz],
+        dE_for_qpe=dE,
+        chi=chi,
+    )
+    assert res.toffolis_per_step == 21426
+    assert res.total_toffolis == 52075764444
+    assert res.logical_qubits == 7015
+
+
+if __name__ == "__main__":
+    test_compute_cost()
diff --git a/src/openfermion/resource_estimates/pbc/sparse/generate_costing_table_sparse_test.py b/src/openfermion/resource_estimates/pbc/sparse/generate_costing_table_sparse_test.py
new file mode 100644
index 0000000..56e04be
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/sparse/generate_costing_table_sparse_test.py
@@ -0,0 +1,40 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc import sparse
+    from openfermion.resource_estimates.pbc.testing import (
+        make_diamond_113_szv,)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_generate_costing_table_sparse():
+    mf = make_diamond_113_szv()
+    thresh = np.array([1e-1, 1e-2, 1e-12])
+    table = sparse.generate_costing_table(mf,
+                                          thresholds=thresh,
+                                          chi=17,
+                                          dE_for_qpe=1e-3)
+    assert np.allclose(table.dE, 1e-3)
+    assert np.allclose(table.chi, 17)
+    assert np.allclose(table.cutoff, thresh)
+    assert np.isclose(table.approx_energy.values[2],
+                      table.exact_energy.values[0])
+    assert not np.isclose(table.approx_energy.values[0],
+                          table.exact_energy.values[0])
diff --git a/src/openfermion/resource_estimates/pbc/thc/compute_thc_resources_test.py b/src/openfermion/resource_estimates/pbc/thc/compute_thc_resources_test.py
new file mode 100644
index 0000000..ec0b82c
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/thc/compute_thc_resources_test.py
@@ -0,0 +1,104 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.thc.compute_thc_resources import (
+        _compute_cost,
+        compute_cost,
+    )
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_thc_resources():
+    lam = 307.68
+    dE = 0.001
+    n = 108
+    chi = 10
+    beta = 16
+    M = 350
+
+    res = _compute_cost(n, lam, dE, chi, beta, M, 1, 1, 1, 20_000)
+    # print(res) # 26205, 12664955115, 2069
+    print(res)  # (80098, 38711603694, 17630)
+    assert np.isclose(res[0], 80098)
+    assert np.isclose(res[1], 38711603694)
+    assert np.isclose(res[2], 17630)
+
+    res = _compute_cost(n, lam, dE, chi, beta, M, 3, 3, 3, 20_000)
+    # print(res)  # {205788, 99457957764, 78813
+    print(res)  # (270394, 130682231382, 78815)
+    assert np.isclose(res[0], 270394)
+    assert np.isclose(res[1], 130682231382)
+    assert np.isclose(res[2], 78815)
+
+    res = _compute_cost(n, lam, dE, chi, beta, M, 3, 5, 1, 20_000)
+    # print(res)  # 151622, 73279367466, 39628
+    print(res)  #  (202209, 97728216327, 77517)
+    assert np.isclose(res[0], 202209)
+    assert np.isclose(res[1], 97728216327)
+    assert np.isclose(res[2], 77517)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+def test_thc_resources_helper():
+    lam = 307.68
+    dE = 0.001
+    n = 108
+    chi = 10
+    beta = 16
+    M = 350
+
+    res = compute_cost(
+        num_spin_orbs=n,
+        lambda_tot=lam,
+        thc_dim=M,
+        kmesh=[1, 1, 1],
+        dE_for_qpe=dE,
+        chi=chi,
+        beta=beta,
+    )
+    assert np.isclose(res.toffolis_per_step, 80098)
+    assert np.isclose(res.total_toffolis, 38711603694)
+    assert np.isclose(res.logical_qubits, 17630)
+
+    res = compute_cost(
+        num_spin_orbs=n,
+        lambda_tot=lam,
+        thc_dim=M,
+        kmesh=[3, 3, 3],
+        dE_for_qpe=dE,
+        chi=chi,
+        beta=beta,
+    )
+    assert np.isclose(res.toffolis_per_step, 270394)
+    assert np.isclose(res.total_toffolis, 130682231382)
+    assert np.isclose(res.logical_qubits, 78815)
+
+    res = compute_cost(
+        num_spin_orbs=n,
+        lambda_tot=lam,
+        thc_dim=M,
+        kmesh=[3, 5, 1],
+        dE_for_qpe=dE,
+        chi=chi,
+        beta=beta,
+    )
+    assert np.isclose(res.toffolis_per_step, 202209)
+    assert np.isclose(res.total_toffolis, 97728216327)
+    assert np.isclose(res.logical_qubits, 77517)
diff --git a/src/openfermion/resource_estimates/pbc/thc/generate_costing_table_thc_test.py b/src/openfermion/resource_estimates/pbc/thc/generate_costing_table_thc_test.py
new file mode 100644
index 0000000..ab5ad7c
--- /dev/null
+++ b/src/openfermion/resource_estimates/pbc/thc/generate_costing_table_thc_test.py
@@ -0,0 +1,44 @@
+# coverage: ignore
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.pbc.thc.\
+        generate_costing_table_thc import generate_costing_table
+    from openfermion.resource_estimates.pbc.testing import (
+        make_diamond_113_szv,)
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES,
+                    reason='pyscf and/or jax not installed.')
+@pytest.mark.slow
+def test_generate_costing_table_thc():
+    mf = make_diamond_113_szv()
+    thc_rank_params = np.array([2, 4, 6])
+    table = generate_costing_table(
+        mf,
+        thc_rank_params=thc_rank_params,
+        chi=10,
+        beta=22,
+        dE_for_qpe=1e-3,
+        bfgs_maxiter=10,
+        adagrad_maxiter=10,
+        fft_df_mesh=[11] * 3,
+    )
+    assert np.allclose(table.dE, 1e-3)
+    assert np.allclose(table.chi, 10)
+    assert np.allclose(table.beta, 22)
+    assert np.allclose(table.cutoff, thc_rank_params)