diff --git a/exachem/scf/scf_atom.cpp b/exachem/scf/scf_atom.cpp
index 954dfa2..7f14d7a 100644
--- a/exachem/scf/scf_atom.cpp
+++ b/exachem/scf/scf_atom.cpp
@@ -299,11 +299,10 @@ void compute_sad_guess(ExecutionContext& ec, ScalapackInfo& scalapack_info, Syst
     // std::tie(X_atom, Xinv, obs_rank, S_condition_number, XtX_condition_number, n_illcond) =
     //     gensqrtinv(ec, S_atom_eig, false, S_condition_number_threshold);
 
-    Tensor<TensorType> X_alpha, X_beta;
-    std::tie(obs_rank, S_condition_number, XtX_condition_number) =
-      gensqrtinv_atscf(ec, sys_data, scf_vars, scalapack_info, S_atom, X_alpha, X_beta, tAO_atom,
-                       false, S_condition_number_threshold);
-    Matrix X_atom = tamm_to_eigen_matrix(X_alpha);
+    Matrix X_atom;
+    std::tie(X_atom, obs_rank, S_condition_number, XtX_condition_number) =
+      gensqrtinv_atscf(ec, sys_data, scf_vars, scalapack_info, S_atom, tAO_atom, false,
+                       S_condition_number_threshold);
 
     // if(rank == 0) cout << std::setprecision(6) << "X_atom: " << endl << X_atom << endl;
 
diff --git a/exachem/scf/scf_common.cpp b/exachem/scf/scf_common.cpp
index e0f2380..e9a7d2b 100644
--- a/exachem/scf/scf_common.cpp
+++ b/exachem/scf/scf_common.cpp
@@ -290,14 +290,6 @@ compute_AO_tiles(const ExecutionContext& ec, const SystemData& sys_data, libint2
     shell_tile_map.push_back(shells.size() - 1);
   }
 
-  // std::vector<int> vtc(AO_tiles.size());
-  // std::iota (std::begin(vtc), std::end(vtc), 0);
-  // cout << "AO tile indexes = " << vtc;
-  // cout << "orig AO tiles = " << AO_tiles;
-
-  // cout << "print new opt AO tiles = " << AO_opttiles;
-  // cout << "print shell-tile map = " << shell_tile_map;
-
   return std::make_tuple(shell_tile_map, AO_tiles, AO_opttiles);
 }
 
