Initial electrostatics framework

avogadro/CMakeLists.txt

@@ -39,7 +39,8 @@ endfunction()
 
 add_subdirectory(core)
 include_directories(${CMAKE_CURRENT_BINARY_DIR}/core)
-
+add_subdirectory(calc)
+include_directories(${CMAKE_CURRENT_BINARY_DIR}/calc)
 add_subdirectory(io)
 include_directories(${CMAKE_CURRENT_BINARY_DIR}/io)
 add_subdirectory(quantumio)

avogadro/calc/CMakeLists.txt

@@ -0,0 +1,22 @@
+find_package(Eigen3 REQUIRED)
+
+# Add as "system headers" to avoid warnings generated by them with
+# compilers that support that notion.
+include_directories(SYSTEM "${EIGEN3_INCLUDE_DIR}")
+
+set(HEADERS
+  chargemodel.h
+  chargemanager.h
+  defaultmodel.h
+)
+
+set(SOURCES
+  chargemodel.cpp
+  chargemanager.cpp
+  defaultmodel.cpp
+)
+
+avogadro_add_library(AvogadroCalc ${HEADERS} ${SOURCES})
+
+target_link_libraries(AvogadroCalc
+  PUBLIC AvogadroCore)

avogadro/calc/chargemanager.cpp

@@ -0,0 +1,178 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#include "chargemanager.h"
+#include "chargemodel.h"
+#include "defaultmodel.h"
+
+#include <algorithm>
+#include <memory>
+
+using std::unique_ptr;
+
+namespace Avogadro {
+namespace Calc {
+
+ChargeManager& ChargeManager::instance()
+{
+  static ChargeManager instance;
+  return instance;
+}
+
+void ChargeManager::appendError(const std::string& errorMessage)
+{
+  m_error += errorMessage + "\n";
+}
+
+bool ChargeManager::registerModel(ChargeModel* model)
+{
+  return instance().addModel(model);
+}
+
+bool ChargeManager::unregisterModel(const std::string& identifier)
+{
+  return instance().removeModel(identifier);
+}
+
+bool ChargeManager::addModel(ChargeModel* model)
+{
+  if (model == nullptr) {
+    appendError("Supplied model was null.");
+    return false;
+  }
+
+  if (m_identifiers.find(model->identifier()) != m_identifiers.end()) {
+    appendError("Model " + model->identifier() + " already loaded.");
+    return false;
+  }
+
+  // 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;
+
+  return true;
+}
+
+bool ChargeManager::removeModel(const std::string& identifier)
+{
+  auto ids = m_identifiers[identifier];
+  m_identifiers.erase(identifier);
+
+  ChargeModel* model = m_models[ids];
+
+  if (model != nullptr) {
+    m_models[ids] = nullptr;
+    delete model;
+  }
+
+  return true;
+}
+
+ChargeManager::ChargeManager()
+{
+  // add any default models here (EEM maybe?)
+}
+
+ChargeManager::~ChargeManager()
+{
+  // Delete the models that were loaded.
+  for (std::vector<ChargeModel*>::const_iterator it = m_models.begin();
+       it != m_models.end(); ++it) {
+    delete (*it);
+  }
+  m_models.clear();
+}
+
+std::set<std::string> ChargeManager::identifiersForMolecule(
+  const Core::Molecule& molecule) const
+{
+  // start with the types already in the molecule
+  std::set<std::string> identifiers = molecule.partialChargeTypes();
+
+  // check our models for compatibility
+  for (std::vector<ChargeModel*>::const_iterator it = m_models.begin();
+       it != m_models.end(); ++it) {
+
+    // We check that every element in the molecule
+    // is handled by the model
+    auto mask = (*it)->elements() & molecule.elements();
+    if (mask.count() == molecule.elements().count())
+      identifiers.insert((*it)->identifier()); // this one will work
+  }
+
+  return identifiers;
+}
+
+const MatrixX ChargeManager::partialCharges(
+  const std::string& identifier, const Core::Molecule& molecule) const
+{
+  // first check if the type is found in the molecule
+  // (i.e., read from a file not computed dynamically)
+  auto molIdentifiers = molecule.partialChargeTypes();
+
+  if (molIdentifiers.find(identifier) != molIdentifiers.end()) {
+    return molecule.partialCharges(identifier);
+  }
+
+  // otherwise go through our list
+  if (m_identifiers.find(identifier) == m_identifiers.end()) {
+    MatrixX charges(molecule.atomCount(),
+                    1); // we have to return something, so zeros
+    return charges;
+  }
+
+  const auto id = m_identifiers[identifier];
+  const ChargeModel* model = m_models[id];
+  return model->partialCharges(molecule);
+}
+
+double ChargeManager::potential(const std::string& identifier,
+                                const Core::Molecule& molecule,
+                                const Vector3& point) const
+{
+  // If the type is found in the molecule
+  // we'll use the DefaultModel to handle the potential
+  auto molIdentifiers = molecule.partialChargeTypes();
+
+  if (molIdentifiers.find(identifier) != molIdentifiers.end()) {
+    DefaultModel model(identifier); // so it knows which charges to use
+    return model.potential(molecule, point);
+  }
+
+  // otherwise go through our list
+  if (m_identifiers.find(identifier) == m_identifiers.end()) {
+    return 0.0;
+  }
+
+  const auto id = m_identifiers[identifier];
+  const ChargeModel* model = m_models[id];
+  return model->potential(molecule, point);
+}
+
+Core::Array<double> ChargeManager::potentials(
+  const std::string& identifier, const Core::Molecule& molecule,
+  const Core::Array<Vector3>& points) const
+{
+  // As above
+  auto molIdentifiers = molecule.partialChargeTypes();
+
+  if (molIdentifiers.find(identifier) != molIdentifiers.end()) {
+    DefaultModel model(identifier); // so it knows which charges to use
+    return model.potentials(molecule, points);
+  }
+
+  if (m_identifiers.find(identifier) == m_identifiers.end()) {
+    Core::Array<double> potentials(points.size(), 0.0);
+    return potentials;
+  }
+
+  const auto id = m_identifiers[identifier];
+  const ChargeModel* model = m_models[id];
+  return model->potentials(molecule, points);
+}
+
+} // namespace Calc
+} // namespace Avogadro

