diff --git a/src/openfermion/resource_estimates/df/__init__.py b/src/openfermion/resource_estimates/df/__init__.py
index 2e25d980e..7f31c73a2 100644
--- a/src/openfermion/resource_estimates/df/__init__.py
+++ b/src/openfermion/resource_estimates/df/__init__.py
@@ -18,5 +18,6 @@
 if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
     from .compute_cost_df import compute_cost
     from .compute_lambda_df import compute_lambda
+    from .compute_lambda_df_with_bliss import compute_lambda_df_with_bliss_post_processing
     from .factorize_df import factorize
     from .generate_costing_table_df import generate_costing_table
diff --git a/src/openfermion/resource_estimates/df/compute_lambda_df_with_bliss.py b/src/openfermion/resource_estimates/df/compute_lambda_df_with_bliss.py
new file mode 100644
index 000000000..d0d0ac90a
--- /dev/null
+++ b/src/openfermion/resource_estimates/df/compute_lambda_df_with_bliss.py
@@ -0,0 +1,40 @@
+# coverage:ignore
+""" Compute lambda for double low rank factoriz. method of von Burg, et al, with BLISS post-processing of DF fragments"""
+import numpy as np
+from openfermion.resource_estimates.molecule import pyscf_to_cas
+
+
+def compute_lambda_df_with_bliss_post_processing(pyscf_mf, df_factors):
+    """Compute lambda for Hamiltonian using DF method of von Burg, et al. where DF is post-processed with a low-rank-preserving BLISS operator
+
+    Args:
+        pyscf_mf - Pyscf mean field object
+        df_factors (ndarray) - (N x N x rank) array of DF factors from ERI
+
+    Returns:
+        lambda_tot (float) - lambda value for the double factorized Hamiltonian
+    """
+    # grab tensors from pyscf_mf object
+    h1, eri_full, _, num_alpha, num_beta = pyscf_to_cas(pyscf_mf)
+    num_elec = num_alpha + num_beta
+    num_orb = h1.shape[0]
+
+    # two body contributions
+    lambda_F = 0.0
+    h1_correction_BLISS = np.zeros([num_orb, num_orb])
+    for vector in range(df_factors.shape[2]):
+        Lij = df_factors[:, :, vector]
+        e = np.linalg.eigvalsh(Lij)
+        s = np.median(e)
+        lambda_F += 0.25 * np.sum(np.abs(e - s)) ** 2
+        h1_correction_BLISS += s * (num_elec - num_orb) * Lij
+
+    # one body contributions
+    T = h1 - 0.5 * np.einsum("illj->ij", eri_full) + np.einsum("llij->ij", eri_full) + h1_correction_BLISS
+    e, _ = np.linalg.eigh(T)
+    s = np.median(e)
+    lambda_T = np.sum(np.abs(e - s))
+
+    lambda_tot = lambda_T + lambda_F
+
+    return lambda_tot
\ No newline at end of file
diff --git a/src/openfermion/resource_estimates/df/compute_lambda_df_with_bliss_test.py b/src/openfermion/resource_estimates/df/compute_lambda_df_with_bliss_test.py
new file mode 100644
index 000000000..9e11c4e60
--- /dev/null
+++ b/src/openfermion/resource_estimates/df/compute_lambda_df_with_bliss_test.py
@@ -0,0 +1,25 @@
+# coverage:ignore
+"""Test cases for compute_lambda_df.py
+"""
+from os import path
+
+import numpy as np
+import pytest
+
+from openfermion.resource_estimates import HAVE_DEPS_FOR_RESOURCE_ESTIMATES, df
+
+if HAVE_DEPS_FOR_RESOURCE_ESTIMATES:
+    from openfermion.resource_estimates.molecule import load_casfile_to_pyscf
+
+
+@pytest.mark.skipif(not HAVE_DEPS_FOR_RESOURCE_ESTIMATES, reason="pyscf and/or jax not installed.")
+def test_reiher_df_bliss_lambda():
+    """Reproduce Reiher et al orbital DF lambda for DF+LRPS method of J. Chem. Theory Comput. 2025, 21, 2, 703–713"""
+
+    THRESH = 1e-8
+    NAME = path.join(path.dirname(__file__), '../integrals/eri_reiher.h5')
+    _, mf = load_casfile_to_pyscf(NAME, num_alpha=27, num_beta=27)
+    eri_rr, LR, L, Lxi = df.factorize(mf._eri, thresh=THRESH)
+    total_lambda = df.compute_lambda_df_with_bliss_post_processing(mf, LR)
+    assert eri_rr.shape[0] * 2 == 108
+    assert np.isclose(np.round(total_lambda, decimals=1), 169.4)
\ No newline at end of file