Merge pull request #1803 from ghutchis/homo-lumo-properties

Add HOMO and LUMO energies to properties
OpenChemistry · Nov 21, 2024 · 9e353ed · 9e353ed
2 parents 3893211 + 07b10b7
commit 9e353ed
Show file tree

Hide file tree

Showing 5 changed files with 34 additions and 11 deletions.
diff --git a/avogadro/qtplugins/dipole/dipole.h b/avogadro/qtplugins/dipole/dipole.h
@@ -39,7 +39,7 @@ class Dipole : public QtGui::ScenePlugin
   }
 
   QWidget* setupWidget() override;
-  bool hasSetupWidget() const override { return true; }
+  bool hasSetupWidget() const override { return false; }
 
 private:
   std::string m_name = "Dipole Moment";

diff --git a/avogadro/qtplugins/molecularproperties/molecularmodel.cpp b/avogadro/qtplugins/molecularproperties/molecularmodel.cpp
@@ -32,6 +32,10 @@ using Avogadro::Core::GaussianSet;
 using Avogadro::QtGui::Molecule;
 using QtGui::Molecule;
 
+// CODATA 2022
+// https://physics.nist.gov/cgi-bin/cuu/Value?hrev
+#define AU_TO_EV 27.211386245981
+
 MolecularModel::MolecularModel(QObject* parent)
   : QAbstractTableModel(parent), m_molecule(nullptr)
 {
@@ -240,7 +244,7 @@ QVariant MolecularModel::data(const QModelIndex& index, int role) const
     return QVariant::fromValue(m_molecule->residueCount());
   else if (key == " 9totalCharge")
     return QVariant::fromValue(m_molecule->totalCharge());
-  else if (key == " 10totalSpinMultiplicity")
+  else if (key == " 9totalSpinMultiplicity")
     return QVariant::fromValue(m_molecule->totalSpinMultiplicity());
 
   return QString::fromStdString(it->second.toString());
@@ -290,7 +294,7 @@ QVariant MolecularModel::headerData(int section, Qt::Orientation orientation,
       return tr("Number of Chains");
     else if (it->first == " 9totalCharge")
       return tr("Net Charge");
-    else if (it->first == " 10totalSpinMultiplicity")
+    else if (it->first == " 9totalSpinMultiplicity")
       return tr("Net Spin Multiplicity");
     else if (it->first == "dipoleMoment")
       return tr("Dipole Moment (Debye)");
@@ -390,16 +394,33 @@ void MolecularModel::updateTable(unsigned int flags)
   if (m_molecule->totalCharge() != 0)
     m_propertiesCache.setValue(" 9totalCharge", m_molecule->totalCharge());
   if (m_molecule->totalSpinMultiplicity() != 1)
-    m_propertiesCache.setValue(" 10totalSpinMultiplicity",
+    m_propertiesCache.setValue(" 9totalSpinMultiplicity",
                                m_molecule->totalSpinMultiplicity());
   if (m_molecule->hasData("dipoleMoment")) {
     auto dipole = m_molecule->data("dipoleMoment").toVector3();
-    m_propertiesCache.setValue("dipoleMoment", dipole.norm());
+    QString moment = QString::number(dipole.norm(), 'f', 3);
+    m_propertiesCache.setValue("dipoleMoment", moment.toStdString());
   }
 
   // TODO check for homo, lumo, or somo energies
-  // m_propertiesCache.setValue("homoEnergy", energy);
-  // m_propertiesCache.setValue("lumoEnergy", energy);
+  const auto* basis = m_molecule->basisSet();
+  const GaussianSet* gaussianSet = dynamic_cast<const GaussianSet*>(basis);
+  if (gaussianSet != nullptr && gaussianSet->scfType() == Core::Rhf) {
+    unsigned int homo = gaussianSet->homo();
+    unsigned int lumo = gaussianSet->lumo();
+    const auto moEnergies = gaussianSet->moEnergy();
+    if (moEnergies.size() > homo) {
+      m_propertiesCache.setValue("homoEnergy", moEnergies[homo] * AU_TO_EV);
+    }
+    // look for the lumo if there's a degenerate HOMO
+    const double threshold = 0.01 / AU_TO_EV; // 0.01 eV minimal separation
+    while (moEnergies.size() > lumo &&
+           std::abs(moEnergies[lumo] - moEnergies[homo]) < threshold) {
+      lumo += 1;
+    }
+    if (moEnergies.size() > lumo)
+      m_propertiesCache.setValue("lumoEnergy", moEnergies[lumo] * AU_TO_EV);
+  }
   // m_propertiesCache.setValue("somoEnergy", energy);
 
   // ignore potentially duplicate properties

diff --git a/avogadro/qtplugins/openbabel/obfileformat.cpp b/avogadro/qtplugins/openbabel/obfileformat.cpp
@@ -240,7 +240,9 @@ bool OBFileFormat::write(std::ostream& out, const Core::Molecule& molecule)
     }
   }
 
+#ifndef NDEBUG
   qDebug() << " writing to " << m_defaultFormat.c_str();
+#endif
 
   // Generate CML or CJSON to give to OpenBabel
   std::string outputString;

diff --git a/avogadro/quantumio/gaussianfchk.cpp b/avogadro/quantumio/gaussianfchk.cpp
@@ -115,9 +115,9 @@ void GaussianFchk::processLine(std::istream& in)
   } else if (key == "Number of atoms" && list.size() > 1) {
     m_numAtoms = Core::lexicalCast<int>(list[1]);
   } else if (key == "Charge" && list.size() > 1) {
-    m_charge = Core::lexicalCast<signed char>(list[1]);
+    m_charge = Core::lexicalCast<int>(list[1]);
   } else if (key == "Multiplicity" && list.size() > 1) {
-    m_spin = Core::lexicalCast<char>(list[1]);
+    m_spin = Core::lexicalCast<int>(list[1]);
   } else if (key == "Dipole Moment" && list.size() > 2) {
     vector<double> dipole = readArrayD(in, Core::lexicalCast<int>(list[2]));
     m_dipoleMoment = Vector3(dipole[0], dipole[1], dipole[2]);

diff --git a/avogadro/quantumio/gaussianfchk.h b/avogadro/quantumio/gaussianfchk.h
@@ -71,8 +71,8 @@ class AVOGADROQUANTUMIO_EXPORT GaussianFchk : public Io::FileFormat
   int m_electronsBeta;
   int m_normalModes;
   int m_numAtoms;
-  unsigned char m_spin;
-  signed char m_charge;
+  int m_spin;
+  int m_charge;
   unsigned int m_numBasisFunctions;
   std::vector<int> m_aNums;
   std::vector<double> m_aPos;