@@ -312,10 +304,7 @@ void compute_orthogonalizer(ExecutionContext& ec, SystemData& sys_data, SCFVars&
   auto rank  = ec.pg().rank();
 
   // compute orthogonalizer X such that X.transpose() . S . X = I
-  double XtX_condition_number; // condition number of "re-conditioned"
-                               // overlap obtained as Xinv.transpose() . Xinv
-  // by default assume can manage to compute with condition number of S <= 1/eps
-  // this is probably too optimistic, but in well-behaved cases even 10^11 is OK
+  double XtX_condition_number;
   size_t obs_rank;
   double S_condition_number;
   double S_condition_number_threshold = sys_data.options_map.scf_options.tol_lindep;
@@ -366,8 +355,6 @@ void compute_hamiltonian(ExecutionContext& ec, const SCFVars& scf_vars,
     (ttensors.H1(mu, nu)  =  ttensors.T1(mu, nu))
     (ttensors.H1(mu, nu) +=  ttensors.V1(mu, nu)).execute();
   // clang-format on
-
-  // tamm::scale_ip(ttensors.H1(),2.0);
 }
 
 template<typename TensorType>
@@ -425,19 +412,9 @@ void compute_density(ExecutionContext& ec, const SystemData& sys_data, const SCF
   sch.execute();
 
   // compute D in eigen for subsequent fock build
-  {
-    // Matrix& C_occ_a = etensors.G_alpha;
-    // Matrix& C_occ_b = etensors.G_beta;
-    // C_occ_a.resize(sys_data.nbf_orig, sys_data.nelectrons_alpha);
-    // if(is_uhf) C_occ_b.resize(sys_data.nbf_orig, sys_data.nelectrons_beta);
-
-    if(!scf_vars.do_dens_fit) {
-      tamm_to_eigen_tensor(ttensors.D_alpha, D_alpha);
-      if(is_uhf) tamm_to_eigen_tensor(ttensors.D_beta, etensors.D_beta);
-    }
-
-    // D_alpha = dfac * C_occ_a * C_occ_a.transpose();
-    // if(is_uhf) etensors.D_beta = C_occ_b * C_occ_b.transpose();
+  if(!scf_vars.do_dens_fit) {
+    tamm_to_eigen_tensor(ttensors.D_alpha, D_alpha);
+    if(is_uhf) tamm_to_eigen_tensor(ttensors.D_beta, etensors.D_beta);
   }
 
   ec.pg().barrier();
@@ -665,8 +642,9 @@ std::tuple<size_t, double, double> gensqrtinv(ExecutionContext& ec, SystemData&
     scalapackpp::BlockCyclicDist2D* blockcyclic_dist = scalapack_info.blockcyclic_dist.get();
     const tamm::Tile                mb               = blockcyclic_dist->mb();
 
-    TiledIndexSpace tN_bc{IndexSpace{range(sys_data.nbf_orig)}, mb};
-    Tensor<T>       S_BC{tN_bc, tN_bc};
+    scf_vars.tN_bc         = TiledIndexSpace{IndexSpace{range(sys_data.nbf_orig)}, mb};
+    TiledIndexSpace& tN_bc = scf_vars.tN_bc;
+    Tensor<T>        S_BC{tN_bc, tN_bc};
     V_sca = {tN_bc, tN_bc};
     S_BC.set_block_cyclic({scalapack_info.npr, scalapack_info.npc});
     V_sca.set_block_cyclic({scalapack_info.npr, scalapack_info.npc});
@@ -682,15 +660,10 @@ std::tuple<size_t, double, double> gensqrtinv(ExecutionContext& ec, SystemData&
     if(grid.ipr() >= 0 and grid.ipc() >= 0) {
       auto desc_S = desc_lambda(N, N);
       auto desc_V = desc_lambda(N, N);
-      // scalapackpp::BlockCyclicMatrix<T> V_sca(grid, N, N, mb, mb);
 
       /*info=*/scalapackpp::hereig(scalapackpp::Job::Vec, scalapackpp::Uplo::Lower, desc_S[2],
                                    S_BC.access_local_buf(), 1, 1, desc_S, eps.data(),
                                    V_sca.access_local_buf(), 1, 1, desc_V);
-
-      // Gather results
-      // if( scalapack_info.pg.rank() == 0 ) V.resize(N,N);
-      // V_sca.gather_from( N, N, V.data(), N, 0, 0 );
     }
 
     Tensor<T>::deallocate(S_BC);
@@ -735,20 +708,45 @@ std::tuple<size_t, double, double> gensqrtinv(ExecutionContext& ec, SystemData&
     }
   }
 
-  if(world_size > 1) {
-    // TODO: Should buffer this
-    // ec.pg().broadcast(&n_cond, 0);
-    ec.pg().broadcast(&n_illcond, 0);
-    // ec.pg().broadcast( &condition_number,        0 );
-    // ec.pg().broadcast(&result_condition_number, 0);
-  }
+  if(world_size > 1) ec.pg().broadcast(&n_illcond, 0);
   n_cond = N - n_illcond;
 
+  sys_data.n_lindep = n_illcond;
+  sys_data.nbf      = n_cond;
+
+  scf_vars.tAO_ortho = TiledIndexSpace{IndexSpace{range((size_t) sys_data.nbf)},
+                                       sys_data.options_map.scf_options.AO_tilesize};
+
+  Tensor<T> X_tmp{scf_vars.tAO, scf_vars.tAO_ortho};
+  Tensor<T> eps_tamm{scf_vars.tAO_ortho};
+  Tensor<T>::allocate(&ec, X_tmp, eps_tamm);
+
+  if(world_rank == 0) {
+    std::vector<T> epso(first_above_thresh, eps.end());
+    std::transform(epso.begin(), epso.end(), epso.begin(),
+                   [](auto& c) { return 1.0 / std::sqrt(c); });
+    tamm::vector_to_tamm_tensor(eps_tamm, epso);
+  }
+  ec.pg().barrier();
+
+#if defined(USE_SCALAPACK)
+  if(scalapack_info.comm != MPI_COMM_NULL) {
+    const tamm::Tile _mb = (scalapack_info.blockcyclic_dist.get())->mb();
+    scf_vars.tNortho_bc  = TiledIndexSpace{IndexSpace{range(sys_data.nbf)}, _mb};
+    ttensors.X_alpha     = {scf_vars.tN_bc, scf_vars.tNortho_bc};
+    ttensors.X_alpha.set_block_cyclic({scalapack_info.npr, scalapack_info.npc});
+    Tensor<TensorType>::allocate(&scalapack_info.ec, ttensors.X_alpha);
+  }
+#else
+  ttensors.X_alpha = {scf_vars.tAO, scf_vars.tAO_ortho};
+  sch.allocate(ttensors.X_alpha).execute();
+#endif
+
 #if defined(USE_SCALAPACK)
   if(scalapack_info.comm != MPI_COMM_NULL) {
     Tensor<T> V_t = from_block_cyclic_tensor(V_sca);
     Tensor<T> X_t = tensor_block(V_t, {n_illcond, 0}, {N, N}, {1, 0});
-    if(scalapack_info.pg.rank() == 0) X = tamm_to_eigen_matrix<T>(X_t);
+    tamm::from_dense_tensor(X_t, X_tmp);
     Tensor<T>::deallocate(V_sca, V_t, X_t);
   }
 #else
@@ -756,63 +754,34 @@ std::tuple<size_t, double, double> gensqrtinv(ExecutionContext& ec, SystemData&
     // auto* V_cond = Vbuf + n_illcond * N;
     Matrix V_cond = V.block(n_illcond, 0, N - n_illcond, N);
     V.resize(0, 0);
-    X.resize(N, n_cond); // Xinv.resize( N, n_cond );
-    // Matrix V_cond(n_cond,N);
-    // V_cond = Eigen::Map<Matrix>(Vbuf + n_illcond * N,n_cond,N);
+    X.resize(N, n_cond);
     X = V_cond.transpose();
     V_cond.resize(0, 0);
+    eigen_to_tamm_tensor(X_tmp, X);
+    X.resize(0, 0);
   }
+  ec.pg().barrier();
 #endif
 
-  if(world_rank == 0) {
-    // Form canonical X/Xinv
-    for(auto i = 0; i < n_cond; ++i) {
-      const double srt = std::sqrt(*(first_above_thresh + i));
-
-      // X is row major...
-      auto* X_col = X.data() + i;
-      // auto* Xinv_col = Xinv.data() + i;
-
-      blas::scal(N, 1. / srt, X_col, n_cond);
-      // blas::scal( N, srt, Xinv_col, n_cond );
-    }
-
-    if(symmetric) {
-      assert(not symmetric);
-      /*
-      // X is row major, thus we need to form X**T = V_cond * X**T
-      Matrix TMP = X;
-      X.resize( N, N );
-      blas::gemm( blas::Op::NoTrans, blas::Op::NoTrans, N, N, n_cond, 1., V_cond, N, TMP.data(),
-      n_cond, 0., X.data(), N );
-      */
-    }
-  } // compute on root
+  Tensor<T> X_comp{scf_vars.tAO, scf_vars.tAO_ortho};
+  auto      mu   = scf_vars.tAO.label("all");
+  auto      mu_o = scf_vars.tAO_ortho.label("all");
 
-  sys_data.n_lindep = n_illcond;
-  sys_data.nbf      = sys_data.nbf_orig - sys_data.n_lindep;
+#if defined(USE_SCALAPACK)
+  sch.allocate(X_comp).execute();
+#else
+  X_comp = ttensors.X_alpha;
+#endif
 
-  scf_vars.tAO_ortho = TiledIndexSpace{IndexSpace{range(0, (size_t) (sys_data.nbf))},
-                                       sys_data.options_map.scf_options.AO_tilesize};
+  sch(X_comp(mu, mu_o) = X_tmp(mu, mu_o) * eps_tamm(mu_o)).deallocate(X_tmp, eps_tamm).execute();
 
 #if defined(USE_SCALAPACK)
   if(scalapack_info.comm != MPI_COMM_NULL) {
-    const tamm::Tile _mb = (scalapack_info.blockcyclic_dist.get())->mb();
-    scf_vars.tN_bc       = TiledIndexSpace{IndexSpace{range(sys_data.nbf_orig)}, _mb};
-    scf_vars.tNortho_bc  = TiledIndexSpace{IndexSpace{range(sys_data.nbf)}, _mb};
-    ttensors.X_alpha     = {scf_vars.tN_bc, scf_vars.tNortho_bc};
-    ttensors.X_alpha.set_block_cyclic({scalapack_info.npr, scalapack_info.npc});
-    Tensor<TensorType>::allocate(&scalapack_info.ec, ttensors.X_alpha);
+    tamm::to_block_cyclic_tensor(X_comp, ttensors.X_alpha);
   }
-#else
-  ttensors.X_alpha = {scf_vars.tAO, scf_vars.tAO_ortho};
-  sch.allocate(ttensors.X_alpha).execute();
+  sch.deallocate(X_comp).execute();
 #endif
 
-  if(world_rank == 0) eigen_to_tamm_tensor(ttensors.X_alpha, X);
-
-  ec.pg().barrier();
-
   return std::make_tuple(size_t(n_cond), condition_number, result_condition_number);
 }
 
@@ -972,13 +941,6 @@ Matrix compute_schwarz_ints(ExecutionContext& ec, const SCFVars& scf_vars,
   K = tamm_to_eigen_matrix<TensorType>(schwarz_mat);
   Tensor<TensorType>::deallocate(schwarz_mat);
 
-  // auto hf_t2   = std::chrono::high_resolution_clock::now();
-  // auto hf_time = std::chrono::duration_cast<std::chrono::duration<double>>((hf_t2 -
-  // hf_t1)).count(); if(ec.pg().rank() == 0)
-  //   std::cout << std::fixed << std::setprecision(2) << "Time to compute schwarz matrix: " <<
-  //   hf_time
-  //             << " secs" << endl;
-
   return K;
 }
 
@@ -1290,11 +1252,10 @@ template double gauxc_util::compute_xcf<double>(ExecutionContext& ec, const Syst
 
 #endif
 
-std::tuple<size_t, double, double>
+std::tuple<Matrix, size_t, double, double>
 gensqrtinv_atscf(ExecutionContext& ec, SystemData& sys_data, SCFVars& scf_vars,
-                 ScalapackInfo& scalapack_info, Tensor<double> S1, Tensor<double>& X_alpha,
-                 Tensor<double>& X_beta, TiledIndexSpace& tao_atom, bool symmetric,
-                 double threshold) {
+                 ScalapackInfo& scalapack_info, Tensor<double> S1, TiledIndexSpace& tao_atom,
+                 bool symmetric, double threshold) {
   using T = double;
 
   Scheduler sch{ec};
@@ -1345,21 +1306,14 @@ gensqrtinv_atscf(ExecutionContext& ec, SystemData& sys_data, SCFVars& scf_vars,
     }
   }
 
-  if(world_size > 1) {
-    // TODO: Should buffer this
-    ec.pg().broadcast(&n_cond, 0);
-    ec.pg().broadcast(&n_illcond, 0);
-    // ec.pg().broadcast( &condition_number,        0 );
-    ec.pg().broadcast(&result_condition_number, 0);
-  }
+  if(world_size > 1) { ec.pg().broadcast(&n_illcond, 0); }
+  n_cond = N - n_illcond;
 
   if(world_rank == 0) {
     // auto* V_cond = Vbuf + n_illcond * N;
     Matrix V_cond = V.block(n_illcond, 0, N - n_illcond, N);
     V.resize(0, 0);
-    X.resize(N, n_cond); // Xinv.resize( N, n_cond );
-    // Matrix V_cond(n_cond,N);
-    // V_cond = Eigen::Map<Matrix>(Vbuf + n_illcond * N,n_cond,N);
+    X.resize(N, n_cond);
     X = V_cond.transpose();
     V_cond.resize(0, 0);
   }
@@ -1377,31 +1331,21 @@ gensqrtinv_atscf(ExecutionContext& ec, SystemData& sys_data, SCFVars& scf_vars,
       // blas::scal( N, srt, Xinv_col, n_cond );
     }
 
-    if(symmetric) {
-      assert(not symmetric);
-      /*
-      // X is row major, thus we need to form X**T = V_cond * X**T
-      Matrix TMP = X;
-      X.resize( N, N );
-      blas::gemm( blas::Op::NoTrans, blas::Op::NoTrans, N, N, n_cond, 1., V_cond, N, TMP.data(),
-      n_cond, 0., X.data(), N );
-      */
-    }
   } // compute on root
 
-  auto nbf_new = N - n_illcond;
-
-  auto tAO_ortho = TiledIndexSpace{IndexSpace{range(0, (size_t) nbf_new)},
-                                   sys_data.options_map.scf_options.AO_tilesize};
+  TiledIndexSpace tAO_atom_ortho{IndexSpace{range((size_t) n_cond)},
+                                 sys_data.options_map.scf_options.AO_tilesize};
 
-  X_alpha = {tao_atom, tAO_ortho};
-  sch.allocate(X_alpha).execute();
-
-  if(world_rank == 0) eigen_to_tamm_tensor(X_alpha, X);
+  Tensor<T> x_tamm{tao_atom, tAO_atom_ortho};
+  sch.allocate(x_tamm).execute();
 
+  if(world_rank == 0) eigen_to_tamm_tensor(x_tamm, X);
   ec.pg().barrier();
 
-  return std::make_tuple(size_t(n_cond), condition_number, result_condition_number);
+  X = tamm_to_eigen_matrix(x_tamm);
+  sch.deallocate(x_tamm).execute();
+
+  return std::make_tuple(X, size_t(n_cond), condition_number, result_condition_number);
 }
 
 template Matrix read_scf_mat<double>(std::string matfile);
@@ -1433,14 +1377,6 @@ template void compute_density<double>(ExecutionContext& ec, const SystemData& sy
 template double tt_trace<double>(ExecutionContext& ec, Tensor<TensorType>& T1,
                                  Tensor<TensorType>& T2);
 
-// template Matrix compute_schwarz_ints<libint2::Operator>(ExecutionContext& ec, const SCFVars&
-// scf_vars,
-//                             const libint2::BasisSet& bs1, const libint2::BasisSet& _bs2,
-//                             bool use_2norm,
-//                             typename
-//                             libint2::operator_traits<libint2::Operator>::oper_params_type
-//                             params);
-
 template std::tuple<std::vector<int>, std::vector<int>, std::vector<int>>
 gather_task_vectors<double>(ExecutionContext& ec, std::vector<int>& s1vec, std::vector<int>& s2vec,
                             std::vector<int>& ntask_vec);
diff --git a/exachem/scf/scf_common.hpp b/exachem/scf/scf_common.hpp
index 66f442a..bca7ace 100644
--- a/exachem/scf/scf_common.hpp
+++ b/exachem/scf/scf_common.hpp
@@ -195,12 +195,6 @@ std::tuple<size_t, double, double> gensqrtinv(ExecutionContext& ec, SystemData&
                                               TAMMTensors& ttensors, bool symmetric = false,
                                               double threshold = 1e-5);
 
-// size_t nbasis(const libint2::BasisSet& shells) {
-//     size_t n = 0;
-//     for(const auto& shell : shells) n += shell.size();
-//     return n;
-// }
-
 inline size_t max_nprim(const libint2::BasisSet& shells) {
   size_t n = 0;
   for(auto shell: shells) n = std::max(shell.nprim(), n);
@@ -329,11 +323,10 @@ std::tuple<size_t, double, double> gensqrtinv(ExecutionContext& ec, SystemData&
                                               TAMMTensors& ttensors, bool symmetric,
                                               double threshold);
 
-std::tuple<size_t, double, double>
+std::tuple<Matrix, size_t, double, double>
 gensqrtinv_atscf(ExecutionContext& ec, SystemData& sys_data, SCFVars& scf_vars,
-                 ScalapackInfo& scalapack_info, Tensor<double> S1, Tensor<double>& X_alpha,
-                 Tensor<double>& X_beta, TiledIndexSpace& tao_atom, bool symmetric,
-                 double threshold);
+                 ScalapackInfo& scalapack_info, Tensor<double> S1, TiledIndexSpace& tao_atom,
+                 bool symmetric, double threshold);
 
 std::tuple<shellpair_list_t, shellpair_data_t> compute_shellpairs(const libint2::BasisSet& bs1,
                                                                   const libint2::BasisSet& _bs2,
@@ -357,23 +350,6 @@ gather_task_vectors(ExecutionContext& ec, std::vector<int>& s1vec, std::vector<i
 
 #if defined(USE_GAUXC)
 
-// template<typename T=double>
-// std::ostream& operator<<(std::ostream& os, const GauXC::Shell<T>& sh) {
-//   os << "GauXC::Shell:( O={" << sh.O()[0] << "," << sh.O()[1] << "," << sh.O()[2] << "}"
-//      << std::endl;
-//   os << "  ";
-//   os << " {l=" << sh.l() << ",sph=" << sh.pure() << "}";
-//   os << std::endl;
-
-//   for(auto i = 0ul; i < sh.nprim(); ++i) {
-//     os << "  " << sh.alpha()[i];
-//     os << " " << sh.coeff().at(i);
-//     os << std::endl;
-//   }
-
-//   return os;
-// }
-
 namespace gauxc_util {
 
 GauXC::Molecule make_gauxc_molecule(const std::vector<libint2::Atom>& atoms);
diff --git a/exachem/scf/scf_guess.cpp b/exachem/scf/scf_guess.cpp
index 2fcd552..f5a96d4 100644
--- a/exachem/scf/scf_guess.cpp
+++ b/exachem/scf/scf_guess.cpp
@@ -127,7 +127,6 @@ void compute_dipole_ints(ExecutionContext& ec, const SCFVars& spvars, Tensor<Ten
   Engine engine(otype, max_nprim(shells), max_l(shells), 0);
 
   // engine.set(otype);
-
   // auto& buf = (engine.results());
 
   auto compute_dipole_ints_lambda = [&](const IndexVector& blockid) {
@@ -142,9 +141,6 @@ void compute_dipole_ints(ExecutionContext& ec, const SCFVars& spvars, Tensor<Ten
 
     auto bd1 = block_dims[1];
 
-    // cout << "blockid: [" << blockid[0] <<"," << blockid[1] << "], dims(0,1) = " <<
-    //  block_dims[0] << ", " << block_dims[1] << endl;
-
     // auto s1 = blockid[0];
     auto                  s1range_end   = shell_tile_map[bi0];
     decltype(s1range_end) s1range_start = 0l;
@@ -191,10 +187,6 @@ void compute_dipole_ints(ExecutionContext& ec, const SCFVars& spvars, Tensor<Ten
         // "map" buffer to a const Eigen Matrix, and copy it to the
         // corresponding blocks of the result
 
-        // cout << buf[1].size() << endl;
-        // cout << buf[2].size() << endl;
-        // cout << buf[3].size() << endl;
-
         Eigen::Map<const Matrix> buf_mat_X(buf[1], n1, n2);
         Eigen::Map<Matrix>(&tbufX[0], n1, n2) = buf_mat_X;
         Eigen::Map<const Matrix> buf_mat_Y(buf[2], n1, n2);
@@ -202,10 +194,6 @@ void compute_dipole_ints(ExecutionContext& ec, const SCFVars& spvars, Tensor<Ten
         Eigen::Map<const Matrix> buf_mat_Z(buf[3], n1, n2);
         Eigen::Map<Matrix>(&tbufZ[0], n1, n2) = buf_mat_Z;
 
-        // cout << "buf_mat_X :" << buf_mat_X << endl;
-        // cout << "buf_mat_Y :" << buf_mat_Y << endl;
-        // cout << "buf_mat_Z :" << buf_mat_Z << endl;
-
         auto curshelloffset_i = 0U;
         auto curshelloffset_j = 0U;
         for(auto x = s1range_start; x < s1; x++) curshelloffset_i += AO_tiles[x];
@@ -270,9 +258,6 @@ void compute_1body_ints(ExecutionContext& ec, const SCFVars& scf_vars, Tensor<Te
 
     auto bd1 = block_dims[1];
 
-    // cout << "blockid: [" << blockid[0] <<"," << blockid[1] << "], dims(0,1) = " <<
-    //  block_dims[0] << ", " << block_dims[1] << endl;
-
     // auto s1 = blockid[0];
     auto                  s1range_end   = shell_tile_map[bi0];
     decltype(s1range_end) s1range_start = 0l;
@@ -309,8 +294,6 @@ void compute_1body_ints(ExecutionContext& ec, const SCFVars& scf_vars, Tensor<Te
         auto n2 = shells[s2].size();
 
         std::vector<TensorType> tbuf(n1 * n2);
-        // cout << "s1,s2,n1,n2 = "  << s1 << "," << s2 <<
-        //       "," << n1 <<"," << n2 <<endl;
 
         // compute shell pair; return is the pointer to the buffer
         engine.compute(shells[s1], shells[s2]);
@@ -330,10 +313,6 @@ void compute_1body_ints(ExecutionContext& ec, const SCFVars& scf_vars, Tensor<Te
         auto   dimi = curshelloffset_i + AO_tiles[s1];
         auto   dimj = curshelloffset_j + AO_tiles[s2];
 
-        // cout << "curshelloffset_i,curshelloffset_j,dimi,dimj = "  << curshelloffset_i << "," <<
-        // curshelloffset_j <<
-        //       "," << dimi <<"," << dimj <<endl;
-
         for(size_t i = curshelloffset_i; i < dimi; i++) {
           for(size_t j = curshelloffset_j; j < dimj; j++, c++) { dbuf[i * bd1 + j] = tbuf[c]; }
         }
@@ -389,9 +368,6 @@ void compute_ecp_ints(ExecutionContext& ec, const SCFVars& scf_vars, Tensor<Tens
 
     auto bd1 = block_dims[1];
 
-    // cout << "blockid: [" << blockid[0] <<"," << blockid[1] << "], dims(0,1) = " <<
-    //  block_dims[0] << ", " << block_dims[1] << endl;
-
     // auto s1 = blockid[0];
     auto                  s1range_end   = shell_tile_map[bi0];
     decltype(s1range_end) s1range_start = 0l;
@@ -519,9 +495,6 @@ void compute_pchg_ints(ExecutionContext& ec, const SCFVars& scf_vars, Tensor<Ten
 
     auto bd1 = block_dims[1];
 
-    // cout << "blockid: [" << blockid[0] <<"," << blockid[1] << "], dims(0,1) = " <<
-    //  block_dims[0] << ", " << block_dims[1] << endl;
-
     // auto s1 = blockid[0];
     auto                  s1range_end   = shell_tile_map[bi0];
     decltype(s1range_end) s1range_start = 0l;
diff --git a/exachem/scf/scf_guess.hpp b/exachem/scf/scf_guess.hpp
index 6f702ce..c43327c 100644
--- a/exachem/scf/scf_guess.hpp
+++ b/exachem/scf/scf_guess.hpp
@@ -25,10 +25,9 @@ void subshell_occvec(Real& occvec, size_t size, size_t& ne);
 
 /// @brief computes average orbital occupancies in the ground state of a neutral
 ///        atoms
-/// @throw if Z > 53
 /// @return occupation vector corresponding to the ground state electronic
 ///         configuration of a neutral atom with atomic number \c Z
-///         corresponding to the orbital ordering in STO-3G basis
+///         corresponding to the orbital ordering.
 template<typename Real = double>
 const std::vector<Real> compute_ao_occupation_vector(size_t Z);
 
diff --git a/exachem/scf/scf_iter.cpp b/exachem/scf/scf_iter.cpp
index f1766c7..f7a4c61 100644
--- a/exachem/scf/scf_iter.cpp
+++ b/exachem/scf/scf_iter.cpp
@@ -515,11 +515,15 @@ void compute_2bf_ri(ExecutionContext& ec, ScalapackInfo& scalapack_info, const S
 
     block_for(ec, K_tamm(), compute_2body_2index_ints_lambda);
 
-    Matrix          V;
-    Eigen::VectorXd eps(ndf);
+    Matrix                  V;
+    std::vector<TensorType> eps(ndf);
 
     auto ig1 = std::chrono::high_resolution_clock::now();
 
+    Tensor<TensorType> v_tmp{scf_vars.tdfAO, scf_vars.tdfAO};
+    Tensor<TensorType> eps_tamm{scf_vars.tdfAO};
+    Tensor<TensorType>::allocate(&ec, v_tmp, eps_tamm);
+
 #if defined(USE_SCALAPACK)
     Tensor<TensorType> V_sca;
     if(scalapack_info.comm != MPI_COMM_NULL) {
@@ -545,7 +549,6 @@ void compute_2bf_ri(ExecutionContext& ec, ScalapackInfo& scalapack_info, const S
       if(grid.ipr() >= 0 and grid.ipc() >= 0) {
         auto desc_S = desc_lambda(ndf, ndf);
         auto desc_V = desc_lambda(ndf, ndf);
-        // scalapackpp::BlockCyclicMatrix<TensorType> V_sca(grid, N, N, mb, mb);
 
         /*info=*/scalapackpp::hereig(scalapackpp::Job::Vec, scalapackpp::Uplo::Lower, desc_S[2],
                                      S_BC.access_local_buf(), 1, 1, desc_S, eps.data(),
@@ -553,8 +556,7 @@ void compute_2bf_ri(ExecutionContext& ec, ScalapackInfo& scalapack_info, const S
       }
 
       Tensor<TensorType>::deallocate(S_BC);
-      if(scalapack_info.pg.rank() == 0) V = tamm_to_eigen_matrix<TensorType>(V_sca);
-      Tensor<TensorType>::deallocate(V_sca);
+      tamm::from_block_cyclic_tensor(V_sca, v_tmp);
     }
 
 #else
@@ -563,29 +565,20 @@ void compute_2bf_ri(ExecutionContext& ec, ScalapackInfo& scalapack_info, const S
       tamm_to_eigen_tensor(K_tamm, V);
 
       lapack::syevd(lapack::Job::Vec, lapack::Uplo::Lower, ndf, V.data(), ndf, eps.data());
+      tamm::eigen_to_tamm_tensor(v_tmp, V);
     }
 #endif
 
     if(rank == 0) {
-      Matrix Vp = V;
-
-      for(int j = 0; j < ndf; ++j) {
-        double tmp = 1.0 / sqrt(eps(j));
-        blas::scal(ndf, tmp, Vp.data() + j * ndf, 1);
-      }
-
-      Matrix ke(ndf, ndf);
-      blas::gemm(blas::Layout::ColMajor, blas::Op::NoTrans, blas::Op::Trans, ndf, ndf, ndf, 1,
-                 Vp.data(), ndf, V.data(), ndf, 0, ke.data(), ndf);
-      eigen_to_tamm_tensor(K_tamm, ke);
-      V.resize(0, 0);
+      std::transform(eps.begin(), eps.end(), eps.begin(),
+                     [](auto& c) { return 1.0 / std::sqrt(c); });
+      tamm::vector_to_tamm_tensor(eps_tamm, eps);
     }
-    eps.resize(0);
-    // sch
-    // (k_tmp(d_mu,d_nu) = V_tamm(d_ku,d_mu) * V_tamm(d_ku,d_nu))
-    // (K_tamm(d_mu,d_nu) = eps_tamm(d_mu,d_ku) * k_tmp(d_ku,d_nu)) .execute();
+    ec.pg().barrier();
 
-    // Tensor<TensorType>::deallocate(k_tmp, eps_tamm, V_tamm);
+    sch(K_tamm(d_mu, d_nu) = v_tmp(d_nu, d_mu) * eps_tamm(d_nu))
+      .deallocate(v_tmp, eps_tamm)
+      .execute();
 
     auto ig2    = std::chrono::high_resolution_clock::now();
     auto igtime = std::chrono::duration_cast<std::chrono::duration<double>>((ig2 - ig1)).count();
@@ -593,8 +586,6 @@ void compute_2bf_ri(ExecutionContext& ec, ScalapackInfo& scalapack_info, const S
     if(rank == 0 && debug)
       std::cout << std::fixed << std::setprecision(2) << "V^-1/2: " << igtime << "s, ";
 
-    // contract(1.0, Zxy, {1, 2, 3}, K, {1, 4}, 0.0, xyK, {2, 3, 4});
-    // Tensor3D xyK = Zxy.contract(K,aidx_00);
     sch(xyK(mu, nu, d_nu) = Zxy(d_mu, mu, nu) * K_tamm(d_mu, d_nu))
       .deallocate(K_tamm, Zxy) // release memory Zxy
       .execute(exhw);
@@ -711,7 +702,6 @@ void compute_2bf(ExecutionContext& ec, ScalapackInfo& scalapack_info, const Syst
   engine.set_precision(engine_precision);
   const auto& buf = engine.results();
 
-#if 1
   auto comp_2bf_lambda = [&](IndexVector blockid) {
     auto s1        = blockid[0];
     auto bf1_first = shell2bf[s1];
@@ -847,7 +837,6 @@ void compute_2bf(ExecutionContext& ec, ScalapackInfo& scalapack_info, const Syst
       }
     }
   };
-#endif
 
   decltype(do_t1) do_t2;
   double          do_time;
diff --git a/exachem/scf/scf_main.cpp b/exachem/scf/scf_main.cpp
index 28a8c10..f31d234 100644
--- a/exachem/scf/scf_main.cpp
+++ b/exachem/scf/scf_main.cpp
@@ -600,14 +600,17 @@ hartree_fock(ExecutionContext& exc, const string filename, OptionsMap options_ma
     // pre-compute data for Schwarz bounds
     std::string schwarz_matfile = files_prefix + ".schwarz";
     Matrix      SchwarzK;
-    if(N >= restart_size && fs::exists(schwarz_matfile)) {
-      if(rank == 0) cout << "Read Schwarz matrix from disk ... " << endl;
-      SchwarzK = read_scf_mat<TensorType>(schwarz_matfile);
-    }
-    else {
-      // if(rank == 0) cout << "pre-computing data for Schwarz bounds... " << endl;
-      SchwarzK = compute_schwarz_ints<>(ec, scf_vars, shells);
-      if(rank == 0) write_scf_mat<TensorType>(SchwarzK, schwarz_matfile);
+
+    if(!do_density_fitting) {
+      if(N >= restart_size && fs::exists(schwarz_matfile)) {
+        if(rank == 0) cout << "Read Schwarz matrix from disk ... " << endl;
+        SchwarzK = read_scf_mat<TensorType>(schwarz_matfile);
+      }
+      else {
+        // if(rank == 0) cout << "pre-computing data for Schwarz bounds... " << endl;
+        SchwarzK = compute_schwarz_ints<>(ec, scf_vars, shells);
+        if(rank == 0) write_scf_mat<TensorType>(SchwarzK, schwarz_matfile);
+      }
     }
 
     hf_t1 = std::chrono::high_resolution_clock::now();
diff --git a/methods/CMakeLists.txt b/methods/CMakeLists.txt
index 7f98a65..3e2ebbe 100644
--- a/methods/CMakeLists.txt
+++ b/methods/CMakeLists.txt
@@ -6,6 +6,8 @@ include(TargetMacros)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/..)
 set(EXACHEM_SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/../exachem)
 
+option(EC_COMPLEX "enable complex codes" ON)
+
 #TODO: cleanup
 include_directories(${EXACHEM_SRC_DIR})
 include_directories(${EXACHEM_SRC_DIR}/util/external)
@@ -22,9 +24,9 @@ include(${EXACHEM_SRC_DIR}/cc/ccsd_t/ccsd_t.cmake)
 include(${EXACHEM_SRC_DIR}/cc/lambda/lambda.cmake)
 include(${EXACHEM_SRC_DIR}/cc/ducc/ducc.cmake)
 
-if(NOT USE_UPCXX)
-include(${EXACHEM_SRC_DIR}/cc/gfcc/gfcc.cmake)
-include(${EXACHEM_SRC_DIR}/cc/rteom/rteom.cmake)
+if(NOT USE_UPCXX AND EC_COMPLEX)
+  include(${EXACHEM_SRC_DIR}/cc/gfcc/gfcc.cmake)
+  include(${EXACHEM_SRC_DIR}/cc/rteom/rteom.cmake)
 endif()
 
 set(CoupledCluster_SRCS ${CC_SRCS} ${CC2_SRCS} ${CCSD_SRCS} ${CCSD_T_SRCS} 
diff --git a/methods/ExaChem.cpp b/methods/ExaChem.cpp
index 25c1f43..be059be 100644
--- a/methods/ExaChem.cpp
+++ b/methods/ExaChem.cpp
@@ -13,7 +13,9 @@
 #include "exachem/cc/eom/eomccsd_opt.hpp"
 // clang-format on
 
-#if !defined(USE_UPCXX)
+#define EC_COMPLEX
+
+#if !defined(USE_UPCXX) and defined(EC_COMPLEX)
 void gfccsd_driver(std::string filename, OptionsMap options_map);
 void rt_eom_cd_ccsd_driver(std::string filename, OptionsMap options_map);
 #include "exachem/fci/fci.hpp"
@@ -107,13 +109,15 @@ int main(int argc, char* argv[]) {
   else if(task.ccsd_lambda) ccsd_lambda_driver(filename, options_map);
   else if(task.eom_ccsd) eom_ccsd_driver(filename, options_map);
   else if(task.ducc) ducc_driver(filename, options_map);
-#if !defined(USE_UPCXX)
+#if !defined(USE_UPCXX) and defined(EC_COMPLEX)
   else if(task.fci || task.fcidump) fci_driver(filename, options_map);
   else if(task.gfccsd) gfccsd_driver(filename, options_map);
   else if(task.rteom_ccsd) rt_eom_cd_ccsd_driver(filename, options_map);
 #endif
 
-  else tamm_terminate("[INPUT FILE ERROR] No task specified!");
+  else
+    tamm_terminate(
+      "[ERROR] Unsupported task specified (or) code for the specified task is not built");
 
   tamm::finalize();