From 94f8a95fe9ee21eb960d6b527c07845249118261 Mon Sep 17 00:00:00 2001
From: Geoff Hutchison <geoff.hutchison@gmail.com>
Date: Thu, 28 Dec 2023 00:25:20 -0500
Subject: [PATCH] Estimate cutoff distance based on current basis set

Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
---
 avogadro/core/gaussiansettools.cpp | 58 +++++++++++++++++++++++++-----
 avogadro/core/gaussiansettools.h   | 13 +++++--
 2 files changed, 61 insertions(+), 10 deletions(-)

diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp
index 98b446d615..4ab0ec5c38 100644
--- a/avogadro/core/gaussiansettools.cpp
+++ b/avogadro/core/gaussiansettools.cpp
@@ -17,8 +17,11 @@ namespace Avogadro::Core {
 
 GaussianSetTools::GaussianSetTools(Molecule* mol) : m_molecule(mol)
 {
-  if (m_molecule)
+  if (m_molecule) {
     m_basis = dynamic_cast<GaussianSet*>(m_molecule->basisSet());
+    m_cutoffDistances.resize(7, 0.0); // s, p, d, f, g, h, i (for now)
+    calculateCutoffs();
+  }
 }
 
 GaussianSetTools::~GaussianSetTools() {}
@@ -132,10 +135,50 @@ bool GaussianSetTools::isValid() const
 
 inline bool GaussianSetTools::isSmall(double val) const
 {
-  if (val > -1e-12 && val < 1e-12)
-    return true;
-  else
-    return false;
+  return std::abs(val) < 1e-12;
+}
+
+inline void GaussianSetTools::calculateCutoffs()
+{
+  // Guesstimate a distance we can ignore the exp(-alpha * r^2) term
+  // .. because it's negligible
+  // This will depend on the angular momentum of the basis function
+  // .. so we calculate it for whatever L values in this basis set
+
+  const double threshold = 0.03 * 0.001; // 0.1% of a typical isovalue
+  const double maxDistance = 15.0;       // 15 Angstroms
+
+  // get the exponents and normalized coefficients
+  const std::vector<double>& exponents = m_basis->gtoA();
+  const std::vector<double>& coefficients = m_basis->gtoCN();
+  const std::vector<int>& sym = m_basis->symmetry();
+
+  // we loop through the "symmetry" (i.e., L values in this basis set)
+  for (size_t i = 0; i < sym.size(); ++i) {
+    int L = symToL[sym[i]];
+
+    // this is a hack, since not all coefficients will be the same
+    // .. but it's a good approximation since they'll be similar
+    unsigned int cIndex = m_basis->cIndices()[i];
+    const double coeff = std::abs(coefficients[cIndex]);
+    const double preFactor = coeff * std::pow(r, L);
+
+    // now loop through all exponents for this L value
+    // (e.g., multiple terms - we don't know which is the most diffuse)
+    for (unsigned int j = m_basis->gtoIndices()[i];
+         j < m_basis->gtoIndices()[i + 1]; ++j) {
+      double alpha = exponents[j];
+      double r = std::sqrt(L / (2 * alpha));
+      double value = prefactor * std::exp(-alpha * r * r);
+
+      while (value > threshold && r < maxDistance) {
+        r += 0.25;
+        value = prefactor * std::exp(-alpha * r * r);
+      }
+
+      m_cutoffDistances[L] = std::max(m_cutoffDistances[L], r * r);
+    }
+  }
 }
 
 inline vector<double> GaussianSetTools::calculateValues(
@@ -169,9 +212,8 @@ inline vector<double> GaussianSetTools::calculateValues(
   // Now calculate the values at this point in space
   for (unsigned int i = 0; i < basisSize; ++i) {
     // bail early if the distance is too big
-    // TODO: this should be smarter and use the exponents
-    //   .. and angular momentum
-    if (dr2[atomIndices[i]] > 100.0)
+    double cutoff = m_cutoffDistances[symToL[basis[i]]];
+    if (dr2[atomIndices[i]] > cutoff)
       continue;
 
     switch (basis[i]) {
diff --git a/avogadro/core/gaussiansettools.h b/avogadro/core/gaussiansettools.h
index 0e3b69be5d..7da89d44c8 100644
--- a/avogadro/core/gaussiansettools.h
+++ b/avogadro/core/gaussiansettools.h
@@ -99,9 +99,14 @@ class AVOGADROCORE_EXPORT GaussianSetTools
   Molecule* m_molecule;
   GaussianSet* m_basis;
   BasisSet::ElectronType m_type = BasisSet::Paired;
+  std::vector<double> m_cutoffDistances;
 
   bool isSmall(double value) const;
 
+  // get the cutoff distance for the given angular momentum
+  // .. and the current basis set (exponents)
+  void calculateCutoffs();
+
   /**
    * @brief Calculate the values at this position in space. The public calculate
    * functions call this function to prepare values before multiplying by the
@@ -122,9 +127,13 @@ class AVOGADROCORE_EXPORT GaussianSetTools
               std::vector<double>& values) const;
   void pointF7(unsigned int index, const Vector3& delta, double dr2,
                std::vector<double>& values) const;
+
+  // map from symmetry to angular momentum
+  // S, SP, P, D, D5, F, F7, G, G9, etc.
+  const int symToL[13] = { 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6 };
 };
 
-} // End Core namespace
-} // End Avogadro namespace
+} // namespace Core
+} // namespace Avogadro
 
 #endif // AVOGADRO_CORE_GAUSSIANSETTOOLS_H