From 06471ba37b3754eba34feedaf97b254d3e0a74b1 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 28 Dec 2023 20:28:23 -0500 Subject: [PATCH] Fix quantum surface max cutoff for diffuse functions Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 5 +++-- avogadro/quantumio/gaussianfchk.cpp | 14 ++++++++++---- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index 1149a60ea3..08d255c235 100644 --- 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& 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 index 2daf34a2f5..a8b881f326 100644 --- 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"; }