Merge pull request #1556 from ghutchis/max-quantum-cutoff

Fix quantum surface max cutoff for diffuse functions
OpenChemistry · Dec 29, 2023 · a68514a · a68514a
2 parents 52aaf9f + 06471ba
commit a68514a
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 6 deletions.
diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp
@@ -146,7 +146,7 @@ inline void GaussianSetTools::calculateCutoffs()
   // .. 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
+  const double maxDistance = 100.0;      // just in case (diffuse functions)
 
   // get the exponents and normalized coefficients
   const std::vector<double>& exponents = m_basis->gtoA();
@@ -167,7 +167,8 @@ inline void GaussianSetTools::calculateCutoffs()
     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));
+      // except for S, we don't want to start at the origin
+      double r = std::min(maxDistance, std::sqrt(L / (2 * alpha)));
       double value = coeff * std::pow(r, L) * std::exp(-alpha * r * r);
 
       while (value > threshold && r < maxDistance) {

diff --git a/avogadro/quantumio/gaussianfchk.cpp b/avogadro/quantumio/gaussianfchk.cpp
@@ -311,14 +311,20 @@ void GaussianFchk::load(GaussianSet* basis)
       basis->setMolecularOrbitals(m_alphaMOcoeffs, BasisSet::Alpha);
     if (m_betaMOcoeffs.size())
       basis->setMolecularOrbitals(m_betaMOcoeffs, BasisSet::Beta);
+
     if (m_density.rows())
       basis->setDensityMatrix(m_density);
     if (m_spinDensity.rows())
       basis->setSpinDensityMatrix(m_spinDensity);
-    if (m_alphaOrbitalEnergy.size())
-      basis->setMolecularOrbitalEnergy(m_alphaOrbitalEnergy, BasisSet::Alpha);
-    if (m_betaOrbitalEnergy.size())
-      basis->setMolecularOrbitalEnergy(m_betaOrbitalEnergy, BasisSet::Beta);
+
+    if (m_orbitalEnergy.size()) // restricted calculation
+      basis->setMolecularOrbitalEnergy(m_orbitalEnergy);
+    else {
+      if (m_alphaOrbitalEnergy.size())
+        basis->setMolecularOrbitalEnergy(m_alphaOrbitalEnergy, BasisSet::Alpha);
+      if (m_betaOrbitalEnergy.size())
+        basis->setMolecularOrbitalEnergy(m_betaOrbitalEnergy, BasisSet::Beta);
+    }
   } else {
     cout << "Basis set is not valid!\n";
   }