Calculate and render dipole moments (#1801)

* Add methods to calculate dipole moments from charge models

Default is to use atomic point charges and convert to Debye

* Render dipole moments - works for organic molecules using MMFF94

* Add support for reading dipole moments from fchk / orca

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

---------

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

avogadro/calc/chargemanager.cpp

@@ -9,6 +9,15 @@
 
 namespace Avogadro::Calc {
 
+// Helper function to convert a string to lowercase
+// to register all lower-case identifiers
+std::string toLower(const std::string& str)
+{
+  std::string result = str;
+  std::transform(result.begin(), result.end(), result.begin(), ::tolower);
+  return result;
+}
+
 ChargeManager& ChargeManager::instance()
 {
   static ChargeManager instance;
@@ -45,17 +54,20 @@ bool ChargeManager::addModel(ChargeModel* model)
   // If we got here then the format is unique enough to be added.
   size_t index = m_models.size();
   m_models.push_back(model);
-  m_identifiers[model->identifier()] = index;
-  m_identifierToName[model->identifier()] = model->name();
+  std::string lowerId = toLower(model->identifier());
+  m_identifiers[lowerId] = index;
+  m_identifierToName[lowerId] = model->name();
 
   return true;
 }
 
 bool ChargeManager::removeModel(const std::string& identifier)
 {
-  auto ids = m_identifiers[identifier];
-  m_identifiers.erase(identifier);
-  m_identifierToName.erase(identifier);
+  std::string lowerId = toLower(identifier);
+
+  auto ids = m_identifiers[lowerId];
+  m_identifiers.erase(lowerId);
+  m_identifierToName.erase(lowerId);
 
   ChargeModel* model = m_models[ids];
 
@@ -69,7 +81,9 @@ bool ChargeManager::removeModel(const std::string& identifier)
 
 std::string ChargeManager::nameForModel(const std::string& identifier) const
 {
-  auto it = m_identifierToName.find(identifier);
+  std::string lowerId = toLower(identifier);
+
+  auto it = m_identifierToName.find(lowerId);
   if (it == m_identifierToName.end()) {
     return identifier;
   }
@@ -114,23 +128,47 @@ MatrixX ChargeManager::partialCharges(const std::string& identifier,
   // first check if the type is found in the molecule
   // (i.e., read from a file not computed dynamically)
   auto molIdentifiers = molecule.partialChargeTypes();
+  std::string lowerId = toLower(identifier);
 
-  if (molIdentifiers.find(identifier) != molIdentifiers.end()) {
-    return molecule.partialCharges(identifier);
+  if (molIdentifiers.find(lowerId) != molIdentifiers.end()) {
+    return molecule.partialCharges(lowerId);
   }
 
   // otherwise go through our list
-  if (m_identifiers.find(identifier) == m_identifiers.end()) {
+  if (m_identifiers.find(lowerId) == m_identifiers.end()) {
     MatrixX charges(molecule.atomCount(),
                     1); // we have to return something, so zeros
     return charges;
   }
 
-  const auto id = m_identifiers[identifier];
+  const auto id = m_identifiers[lowerId];
   const ChargeModel* model = m_models[id];
   return model->partialCharges(molecule);
 }
 
+Vector3 ChargeManager::dipoleMoment(const std::string& identifier,
+                                    const Core::Molecule& molecule) const
+{
+  // If the type is found in the molecule
+  // we'll use the DefaultModel to handle the dipole moment
+  auto molIdentifiers = molecule.partialChargeTypes();
+  std::string lowerId = toLower(identifier);
+
+  if (molIdentifiers.find(lowerId) != molIdentifiers.end()) {
+    DefaultModel model(lowerId); // so it knows which charges to use
+    return model.dipoleMoment(molecule);
+  }
+
+  // otherwise go through our list
+  if (m_identifiers.find(lowerId) == m_identifiers.end()) {
+    return Vector3(0.0, 0.0, 0.0);
+  }
+
+  const auto id = m_identifiers[lowerId];
+  const ChargeModel* model = m_models[id];
+  return model->dipoleMoment(molecule);
+}
+
 double ChargeManager::potential(const std::string& identifier,
                                 Core::Molecule& molecule,
                                 const Vector3& point) const

avogadro/calc/chargemanager.h

@@ -97,12 +97,29 @@ class AVOGADROCALC_EXPORT ChargeManager
   std::string nameForModel(const std::string& identifier) const;
 
   /**
+   * This method will calculate atomic partial charges and save them,
+   * if the model is available (i.e., the molecule will change)
    * Note that some models do not have well-defined atomic partial charges
    * @return atomic partial charges for the molecule, or 0.0 if undefined
    */
   MatrixX partialCharges(const std::string& identifier,
                          Core::Molecule& mol) const;
 
+  /**
+   * Calculate the atomic partial charges and leave the molecule unchanged.
+   * Note that some models do not have well-defined atomic partial charges
+   * @return the dipole moment for the molecule, or 0.0 if undefined
+   */
+  MatrixX partialCharges(const std::string& identifier,
+                         const Core::Molecule& mol) const;
+
+  /**
+   * @return the dipole moment for the molecule, or 0.0 if the model is not
+   * available
+   */
+  Vector3 dipoleMoment(const std::string& identifier,
+                       const Core::Molecule& mol) const;
+
   /**
    * @return the potential at the point for the molecule, or 0.0 if the model is
    * not available

avogadro/calc/chargemodel.cpp

@@ -8,6 +8,8 @@
 #include <avogadro/core/array.h>
 #include <avogadro/core/molecule.h>
 
+#include <iostream>
+
 namespace Avogadro {
 
 using Core::Array;
@@ -19,6 +21,27 @@ namespace Calc {
 constexpr double M_PI = 3.14159265358979323846;
 #endif
 
+Vector3 ChargeModel::dipoleMoment(const Molecule& mol) const
+{
+  // default is to get the set of partial atomic charges
+  // (some models might do something more sophisticated)
+  const MatrixX charges = partialCharges(mol);
+  // also get the positions of the atoms
+  const Array<Vector3> positions = mol.atomPositions3d();
+
+  Vector3 dipole(0.0, 0.0, 0.0);
+  if (charges.rows() != positions.size())
+    std::cout << "Error: charges " << charges.rows() << " != positions "
+              << positions.size() << std::endl;
+
+  for (unsigned int i = 0; i < charges.size(); ++i)
+    dipole += charges(i, 0) * positions[i];
+
+  dipole *= 4.80320471257; // convert to Debye from electron-Angstrom
+
+  return dipole;
+}
+
 double ChargeModel::potential(Molecule& mol, const Vector3& point) const
 {
   // default is to get the set of partial atomic charges
@@ -30,6 +53,8 @@ double ChargeModel::potential(Molecule& mol, const Vector3& point) const
   //  calculate the atoms within a cutoff distance
   //  and sum the potentials
 
+  // note this is usually multithreaded by the caller
+  // but more efficient methods can be implemented
   double potential = 0.0;
   for (unsigned int i = 0; i < charges.size(); ++i) {
     double distance = (positions[i] - point).norm();

avogadro/calc/chargemodel.h

@@ -79,6 +79,17 @@ class AVOGADROCALC_EXPORT ChargeModel
 
   virtual MatrixX partialCharges(Core::Molecule& mol) const = 0;
 
+  virtual MatrixX partialCharges(const Core::Molecule& mol) const = 0;
+
+  /**
+   * @brief Calculate the dipole moment of the molecule.
+   *
+   * Defaults to using the partial charges and atomic positions
+   * to calculate the net dipole moment.
+   * @return The dipole moment vector of the molecule
+   */
+  virtual Vector3 dipoleMoment(const Core::Molecule& mol) const;
+
   /**
    * @brief Calculate the electrostatic potential at a particular point in
    * space.

avogadro/calc/defaultmodel.cpp

@@ -24,4 +24,9 @@ MatrixX DefaultModel::partialCharges(Core::Molecule& mol) const
   return mol.partialCharges(m_identifier);
 }
 
+MatrixX DefaultModel::partialCharges(const Core::Molecule& mol) const
+{
+  return mol.partialCharges(m_identifier);
+}
+
 } // namespace Avogadro::Calc

avogadro/calc/defaultmodel.h

@@ -76,6 +76,8 @@ class AVOGADROCALC_EXPORT DefaultModel : public ChargeModel
    */
   virtual MatrixX partialCharges(Core::Molecule& mol) const override;
 
+  virtual MatrixX partialCharges(const Core::Molecule& mol) const override;
+
 protected:
   std::string m_identifier;
   Core::Molecule::ElementMask m_elements;

avogadro/core/variant-inline.h

@@ -8,12 +8,19 @@
 
 #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)
 {
@@ -34,12 +41,28 @@ inline Variant::Variant(const MatrixXf& v) : m_type(Matrix)
   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,11 +82,31 @@ 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)
 {
   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;
@@ -185,6 +228,28 @@ inline bool Variant::setValue(MatrixX matrix)
   return true;
 }
 
+template <>
+inline bool Variant::setValue(Vector3 vector)
+{
+  clear();
+
+  m_type = Vector;
+  m_value.vector = new Vector3(vector);
+
+  return true;
+}
+
+template <>
+inline bool Variant::setValue(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 +396,24 @@ inline const MatrixX& Variant::value() const
   return nullMatrix;
 }
 
+template <>
+inline Vector3 Variant::value() const
+{
+  if (m_type == Vector)
+    return *m_value.vector;
+
+  return Vector3();
+}
+
+template <>
+inline const Vector3& Variant::value() const
+{
+  if (m_type == Vector)
+    return *m_value.vector;
+
+  return Vector3::Zero();
+}
+
 inline void Variant::clear()
 {
   if (m_type == String) {
@@ -424,6 +507,11 @@ inline const MatrixX& Variant::toMatrixRef() const
   return value<const MatrixX&>();
 }
 
+inline Vector3 Variant::toVector3() const
+{
+  return value<Vector3>();
+}
+
 // --- Operators ----------------------------------------------------------- //
 inline Variant& Variant::operator=(const Variant& variant)
 {