Add support for reading dipole moments from fchk / orca

Signed-off-by: Geoff Hutchison <[email protected]>

avogadro/core/variant-inline.h

@@ -8,14 +8,21 @@
 
 #include "variant.h"
 
+#include <iostream>
 #include <sstream>
 
 namespace Avogadro::Core {
 
 inline Variant::Variant() : m_type(Null) {}
 
+inline Variant::Variant(double x, double y, double z) : m_type(Vector)
+{
+  Vector3* v = new Vector3(x, y, z);
+  m_value.vector = v;
+}
+
 template <typename T>
-inline Variant::Variant(T v) : m_type(Null)
+inline Variant::Variant(const T v) : m_type(Null)
 {
   setValue(v);
 }
@@ -34,12 +41,36 @@ inline Variant::Variant(const MatrixXf& v) : m_type(Matrix)
   m_value.matrix = m;
 }
 
+template <>
+inline Variant::Variant(const MatrixX& v) : m_type(Matrix)
+{
+  MatrixX* m = new MatrixX(v.rows(), v.cols());
+  *m = v;
+  m_value.matrix = m;
+}
+
+template <>
+inline Variant::Variant(const Vector3& v) : m_type(Vector)
+{
+  Vector3* _v = new Vector3(v);
+  m_value.vector = _v;
+}
+
+template <>
+inline Variant::Variant(const Vector3f& v) : m_type(Vector)
+{
+  Vector3* _v = new Vector3(v.x(), v.y(), v.z());
+  m_value.vector = _v;
+}
+
 inline Variant::Variant(const Variant& variant) : m_type(variant.type())
 {
   if (m_type == String)
     m_value.string = new std::string(variant.toString());
   else if (m_type == Matrix)
     m_value.matrix = new MatrixX(*variant.m_value.matrix);
+  else if (m_type == Vector)
+    m_value.vector = new Vector3(*variant.m_value.vector);
   else if (m_type != Null)
     m_value = variant.m_value;
 }
@@ -59,18 +90,38 @@ inline bool Variant::isNull() const
   return m_type == Null;
 }
 
+inline bool Variant::setValue(double x, double y, double z)
+{
+  clear();
+
+  m_type = Vector;
+  m_value.vector = new Vector3(x, y, z);
+
+  return true;
+}
+
 template <typename T>
-inline bool Variant::setValue(T v)
+inline bool Variant::setValue(const T v)
 {
   AVO_UNUSED(v);
 
+#ifndef NDEBUG
+#if defined(_MSC_VER)
+  std::cerr << " Variant::setValue() not implemented for " << __FUNCSIG__
+            << std::endl;
+#else
+  std::cerr << " Variant::setValue() not implemented for "
+            << __PRETTY_FUNCTION__ << std::endl;
+#endif
+#endif
+
   clear();
 
   return false;
 }
 
 template <>
-inline bool Variant::setValue(bool v)
+inline bool Variant::setValue(const bool v)
 {
   clear();
 
@@ -81,7 +132,7 @@ inline bool Variant::setValue(bool v)
 }
 
 template <>
-inline bool Variant::setValue(char v)
+inline bool Variant::setValue(const char v)
 {
   clear();
 
@@ -92,7 +143,7 @@ inline bool Variant::setValue(char v)
 }
 
 template <>
-inline bool Variant::setValue(short v)
+inline bool Variant::setValue(const short v)
 {
   clear();
 
@@ -103,7 +154,7 @@ inline bool Variant::setValue(short v)
 }
 
 template <>
-inline bool Variant::setValue(int v)
+inline bool Variant::setValue(const int v)
 {
   clear();
 
@@ -114,7 +165,7 @@ inline bool Variant::setValue(int v)
 }
 
 template <>
-inline bool Variant::setValue(long v)
+inline bool Variant::setValue(const long v)
 {
   clear();
 
@@ -125,7 +176,7 @@ inline bool Variant::setValue(long v)
 }
 
 template <>
-inline bool Variant::setValue(float v)
+inline bool Variant::setValue(const float v)
 {
   clear();
 
@@ -136,7 +187,7 @@ inline bool Variant::setValue(float v)
 }
 
 template <>
-inline bool Variant::setValue(double v)
+inline bool Variant::setValue(const double v)
 {
   clear();
 
@@ -147,7 +198,7 @@ inline bool Variant::setValue(double v)
 }
 
 template <>
