Generate the density matrix if needed for the electron density surface

Fix OpenChemistry#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 <[email protected]>
ghutchis · Nov 28, 2023 · 4e52a46 · 4e52a46 · github-actions · Nov 28, 2023
1 parent 57e0dd6
commit 4e52a46
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 4 deletions.
diff --git 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<int>(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
@@ -56,8 +56,8 @@ void GaussianSetConcurrent::setMolecule(Core::Molecule* mol)
   if (!mol)
     return;
   m_set = dynamic_cast<GaussianSet*>(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