From d7a5d31b5b4d836f45168f80d681dd5a06de1cdc Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Tue, 28 Nov 2023 12:44:24 -0500 Subject: [PATCH 1/2] Generate the density matrix if needed for the electron density surface Fix #1480 and the bug report on the forum: https://discuss.avogadro.cc/t/visualize-electron-density-from-orca-molden-file/3750 Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 9 ++++++++- .../qtplugins/surfaces/gaussiansetconcurrent.cpp | 12 +++++++++--- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index 72382090f5..b15b8df117 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -55,7 +55,13 @@ double GaussianSetTools::calculateMolecularOrbital(const Vector3& position, } bool GaussianSetTools::calculateElectronDensity(Cube& cube) const -{ +{ + const MatrixX& matrix = m_basis->densityMatrix(); + if (matrix.rows() == 0 || matrix.cols() == 0) { + // we don't have a density matrix, so generate one + m_basis->generateDensityMatrix(); + } + for (size_t i = 0; i < cube.data()->size(); ++i) { Vector3 pos = cube.position(i); cube.setValue(i, calculateElectronDensity(pos)); @@ -67,6 +73,7 @@ double GaussianSetTools::calculateElectronDensity(const Vector3& position) const { const MatrixX& matrix = m_basis->densityMatrix(); int matrixSize(static_cast(m_basis->moMatrix().rows())); + if (matrix.rows() != matrixSize || matrix.cols() != matrixSize) { return 0.0; } diff --git a/avogadro/qtplugins/surfaces/gaussiansetconcurrent.cpp b/avogadro/qtplugins/surfaces/gaussiansetconcurrent.cpp index 7c3771fee1..099fc7dacd 100644 --- a/avogadro/qtplugins/surfaces/gaussiansetconcurrent.cpp +++ b/avogadro/qtplugins/surfaces/gaussiansetconcurrent.cpp @@ -56,8 +56,8 @@ void GaussianSetConcurrent::setMolecule(Core::Molecule* mol) if (!mol) return; m_set = dynamic_cast(mol->basisSet()); - - delete m_tools; + + delete m_tools; m_tools = new GaussianSetTools(mol); } @@ -76,6 +76,12 @@ bool GaussianSetConcurrent::calculateMolecularOrbital(Core::Cube* cube, bool GaussianSetConcurrent::calculateElectronDensity(Core::Cube* cube) { + const MatrixX& matrix = m_set->densityMatrix(); + if (matrix.rows() == 0 || matrix.cols() == 0) { + // we don't have a density matrix, so calculate one + m_set->generateDensityMatrix(); + } + return setUpCalculation(cube, 0, GaussianSetConcurrent::processDensity); } @@ -141,4 +147,4 @@ void GaussianSetConcurrent::processSpinDensity(GaussianShell& shell) Vector3 pos = shell.tCube->position(shell.pos); shell.tCube->setValue(shell.pos, shell.tools->calculateSpinDensity(pos)); } -} +} // namespace Avogadro::QtPlugins From 2a8b6d4d69ff26eb0a4fe64a4761ab6693b14a1f Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Tue, 28 Nov 2023 13:31:31 -0500 Subject: [PATCH 2/2] Fix formatting Signed-off-by: Geoff Hutchison --- avogadro/core/gaussiansettools.cpp | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/avogadro/core/gaussiansettools.cpp b/avogadro/core/gaussiansettools.cpp index b15b8df117..140032f980 100644 --- a/avogadro/core/gaussiansettools.cpp +++ b/avogadro/core/gaussiansettools.cpp @@ -11,7 +11,6 @@ #include - using std::vector; namespace Avogadro::Core { @@ -22,9 +21,7 @@ GaussianSetTools::GaussianSetTools(Molecule* mol) : m_molecule(mol) m_basis = dynamic_cast(m_molecule->basisSet()); } -GaussianSetTools::~GaussianSetTools() -{ -} +GaussianSetTools::~GaussianSetTools() {} bool GaussianSetTools::calculateMolecularOrbital(Cube& cube, int moNumber) const { @@ -55,7 +52,7 @@ double GaussianSetTools::calculateMolecularOrbital(const Vector3& position, } bool GaussianSetTools::calculateElectronDensity(Cube& cube) const -{ +{ const MatrixX& matrix = m_basis->densityMatrix(); if (matrix.rows() == 0 || matrix.cols() == 0) { // we don't have a density matrix, so generate one @@ -160,7 +157,7 @@ inline vector GaussianSetTools::calculateValues( // Calculate the deltas for the position for (Index i = 0; i < atomsSize; ++i) { deltas.emplace_back(pos - - (m_molecule->atom(i).position3d() * ANGSTROM_TO_BOHR)); + (m_molecule->atom(i).position3d() * ANGSTROM_TO_BOHR)); dr2.push_back(deltas[i].squaredNorm()); } @@ -253,7 +250,7 @@ inline void GaussianSetTools::pointD(unsigned int moIndex, const Vector3& delta, i < m_basis->gtoIndices()[moIndex + 1]; ++i) { // Calculate the common factor double tmpGTO = exp(-gtoA[i] * dr2); - for (double & component : components) + for (double& component : components) component += gtoCN[cIndex++] * tmpGTO; } @@ -286,7 +283,7 @@ inline void GaussianSetTools::pointD5(unsigned int moIndex, i < m_basis->gtoIndices()[moIndex + 1]; ++i) { // Calculate the common factor double tmpGTO = exp(-gtoA[i] * dr2); - for (double & component : components) + for (double& component : components) component += gtoCN[cIndex++] * tmpGTO; } @@ -325,7 +322,7 @@ inline void GaussianSetTools::pointF(unsigned int moIndex, const Vector3& delta, i < m_basis->gtoIndices()[moIndex + 1]; ++i) { // Calculate the common factor double tmpGTO = exp(-gtoA[i] * dr2); - for (double & component : components) + for (double& component : components) component += gtoCN[cIndex++] * tmpGTO; } @@ -343,16 +340,7 @@ inline void GaussianSetTools::pointF(unsigned int moIndex, const Vector3& delta, double componentsF[10] = { // molden order // e.g https://gau2grid.readthedocs.io/en/latest/order.html - xxx, - yyy, - zzz, - xyy, - xxy, - xxz, - xzz, - yzz, - yyz, - xyz + xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz }; for (int i = 0; i < 10; ++i) @@ -378,7 +366,7 @@ inline void GaussianSetTools::pointF7(unsigned int moIndex, i < m_basis->gtoIndices()[moIndex + 1]; ++i) { // Calculate the common factor double tmpGTO = exp(-gtoA[i] * dr2); - for (double & component : components) + for (double& component : components) component += gtoCN[cIndex++] * tmpGTO; } @@ -425,4 +413,4 @@ final normalization values[baseIndex + i] += components[i] * componentsF[i]; } -} // End Avogadro namespace +} // namespace Avogadro::Core