diff --git a/src/density_fitting.cpp b/src/density_fitting.cpp
index 5d128370..6a42d6fd 100644
--- a/src/density_fitting.cpp
+++ b/src/density_fitting.cpp
@@ -783,6 +783,15 @@ arma::mat DensityFit::calcJ(const arma::mat & P) const {
 
   // Get the expansion coefficients
   arma::vec c=compute_expansion(P);
+  return digestJ(c);
+}
+
+arma::mat DensityFit::calcJ(const arma::vec & c) const {
+  if(c.n_elem != Naux) {
+    std::ostringstream oss;
+    oss << "Error in DensityFit: Naux = " << Naux << ", c.n_elem = " << c.n_elem << "!\n";
+    throw std::logic_error(oss.str());
+  }
 
   arma::mat J(Nbf,Nbf);
   J.zeros();
diff --git a/src/density_fitting.h b/src/density_fitting.h
index 3445a06a..5d6337d1 100644
--- a/src/density_fitting.h
+++ b/src/density_fitting.h
@@ -149,6 +149,8 @@ class DensityFit {
   arma::mat calcJ(const arma::mat & P) const;
   /// Get Coulomb matrix from P
   std::vector<arma::mat> calcJ(const std::vector<arma::mat> & P) const;
+  /// Digest J matrix from computed expansion
+  arma::mat digest(const arma::vec & gamma) const;
 
   /// Calculate force from P
   arma::vec forceJ(const arma::mat & P);