avogadro/calc/chargemanager.h

@@ -0,0 +1,158 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#ifndef AVOGADRO_CALC_CHARGEMANAGER_H
+#define AVOGADRO_CALC_CHARGEMANAGER_H
+
+#include "avogadrocalcexport.h"
+
+#include <avogadro/core/array.h>
+#include <avogadro/core/matrix.h>
+#include <avogadro/core/vector.h>
+
+#include <map>
+#include <set>
+#include <string>
+#include <vector>
+
+namespace Avogadro {
+namespace Core {
+class Molecule;
+}
+namespace Calc {
+
+class ChargeModel;
+
+/**
+ * @class ChargeManager chargemanager.h
+ * <avogadro/calc/chargemanager.h>
+ * @brief Class to manage registration, searching and creation of partial charge
+ * models.
+ * @author Geoffrey R. Hutchison
+ *
+ * The charge manager is a singleton class that handles the runtime
+ * registration, search, creation and eventual destruction of electrostatics
+ * models. It can be used to gain a listing of available models, register new
+ * models, etc.
+ *
+ * All electrostatics can take place independent of this code, but for automated
+ * registration and look up, this is the preferred API. It is possible to use
+ * the convenience API without ever dealing directly with a model class.
+ */
+
+class AVOGADROCALC_EXPORT ChargeManager
+{
+public:
+  /**
+   * Get the singleton instance of the charge manager. This instance should
+   * not be deleted.
+   */
+  static ChargeManager& instance();
+
+  /**
+   * @brief Register a new charge model with the manager.
+   * @param format An instance of the format to manage, the manager assumes
+   * ownership of the object passed in.
+   * @return True on success, false on failure.
+   */
+  static bool registerModel(ChargeModel* model);
+
+  /**
+   * @brief Unregister a charge model from the manager.
+   * @param identifier The identifier for the model to remove.
+   * @return True on success, false on failure.
+   */
+  static bool unregisterModel(const std::string& identifier);
+
+  /**
+   * Add the supplied @p model to the manager, registering its ID and other
+   * relevant data for later lookup. The manager assumes ownership of the
+   * supplied object.
+   * @return True on success, false on failure.
+   */
+  bool addModel(ChargeModel* model);
+
+  /**
+   * Remove the model with the identifier @a identifier from the manager.
+   * @return True on success, false on failure.
+   */
+  bool removeModel(const std::string& identifier);
+
+  /**
+   * New instance of the model for the specified @p identifier. Ownership
+   * is passed to the caller.
+   * @param identifier The unique identifier of the format.
+   * @return Instance of the format, nullptr if not found. Ownership passes to
+   * the
+   * caller.
+   */
+  ChargeModel* newModelFromIdentifier(const std::string& identifier) const;
+
+  /**
+   * Get a list of all loaded identifiers
+   */
+  std::set<std::string> identifiers() const;
+
+  /**
+   * @brief Get a list of models that work for this molecule.
+   *
+   * Includes partial charge types in the molecule itself (e.g., from a file)
+   * This is probably the method you want to get a list for a user
+   */
+  std::set<std::string> identifiersForMolecule(
+    const Core::Molecule& molecule) const;
+
+  /**
+   * Note that some models do not have well-defined atomic partial charges
+   * @return atomic partial charges for the molecule, or 0.0 if undefined
+   */
+  const MatrixX partialCharges(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
+   */
+  double potential(const std::string& identifier, const Core::Molecule& mol,
+                   const Vector3& point) const;
+
+  /**
+   * @return the potentials at the point for the molecule, or an array of 0.0 if
+   * the model is not available
+   */
+  Core::Array<double> potentials(const std::string& identifier,
+                                 const Core::Molecule& mol,
+                                 const Core::Array<Vector3>& points) const;
+
+  /**
+   * Get any errors that have been logged when loading models.
+   */
+  std::string error() const;
+
+private:
+  typedef std::map<std::string, size_t> ChargeIdMap;
+
+  ChargeManager();
+  ~ChargeManager();
+
+  ChargeManager(const ChargeManager&);            // Not implemented.
+  ChargeManager& operator=(const ChargeManager&); // Not implemented.
+
+  /**
+   * @brief Append warnings/errors to the error message string.
+   * @param errorMessage The error message to append.
+   */
+  void appendError(const std::string& errorMessage);
+
+  std::vector<ChargeModel*> m_models;
+  mutable ChargeIdMap m_identifiers;
+
+  std::string m_error;
+};
+
+} // namespace Calc
+} // namespace Avogadro
+
+#endif // AVOGADRO_CALC_CHARGEMANAGER_H