-inline bool Variant::setValue(std::string string)
+inline bool Variant::setValue(const std::string string)
 {
   clear();
 
@@ -175,7 +226,7 @@ inline bool Variant::setValue(void* pointer)
 }
 
 template <>
-inline bool Variant::setValue(MatrixX matrix)
+inline bool Variant::setValue(const MatrixX& matrix)
 {
   clear();
 
@@ -185,6 +236,40 @@ inline bool Variant::setValue(MatrixX matrix)
   return true;
 }
 
+template <>
+inline bool Variant::setValue(const MatrixXf& matrix)
+{
+  clear();
+
+  m_type = Matrix;
+  m_value.matrix = new MatrixX(matrix.rows(), matrix.cols());
+  *m_value.matrix = matrix.cast<double>();
+
+  return true;
+}
+
+template <>
+inline bool Variant::setValue(const Vector3& vector)
+{
+  clear();
+
+  m_type = Vector;
+  m_value.vector = new Vector3(vector);
+
+  return true;
+}
+
+template <>
+inline bool Variant::setValue(const Vector3f& vector)
+{
+  clear();
+
+  m_type = Vector;
+  m_value.vector = new Vector3(vector.x(), vector.y(), vector.z());
+
+  return true;
+}
+
 template <typename T>
 inline T Variant::value() const
 {
@@ -331,6 +416,15 @@ inline const MatrixX& Variant::value() const
   return nullMatrix;
 }
 
+template <>
+inline Vector3 Variant::value() const
+{
+  if (m_type == Vector)
+    return *m_value.vector;
+
+  return Vector3();
+}
+
 inline void Variant::clear()
 {
   if (m_type == String) {
@@ -421,11 +515,7 @@ inline MatrixX Variant::toMatrix() const
 
 inline Vector3 Variant::toVector3() const
 {
-  MatrixX m = toMatrix();
-  if (m.rows() == 3 && m.cols() == 1)
-    return Vector3(m(0), m(1), m(2));
-  else
-    return Vector3(0.0, 0.0, 0.0);
+  return value<Vector3>();
 }
 
 inline const MatrixX& Variant::toMatrixRef() const

avogadro/core/variant.h

@@ -38,6 +38,7 @@ class AVOGADROCORE_EXPORT Variant
     Double,
     Pointer,
     String,
+    Vector,
     Matrix
   };
 
@@ -46,7 +47,10 @@ class AVOGADROCORE_EXPORT Variant
 
   /** Creates a variant to store @p value. */
   template <typename T>
-  Variant(T value);
+  Variant(const T value);
+
+  /** Creates a variant to store a 3D vector */
+  Variant(double x, double y, double z);
 
   /** Creates a new copy of @p variant. */
   inline Variant(const Variant& variant);
@@ -62,7 +66,10 @@ class AVOGADROCORE_EXPORT Variant
 
   /** Sets the value of the variant to @p value. */
   template <typename T>
-  bool setValue(T value);
+  bool setValue(const T value);
+
+  /** Sets the value of the variant to a 3D vector */
+  bool setValue(double x, double y, double z);
 
   /** @return the value of the variant in the type given by \c T. */
   template <typename T>
@@ -145,6 +152,7 @@ class AVOGADROCORE_EXPORT Variant
     double _double;
     void* pointer;
     std::string* string;
+    Vector3* vector;
     MatrixX* matrix;
   } m_value;
 };

avogadro/io/cjsonformat.cpp

