From 84cb4b0dd4eda61065d6998e945d9864b27ff238 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sun, 24 Dec 2023 23:20:16 -0500 Subject: [PATCH 1/4] Skip calculating points too far apart (e.g., negligible) Rough timing on C60 shows significant speedup ~2x on C60 Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index 140032f980..98b446d615 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -132,7 +132,7 @@ bool GaussianSetTools::isValid() const inline bool GaussianSetTools::isSmall(double val) const { - if (val > -1e-20 && val < 1e-20) + if (val > -1e-12 && val < 1e-12) return true; else return false; @@ -168,6 +168,12 @@ inline vector 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) + continue; + switch (basis[i]) { case GaussianSet::S: pointS(i, dr2[atomIndices[i]], values); From 94f8a95fe9ee21eb960d6b527c07845249118261 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 28 Dec 2023 00:25:20 -0500 Subject: [PATCH 2/4] Estimate cutoff distance based on current basis set Signed-off-by: Geoff Hutchison --- 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(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& exponents = m_basis->gtoA(); + const std::vector& coefficients = m_basis->gtoCN(); + const std::vector& 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 GaussianSetTools::calculateValues( @@ -169,9 +212,8 @@ inline vector 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 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& values) const; void pointF7(unsigned int index, const Vector3& delta, double dr2, std::vector& 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 From 86003779e41e28e147e2306effc5f1ad410751d5 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 28 Dec 2023 00:31:32 -0500 Subject: [PATCH 3/4] Fix typo Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index 4ab0ec5c38..ea078859b3 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -161,7 +161,7 @@ inline void GaussianSetTools::calculateCutoffs() // .. 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); + 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) From 521cfa576153006b3e5f6153f44ee54542973996 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 28 Dec 2023 00:32:25 -0500 Subject: [PATCH 4/4] Fix typo again Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index ea078859b3..1149a60ea3 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -161,7 +161,6 @@ inline void GaussianSetTools::calculateCutoffs() // .. 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) @@ -169,11 +168,11 @@ inline void GaussianSetTools::calculateCutoffs() 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); + double value = coeff * std::pow(r, L) * std::exp(-alpha * r * r); while (value > threshold && r < maxDistance) { r += 0.25; - value = prefactor * std::exp(-alpha * r * r); + value = coeff * std::pow(r, L) * std::exp(-alpha * r * r); } m_cutoffDistances[L] = std::max(m_cutoffDistances[L], r * r);