@@ -609,15 +609,22 @@ bool CjsonFormat::deserialize(std::istream& file, Molecule& molecule,
       if (properties.find("totalCharge") != properties.end()) {
         molecule.setData("totalCharge",
                          static_cast<int>(properties["totalCharge"]));
-      }
-      if (properties.find("totalSpinMultiplicity") != properties.end()) {
+      } else if (properties.find("totalSpinMultiplicity") != properties.end()) {
         molecule.setData("totalSpinMultiplicity",
                          static_cast<int>(properties["totalSpinMultiplicity"]));
+      } else if (properties.find("dipoleMoment") != properties.end()) {
+        // read the numeric array
+        json dipole = properties["dipoleMoment"];
+        if (isNumericArray(dipole) && dipole.size() == 3) {
+          Core::Variant dipoleMoment(dipole[0], dipole[1], dipole[2]);
+          molecule.setData("dipoleMoment", dipoleMoment);
+        }
       }
       // iterate through everything else
       for (auto& element : properties.items()) {
         if (element.key() == "totalCharge" ||
-            element.key() == "totalSpinMultiplicity") {
+            element.key() == "totalSpinMultiplicity" ||
+            element.key() == "dipoleMoment") {
           continue;
         }
         if (element.value().type() == json::value_t::array) {

avogadro/qtplugins/dipole/dipole.cpp

@@ -47,9 +47,13 @@ void Dipole::process(const QtGui::Molecule& molecule,
 
   // check if the molecule has the dipole set
   if (!m_customDipole) {
-    if (!molecule.isInteractive()) {
-      m_dipoleVector =
-        Calc::ChargeManager::instance().dipoleMoment(m_type, molecule);
+    if (molecule.hasData("dipoleMoment")) {
+      m_dipoleVector = molecule.data("dipoleMoment").toVector3();
+    } else {
+      if (!molecule.isInteractive()) {
+        m_dipoleVector =
+          Calc::ChargeManager::instance().dipoleMoment(m_type, molecule);
+      }
     }
   }
 

avogadro/qtplugins/molecularproperties/molecularmodel.cpp

@@ -241,6 +241,10 @@ QVariant MolecularModel::data(const QModelIndex& index, int role) const
     return QVariant::fromValue(m_molecule->totalCharge());
   else if (key == " 10totalSpinMultiplicity")
     return QVariant::fromValue(m_molecule->totalSpinMultiplicity());
+  else if (key == "dipoleMoment") {
+    auto dipole = m_molecule->data("dipoleMoment").toVector3();
+    return QString::fromValue(dipole.norm());
+  }
 
   return QString::fromStdString(it->second.toString());
 }
@@ -291,6 +295,8 @@ QVariant MolecularModel::headerData(int section, Qt::Orientation orientation,
       return tr("Net Charge");
     else if (it->first == " 10totalSpinMultiplicity")
       return tr("Net Spin Multiplicity");
+    else if (it->first == "dipoleMoment")
+      return tr("Dipole Moment (Debye)");
     else if (it->first == "homoEnergy")
       return tr("HOMO Energy (eV)", "highest occupied molecular orbital");
     else if (it->first == "lumoEnergy")

avogadro/qtplugins/openbabel/obprocess.cpp

@@ -570,9 +570,11 @@ void OBProcess::executeObabel(const QStringList& options, QObject* receiver,
   }
 
   // Start process
+#ifndef NDEBUG
   qDebug() << "OBProcess::executeObabel: "
               "Running"
            << m_obabelExecutable << options.join(" ");
+#endif
   m_process->start(m_obabelExecutable, options);
   if (!obabelStdin.isNull()) {
     m_process->write(obabelStdin);

avogadro/quantumio/gaussianfchk.cpp

@@ -73,7 +73,10 @@ bool GaussianFchk::read(std::istream& in, Core::Molecule& molecule)
   // set the spin multiplicity
   molecule.setData("totalSpinMultiplicity", m_spin);
   // dipole moment
-  molecule.setData("dipoleMoment", m_dipoleMoment);
+  // TODO: This should be a Vector3d
+  Core::Variant dipole(m_dipoleMoment.x(), m_dipoleMoment.y(),
+                       m_dipoleMoment.z());
+  molecule.setData("dipoleMoment", dipole);
 
   // Do simple bond perception.
   molecule.perceiveBondsSimple();

avogadro/quantumio/orca.cpp

@@ -132,7 +132,10 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
 
   molecule.setData("totalCharge", m_charge);
   molecule.setData("totalSpinMultiplicity", m_spin);
-  molecule.setData("dipoleMoment", m_dipoleMoment);
+  // at the moment, Variant doesn't want to take Vector3d
+  Core::Variant dipole(m_dipoleMoment.x(), m_dipoleMoment.y(),
+                       m_dipoleMoment.z());
+  molecule.setData("dipoleMoment", dipole);
 
   return true;
 }

avogadro/quantumio/orca.h

@@ -62,7 +62,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
   std::vector<int> m_atomNums;
   std::vector<Eigen::Vector3d> m_atomPos;
 
-  Eigen::Vector3d m_dipoleMoment;
+  Vector3 m_dipoleMoment;
 
   std::vector<std::vector<int>> m_bondOrders;