From a621bc39418fe5b15de6d31aa26528f3d0661f16 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 28 May 2022 07:50:18 -0400 Subject: [PATCH 1/8] Add basic support for partial charges, invalidates on modification Signed-off-by: Geoff Hutchison --- avogadro/CMakeLists.txt | 3 +- avogadro/calc/CMakeLists.txt | 20 +++ avogadro/calc/chargemanager.cpp | 125 ++++++++++++++++++ avogadro/calc/chargemanager.h | 124 +++++++++++++++++ avogadro/calc/chargemodel.cpp | 57 ++++++++ avogadro/calc/chargemodel.h | 116 ++++++++++++++++ avogadro/core/molecule.cpp | 39 +++++- avogadro/core/molecule.h | 40 +++--- avogadro/core/variantmap.h | 5 + avogadro/io/cjsonformat.cpp | 15 +++ .../propertytables/propertymodel.cpp | 20 ++- .../qtplugins/propertytables/propertymodel.h | 1 + 12 files changed, 546 insertions(+), 19 deletions(-) create mode 100644 avogadro/calc/CMakeLists.txt create mode 100644 avogadro/calc/chargemanager.cpp create mode 100644 avogadro/calc/chargemanager.h create mode 100644 avogadro/calc/chargemodel.cpp create mode 100644 avogadro/calc/chargemodel.h diff --git a/avogadro/CMakeLists.txt b/avogadro/CMakeLists.txt index b878f5f510..843ef4b6c3 100644 --- a/avogadro/CMakeLists.txt +++ b/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) diff --git a/avogadro/calc/CMakeLists.txt b/avogadro/calc/CMakeLists.txt new file mode 100644 index 0000000000..8c563a586d --- /dev/null +++ b/avogadro/calc/CMakeLists.txt @@ -0,0 +1,20 @@ +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 +) + +set(SOURCES + chargemodel.cpp + chargemanager.cpp +) + +avogadro_add_library(AvogadroCalc ${HEADERS} ${SOURCES}) + +target_link_libraries(AvogadroCalc + PUBLIC AvogadroCore) diff --git a/avogadro/calc/chargemanager.cpp b/avogadro/calc/chargemanager.cpp new file mode 100644 index 0000000000..ae4f55cefa --- /dev/null +++ b/avogadro/calc/chargemanager.cpp @@ -0,0 +1,125 @@ +/****************************************************************************** + 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 +#include + +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) { + appendError("Supplied model was null."); + return false; + } + if (m_identifiers.count(model->identifier()) > 0) { + appendError("Model " + model->identifier() + " already loaded."); + return false; + } + for (std::vector::const_iterator it = m_models.begin(); + it != m_models.end(); ++it) { + if (*it == model) { + appendError("The model object was 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()].push_back(index); + + return true; +} + +namespace { +// Lookup each key from "keys" in "map", and remove "val" from the Map's +// data value (which is a vector of ValueType) +template +void removeFromMap(Map& map, const VectorOfKeys& keys, const ValueType& val) +{ + typedef typename VectorOfKeys::const_iterator KeysIter; + for (KeysIter key = keys.begin(), keyEnd = keys.end(); key != keyEnd; ++key) { + typename Map::iterator mapMatch = map.find(*key); + if (mapMatch == map.end()) + continue; + typename Map::mapped_type& vec = mapMatch->second; + if (vec.size() <= 1) { + map.erase(*key); + } else { + typename Map::mapped_type::iterator newEnd = + std::remove(vec.begin(), vec.end(), val); + vec.resize(newEnd - vec.begin()); + } + } +} +} // namespace + +bool ChargeManager::removeModel(const std::string& identifier) +{ + ChargeIdVector ids = m_identifiers[identifier]; + m_identifiers.erase(identifier); + + if (ids.empty()) + return false; + + for (ChargeIdVector::const_iterator it = ids.begin(), itEnd = ids.end(); + it != itEnd; ++it) { + ChargeModel* model = m_models[*it]; + + if (model == nullptr) + continue; + + m_models[*it] = nullptr; + delete model; + } + + return true; +} + +ChargeManager::ChargeManager() +{ + // add any default models (EEM maybe?) +} + +ChargeManager::~ChargeManager() +{ + // Delete the models that were loaded. + for (std::vector::const_iterator it = m_models.begin(); + it != m_models.end(); ++it) { + delete (*it); + } + m_models.clear(); +} + +} // namespace Calc +} // namespace Avogadro diff --git a/avogadro/calc/chargemanager.h b/avogadro/calc/chargemanager.h new file mode 100644 index 0000000000..e2021cb6b4 --- /dev/null +++ b/avogadro/calc/chargemanager.h @@ -0,0 +1,124 @@ +/****************************************************************************** + 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 +#include +#include + +namespace Avogadro { +namespace Core { +class Molecule; +} +namespace Calc { + +class ChargeModel; + +/** + * @class ChargeManager 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::vector identifiers() const; + + /** + * Get any errors that have been logged when loading models. + */ + std::string error() const; + +private: + typedef std::vector ChargeIdVector; + typedef std::map 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 m_models; + + ChargeIdMap m_identifiers; + + std::string m_error; +}; + +} // namespace Calc +} // namespace Avogadro + +#endif // AVOGADRO_CALC_CHARGEMANAGER_H diff --git a/avogadro/calc/chargemodel.cpp b/avogadro/calc/chargemodel.cpp new file mode 100644 index 0000000000..5a6e61d84a --- /dev/null +++ b/avogadro/calc/chargemodel.cpp @@ -0,0 +1,57 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "chargemodel.h" + +#include +#include + +namespace Avogadro { + +using Core::Array; +using Core::Molecule; + +namespace Calc { + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +ChargeModel::ChargeModel() : m_dielectric(1.0) {} + +ChargeModel::~ChargeModel() {} + +double ChargeModel::potential(const Molecule& mol, const Vector3& point) const +{ + // default is to get the set of partial atomic charges + const Array charges = partialCharges(mol); + // also get the positions of the atoms + const Array positions = mol.atomPositions3d(); + + // @todo: this is naive and inefficient + // calculate the atoms within a cutoff distance + // and sum the potentials + + double potential = 0.0; + for (unsigned i = 0; i < charges.size(); ++i) { + double distance = (positions[i] - point).norm(); + if (distance > 0.01) { + // drop small distances to avoid overflow + potential += charges[i] / distance; + } + } + + return potential / m_dielectric; +} + +void ChargeModel::appendError(const std::string& errorString, bool newLine) +{ + m_error += errorString; + if (newLine) + m_error += "\n"; +} + +} // namespace Calc +} // namespace Avogadro diff --git a/avogadro/calc/chargemodel.h b/avogadro/calc/chargemodel.h new file mode 100644 index 0000000000..2c1e0a040e --- /dev/null +++ b/avogadro/calc/chargemodel.h @@ -0,0 +1,116 @@ +/****************************************************************************** + 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_CHARGEMODEL_H +#define AVOGADRO_CALC_CHARGEMODEL_H + +#include "avogadrocalcexport.h" + +#include +#include +#include + +#include +#include + +namespace Avogadro { + +namespace Core { +class Molecule; +} + +namespace Calc { + +/** + * @class ChargeModel chargemodel.h + * @brief General API for charge / electrostatics models + * @author Geoffrey R. Hutchison + * + * This serves as the common base class for electrostatics models. + * Many use atomic point charges, but we also allow for custom models + * for Slater / Gaussian distributions, distributed electrostatics, + * use of quantum mechanics, etc. + * + * Key methods are to determine either atomic partial charges or + * electrostatic potentials at particular points in space. + * + * There is a default implementation for the electrostatic potential + * at points in space, based on the atomic partial charges. If you + * implement a different mechanism, you should override this method. + */ + +class AVOGADROCALC_EXPORT ChargeModel +{ +public: + ChargeModel(); + virtual ~ChargeModel(); + + /** + * Create a new instance of the file format class. Ownership passes to the + * caller. + */ + virtual ChargeModel* newInstance() const = 0; + + /** + * @brief A unique identifier, used to retrieve models programatically. + * EEM2, NPA, etc. A runtime warning will be generated if the identifier + * is not unique. + */ + virtual std::string identifier() const = 0; + + /** + * @brief The user-visibile name of the model (e.g., "Natural Population + * Analysis") + */ + virtual std::string name() const = 0; + + /** + * A longer description of the model, along with any relevant help text for + * users. This can be used for citations, for example. + */ + virtual std::string description() const = 0; + + /** + * Set the dielectric constant for the model. + * @param dielectric constant. + */ + virtual void setDielectric(double dielectric) { m_dielectric = dielectric; }; + + /** + * @return The dielectric constant. + */ + virtual float dielectric() const { return m_dielectric; } + + virtual Core::Array partialCharges( + const Core::Molecule& mol) const = 0; + + /** + * @brief Calculate the electrostatic potential at a particular point in + * space. + * @param mol The molecule to calculate the potential for. + * @param point The point in space to calculate the potential at. + * @return The electrostatic potential at the point. + */ + virtual double potential(const Core::Molecule& mol, + const Vector3& point) const; + +protected: + /** + * @brief Append an error to the error string for the model. + * @param errorString The error to be added. + * @param newLine Add a new line after the error string? + */ + void appendError(const std::string& errorString, bool newLine = true); + +private: + std::string m_error; + + float m_dielectric; +}; + +} // namespace Calc +} // namespace Avogadro + +#endif // AVOGADRO_CALC_CHARGEMODEL_H diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index c4de67b362..74b4136459 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -31,7 +31,7 @@ Molecule::Molecule() {} Molecule::Molecule(const Molecule& other) - : m_data(other.m_data), m_customElementMap(other.m_customElementMap), + : m_data(other.m_data), m_charges(other.m_charges), m_customElementMap(other.m_customElementMap), m_positions2d(other.m_positions2d), m_positions3d(other.m_positions3d), m_label(other.m_label), m_coordinates3d(other.m_coordinates3d), m_timesteps(other.m_timesteps), m_hybridizations(other.m_hybridizations), @@ -73,6 +73,7 @@ Molecule::Molecule(const Molecule& other) Molecule::Molecule(Molecule&& other) noexcept : m_data(std::move(other.m_data)), + m_charges(std::move(other.m_charges)), m_customElementMap(std::move(other.m_customElementMap)), m_positions2d(std::move(other.m_positions2d)), m_positions3d(std::move(other.m_positions3d)), @@ -114,6 +115,7 @@ Molecule& Molecule::operator=(const Molecule& other) { if (this != &other) { m_data = other.m_data; + m_charges = other.m_charges; m_customElementMap = other.m_customElementMap; m_positions2d = other.m_positions2d; m_positions3d = other.m_positions3d; @@ -172,6 +174,7 @@ Molecule& Molecule::operator=(Molecule&& other) noexcept { if (this != &other) { m_data = std::move(other.m_data); + m_charges = std::move(other.m_charges); m_customElementMap = std::move(other.m_customElementMap); m_positions2d = std::move(other.m_positions2d); m_positions3d = std::move(other.m_positions3d); @@ -237,6 +240,33 @@ const Layer& Molecule::layer() const return m_layers; } +void Molecule::setPartialCharges(const std::string& type, const MatrixX& value) +{ + if (value.size() != atomCount()) + return; + + m_charges[type] = value; +} + +MatrixX Molecule::partialCharges(const std::string& type) const +{ + auto search = m_charges.find(type); + if (search != m_charges.end()) { + return search->second; // value from the map + } else { + MatrixX charges(atomCount(), 1); + return charges; + } +} + +std::vector Molecule::partialChargeTypes() const +{ + std::vector types; + for (auto& it : m_charges) + types.push_back(it.first); + return types; +} + void Molecule::setData(const std::string& name, const Variant& value) { m_data.setValue(name, value); @@ -332,6 +362,7 @@ Molecule::AtomType Molecule::addAtom(unsigned char number) m_graph.addVertex(); m_atomicNumbers.push_back(number); m_layers.addAtomToActiveLayer(atomCount() - 1); + m_charges.clear(); return AtomType(this, static_cast(atomCount() - 1)); } @@ -390,6 +421,7 @@ bool Molecule::removeAtom(Index index) m_selectedAtoms.pop_back(); } + m_charges.clear(); removeBonds(index); m_atomicNumbers.swapAndPop(index); m_graph.removeVertex(index); @@ -415,6 +447,7 @@ void Molecule::clearAtoms() m_atomicNumbers.clear(); m_bondOrders.clear(); m_graph.clear(); + m_charges.clear(); } Molecule::AtomType Molecule::atom(Index index) const @@ -436,6 +469,8 @@ Molecule::BondType Molecule::addBond(Index atom1, Index atom2, } else { m_bondOrders[index] = order; } + // any existing charges are invalidated + m_charges.clear(); return BondType(this, index); } @@ -465,6 +500,7 @@ bool Molecule::removeBond(Index index) return false; m_graph.removeEdge(index); m_bondOrders.swapAndPop(index); + m_charges.clear(); return true; } @@ -488,6 +524,7 @@ void Molecule::clearBonds() m_bondOrders.clear(); m_graph.removeEdges(); m_graph.setSize(atomCount()); + m_charges.clear(); } Molecule::BondType Molecule::bond(Index index) const diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index 3208863d4a..c25bb4f16a 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -69,24 +69,33 @@ class AVOGADROCORE_EXPORT Molecule /** Sets the data value with @p name to @p value. */ void setData(const std::string& name, const Variant& value); - /** Returns the data value for @p name. */ + /** @return the data value for @p name. */ Variant data(const std::string& name) const; /** - * Returns true if the molecule has data with the given key, false otherwise. + * @return true if the molecule has data with the given key, false otherwise. */ bool hasData(const std::string& name) const; /** Set the molecule's variant data to the entries in map. */ void setDataMap(const VariantMap& map); - /** Return the molecule's variant data. */ + /** @return the molecule's variant data. */ const VariantMap& dataMap() const; /** \overload */ VariantMap& dataMap(); - /** Returns a vector of hybridizations for the atoms in the molecule. */ + /** Sets atomic partial charges with @p type to @p value. */ + void setPartialCharges(const std::string& type, const MatrixX& value); + + /** @return the atomic partial charges of type @p type */ + MatrixX partialCharges(const std::string& type) const; + + /** @return the types of partial charges available stored with this molecule. */ + std::vector partialChargeTypes() const; + + /** @return a vector of hybridizations for the atoms in the molecule. */ Array& hybridizations(); /** \overload */ @@ -115,7 +124,7 @@ class AVOGADROCORE_EXPORT Molecule */ bool setHybridization(Index atomId, AtomHybridization hybridization); - /** Returns a vector of formal charges for the atoms in the molecule. */ + /** @return a vector of formal charges for the atoms in the molecule. */ Array& formalCharges(); /** \overload */ @@ -178,7 +187,7 @@ class AVOGADROCORE_EXPORT Molecule bool setLayer(Index atomId, size_t layer); size_t layer(Index atomId) const; - /** Returns a vector of 2d atom positions for the atoms in the molecule. */ + /** @return a vector of 2d atom positions for the atoms in the molecule. */ const Array& atomPositions2d() const; /** \overload */ @@ -207,7 +216,7 @@ class AVOGADROCORE_EXPORT Molecule */ bool setAtomPosition2d(Index atomId, const Vector2& pos); - /** Returns a vector of 2d atom positions for the atoms in the molecule. */ + /** @return a vector of 3d atom positions for the atoms in the molecule. */ const Array& atomPositions3d() const; /** \overload */ @@ -250,7 +259,7 @@ class AVOGADROCORE_EXPORT Molecule */ bool atomSelected(Index atomId) const; - /** Returns whether the selection is empty or not */ + /** @return whether the selection is empty or not */ bool isSelectionEmpty() const; /** A map of custom element atomic numbers to string identifiers. These ids @@ -341,13 +350,13 @@ class AVOGADROCORE_EXPORT Molecule */ virtual void clearBonds(); - /** Returns the bond at @p index in the molecule. */ + /** @return the bond at @p index in the molecule. */ BondType bond(Index index) const; - /** Returns the bond between atoms @p a and @p b. */ + /** @return the bond between atoms @p a and @p b. */ BondType bond(const AtomType& a, const AtomType& b) const; - /** Returns the bond between atomId1 and atomId2. */ + /** @return the bond between atomId1 and atomId2. */ BondType bond(Index atomId1, Index atomId2) const; /** @@ -394,7 +403,7 @@ class AVOGADROCORE_EXPORT Molecule const std::vector cubes() const { return m_cubes; } /** - * Returns the chemical formula of the molecule. + * @return the chemical formula of the molecule. * @param delimiter Delimiter to insert between tokens, defaults to none. * @param showCountsOver Show atom counts above this (defaults to 1). */ @@ -443,7 +452,7 @@ class AVOGADROCORE_EXPORT Molecule void setBasisSet(BasisSet* basis) { m_basisSet = basis; } /** - * Get the basis set (if present) for the molecule. + * @return the basis set (if present) for the molecule. */ BasisSet* basisSet() { return m_basisSet; } const BasisSet* basisSet() const { return m_basisSet; } @@ -556,7 +565,7 @@ class AVOGADROCORE_EXPORT Molecule * @return The number of atoms with the supplied atomic number. */ Index atomCount(unsigned char atomicNumber) const; - /** Returns the number of bonds in the molecule. */ + /** @return the number of bonds in the molecule. */ inline Index bondCount() const; // getters and setters @@ -676,11 +685,12 @@ class AVOGADROCORE_EXPORT Molecule protected: VariantMap m_data; + std::map m_charges; //!< Sets of atomic partial charges CustomElementMap m_customElementMap; Array m_positions2d; Array m_positions3d; Array m_label; - Array> m_coordinates3d; // Used for conformers/trajectories. + Array> m_coordinates3d; //!< Store conformers/trajectories. Array m_timesteps; Array m_hybridizations; Array m_formalCharges; diff --git a/avogadro/core/variantmap.h b/avogadro/core/variantmap.h index 8aa0bf2d96..03c2e687af 100644 --- a/avogadro/core/variantmap.h +++ b/avogadro/core/variantmap.h @@ -58,6 +58,11 @@ class AVOGADROCORE_EXPORT VariantMap */ bool hasValue(const std::string& name) const; + /** + * Clears the map. + */ + void clear() { m_map.clear(); } + /** Return an iterator pointing to the beginning of the map. */ iterator begin(); diff --git a/avogadro/io/cjsonformat.cpp b/avogadro/io/cjsonformat.cpp index 350ed5358b..abf14cc550 100644 --- a/avogadro/io/cjsonformat.cpp +++ b/avogadro/io/cjsonformat.cpp @@ -208,6 +208,21 @@ bool CjsonFormat::read(std::istream& file, Molecule& molecule) } } + // Partial charges are optional, but if present should be loaded. + json partialCharges = atoms["partialCharges"]; + if (partialCharges.is_object()) { + // keys are types, values are arrays of charges + for (auto& kv : partialCharges.items()) { + MatrixX charges(atomCount, 1); + if (isNumericArray(kv.value()) && kv.value().size() == atomCount) { + for (size_t i = 0; i < kv.value().size(); ++i) { + charges(i, 0) = kv.value()[i]; + } + molecule.setPartialCharges(kv.key(), charges); + } + } + } + // Bonds are optional, but if present should be loaded. json bonds = jsonRoot["bonds"]; if (bonds.is_object() && isNumericArray(bonds["connections"]["index"])) { diff --git a/avogadro/qtplugins/propertytables/propertymodel.cpp b/avogadro/qtplugins/propertytables/propertymodel.cpp index bbdfa2a9b2..671486fe75 100644 --- a/avogadro/qtplugins/propertytables/propertymodel.cpp +++ b/avogadro/qtplugins/propertytables/propertymodel.cpp @@ -30,8 +30,8 @@ using std::vector; using SecondaryStructure = Avogadro::Core::Residue::SecondaryStructure; -// element, valence, formal charge, x, y, z, color TODO: partial charge -const int AtomColumns = 7; +// element, valence, formal charge, partial charge, x, y, z, color +const int AtomColumns = 8; // type, atom 1, atom 2, bond order, length const int BondColumns = 5; // type, atom 1, atom 2, atom 3, angle @@ -116,6 +116,18 @@ int PropertyModel::columnCount(const QModelIndex& parent) const return 0; } +QString partialCharge(Molecule* molecule, int atom) +{ + // TODO: we need to track type and/or calling the charge calculator + float charge = 0.0; + std::vector types = molecule->partialChargeTypes(); + if (types.size() > 0) { + MatrixX charges = molecule->partialCharges(types[0]); + charge = charges(atom, 0); + } + return QString("%L1").arg(charge, 0, 'f', 3); +} + // Qt calls this for multiple "roles" across row / columns in the index // we also combine multiple types into this class, so lots of special cases QVariant PropertyModel::data(const QModelIndex& index, int role) const @@ -200,6 +212,8 @@ QVariant PropertyModel::data(const QModelIndex& index, int role) const return QVariant::fromValue(m_molecule->bonds(row).size()); case AtomDataFormalCharge: return m_molecule->formalCharge(row); + case AtomDataPartialCharge: + return partialCharge(m_molecule, row); case AtomDataX: return QString("%L1").arg(m_molecule->atomPosition3d(row).x(), 0, 'f', 4); @@ -357,6 +371,8 @@ QVariant PropertyModel::headerData(int section, Qt::Orientation orientation, return tr("Valence"); case AtomDataFormalCharge: return tr("Formal Charge"); + case AtomDataPartialCharge: + return tr("Partial Charge"); case AtomDataX: return tr("X (Å)"); case AtomDataY: diff --git a/avogadro/qtplugins/propertytables/propertymodel.h b/avogadro/qtplugins/propertytables/propertymodel.h index 6bd8591639..d776dcf0db 100644 --- a/avogadro/qtplugins/propertytables/propertymodel.h +++ b/avogadro/qtplugins/propertytables/propertymodel.h @@ -86,6 +86,7 @@ public slots: AtomDataElement = 0, AtomDataValence, AtomDataFormalCharge, + AtomDataPartialCharge, AtomDataX, AtomDataY, AtomDataZ, From b559ea357f7c9ba28bd66399e6894a87c41cc04f Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sun, 19 Jun 2022 00:18:49 -0400 Subject: [PATCH 2/8] Add an element mask to allow filtering calc methods by elements For example some methods only do CHONF, or other limited types .. so don't show them to the user. Also helpful for some materials potentials (e.g., silver-only) Signed-off-by: Geoff Hutchison --- avogadro/core/elementdata.h | 2 -- avogadro/core/elements.h | 2 ++ avogadro/core/molecule.cpp | 40 ++++++++++++++++++++++++++++--------- avogadro/core/molecule.h | 7 +++++++ 4 files changed, 40 insertions(+), 11 deletions(-) diff --git a/avogadro/core/elementdata.h b/avogadro/core/elementdata.h index a91c27dad5..e55c32a079 100644 --- a/avogadro/core/elementdata.h +++ b/avogadro/core/elementdata.h @@ -4,8 +4,6 @@ namespace Avogadro { namespace Core { -unsigned char element_count = 119; - const char* element_symbols[] = { "Xx", "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", diff --git a/avogadro/core/elements.h b/avogadro/core/elements.h index 87ca6f3e71..6659b9efe6 100644 --- a/avogadro/core/elements.h +++ b/avogadro/core/elements.h @@ -13,6 +13,8 @@ namespace Avogadro { namespace Core { +const unsigned char element_count = 119; //!< from 0 to 118 + /** * @class Elements elements.h * @brief The Elements class stores basic data about chemical elements. diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index 74b4136459..17e7badb25 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -28,10 +28,13 @@ using std::swap; Molecule::Molecule() : m_basisSet(nullptr), m_unitCell(nullptr), m_layers(LayerManager::getMoleculeLayer(this)) -{} +{ + m_elements.reset(); +} Molecule::Molecule(const Molecule& other) - : m_data(other.m_data), m_charges(other.m_charges), m_customElementMap(other.m_customElementMap), + : m_data(other.m_data), m_charges(other.m_charges), + m_customElementMap(other.m_customElementMap), m_elements(other.m_elements), m_positions2d(other.m_positions2d), m_positions3d(other.m_positions3d), m_label(other.m_label), m_coordinates3d(other.m_coordinates3d), m_timesteps(other.m_timesteps), m_hybridizations(other.m_hybridizations), @@ -44,8 +47,7 @@ Molecule::Molecule(const Molecule& other) m_basisSet(other.m_basisSet ? other.m_basisSet->clone() : nullptr), m_unitCell(other.m_unitCell ? new UnitCell(*other.m_unitCell) : nullptr), m_residues(other.m_residues), m_graph(other.m_graph), - m_bondOrders(other.m_bondOrders), - m_atomicNumbers(other.m_atomicNumbers), + m_bondOrders(other.m_bondOrders), m_atomicNumbers(other.m_atomicNumbers), m_hallNumber(other.m_hallNumber), m_layers(LayerManager::getMoleculeLayer(this)) { @@ -72,10 +74,9 @@ Molecule::Molecule(const Molecule& other) } Molecule::Molecule(Molecule&& other) noexcept - : m_data(std::move(other.m_data)), - m_charges(std::move(other.m_charges)), + : m_data(std::move(other.m_data)), m_charges(std::move(other.m_charges)), m_customElementMap(std::move(other.m_customElementMap)), - m_positions2d(std::move(other.m_positions2d)), + m_elements(other.m_elements), m_positions2d(std::move(other.m_positions2d)), m_positions3d(std::move(other.m_positions3d)), m_label(std::move(other.m_label)), m_coordinates3d(std::move(other.m_coordinates3d)), @@ -117,6 +118,7 @@ Molecule& Molecule::operator=(const Molecule& other) m_data = other.m_data; m_charges = other.m_charges; m_customElementMap = other.m_customElementMap; + m_elements = other.m_elements; m_positions2d = other.m_positions2d; m_positions3d = other.m_positions3d; m_label = other.m_label; @@ -176,6 +178,7 @@ Molecule& Molecule::operator=(Molecule&& other) noexcept m_data = std::move(other.m_data); m_charges = std::move(other.m_charges); m_customElementMap = std::move(other.m_customElementMap); + m_elements = std::move(other.m_elements); m_positions2d = std::move(other.m_positions2d); m_positions3d = std::move(other.m_positions3d); m_label = std::move(other.m_label); @@ -244,12 +247,12 @@ void Molecule::setPartialCharges(const std::string& type, const MatrixX& value) { if (value.size() != atomCount()) return; - + m_charges[type] = value; } MatrixX Molecule::partialCharges(const std::string& type) const -{ +{ auto search = m_charges.find(type); if (search != m_charges.end()) { return search->second; // value from the map @@ -361,6 +364,7 @@ Molecule::AtomType Molecule::addAtom(unsigned char number) { m_graph.addVertex(); m_atomicNumbers.push_back(number); + m_elements.set(number); m_layers.addAtomToActiveLayer(atomCount() - 1); m_charges.clear(); return AtomType(this, static_cast(atomCount() - 1)); @@ -423,6 +427,23 @@ bool Molecule::removeAtom(Index index) m_charges.clear(); removeBonds(index); + + // before we remove, check if there's any other atom of this element + // (e.g., we removed the last oxygen) + auto elementToRemove = m_atomicNumbers[index]; + bool foundAnother = false; + for (Index i = 0; i < atomCount(); ++i) { + if (i == index) + continue; + + if (m_atomicNumbers[index] == elementToRemove) { + foundAnother = true; + break; // we're done + } + } + if (!foundAnother) + m_elements.reset(elementToRemove); + m_atomicNumbers.swapAndPop(index); m_graph.removeVertex(index); @@ -448,6 +469,7 @@ void Molecule::clearAtoms() m_bondOrders.clear(); m_graph.clear(); m_charges.clear(); + m_elements.reset(); } Molecule::AtomType Molecule::atom(Index index) const diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index c25bb4f16a..a40313b112 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -48,6 +48,9 @@ class AVOGADROCORE_EXPORT Molecule /** Type for custom element map. */ typedef std::map CustomElementMap; + /** Type for element masks (e.g., does this molecule contain certain elements) */ + typedef std::bitset ElementMask; + /** Creates a new, empty molecule. */ Molecule(); @@ -275,6 +278,9 @@ class AVOGADROCORE_EXPORT Molecule void setCustomElementMap(const CustomElementMap& map); /** @} */ + /** @return the elements currently in this molecule */ + const ElementMask& elements() const; + /** Adds an atom to the molecule. */ virtual AtomType addAtom(unsigned char atomicNumber); AtomType addAtom(unsigned char atomicNumber, Vector3 position3d); @@ -687,6 +693,7 @@ class AVOGADROCORE_EXPORT Molecule VariantMap m_data; std::map m_charges; //!< Sets of atomic partial charges CustomElementMap m_customElementMap; + ElementMask m_elements; //!< Which elements this molecule contains (e.g., for force fields) Array m_positions2d; Array m_positions3d; Array m_label; From 63a71808825a0c71ccbd20b124c68a3a576f9199 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sun, 19 Jun 2022 13:49:32 -0400 Subject: [PATCH 3/8] Add a default model to retrieve partial charges stored in a molecule Class handles electrostatics from atomic partial charges (e.g., Mulliken from a quantum program, or Gasteiger from Open Babel) Signed-off-by: Geoff Hutchison --- avogadro/calc/CMakeLists.txt | 2 + avogadro/calc/chargemodel.cpp | 17 +++++-- avogadro/calc/chargemodel.h | 29 +++++++---- avogadro/calc/defaultmodel.cpp | 36 ++++++++++++++ avogadro/calc/defaultmodel.h | 88 ++++++++++++++++++++++++++++++++++ avogadro/core/molecule.cpp | 2 +- avogadro/core/molecule.h | 2 +- 7 files changed, 161 insertions(+), 15 deletions(-) create mode 100644 avogadro/calc/defaultmodel.cpp create mode 100644 avogadro/calc/defaultmodel.h diff --git a/avogadro/calc/CMakeLists.txt b/avogadro/calc/CMakeLists.txt index 8c563a586d..0e565e8e34 100644 --- a/avogadro/calc/CMakeLists.txt +++ b/avogadro/calc/CMakeLists.txt @@ -7,11 +7,13 @@ 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}) diff --git a/avogadro/calc/chargemodel.cpp b/avogadro/calc/chargemodel.cpp index 5a6e61d84a..98261751d7 100644 --- a/avogadro/calc/chargemodel.cpp +++ b/avogadro/calc/chargemodel.cpp @@ -26,7 +26,7 @@ ChargeModel::~ChargeModel() {} double ChargeModel::potential(const Molecule& mol, const Vector3& point) const { // default is to get the set of partial atomic charges - const Array charges = partialCharges(mol); + const MatrixX charges = partialCharges(mol); // also get the positions of the atoms const Array positions = mol.atomPositions3d(); @@ -35,17 +35,28 @@ double ChargeModel::potential(const Molecule& mol, const Vector3& point) const // and sum the potentials double potential = 0.0; - for (unsigned i = 0; i < charges.size(); ++i) { + for (unsigned int i = 0; i < charges.size(); ++i) { double distance = (positions[i] - point).norm(); if (distance > 0.01) { // drop small distances to avoid overflow - potential += charges[i] / distance; + potential += charges(i,0) / distance; } } return potential / m_dielectric; } +Array& ChargeModel::potentials(const Core::Molecule& mol, + const Array& points) const +{ + // This is naive and slow, but can be re-implemented by methods + // for batching + Array potentials(points.size(), 0.0); + for(unsigned int i = 0; i < points.size(); ++i) + potentials[i] = potential(mol, points[i]); + return potentials; +} + void ChargeModel::appendError(const std::string& errorString, bool newLine) { m_error += errorString; diff --git a/avogadro/calc/chargemodel.h b/avogadro/calc/chargemodel.h index 2c1e0a040e..088c7f7315 100644 --- a/avogadro/calc/chargemodel.h +++ b/avogadro/calc/chargemodel.h @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -17,10 +18,6 @@ namespace Avogadro { -namespace Core { -class Molecule; -} - namespace Calc { /** @@ -48,7 +45,7 @@ class AVOGADROCALC_EXPORT ChargeModel virtual ~ChargeModel(); /** - * Create a new instance of the file format class. Ownership passes to the + * Create a new instance of the model. Ownership passes to the * caller. */ virtual ChargeModel* newInstance() const = 0; @@ -61,16 +58,16 @@ class AVOGADROCALC_EXPORT ChargeModel virtual std::string identifier() const = 0; /** - * @brief The user-visibile name of the model (e.g., "Natural Population + * @brief A user-visibile name of the model (e.g., "Natural Population * Analysis") */ virtual std::string name() const = 0; /** - * A longer description of the model, along with any relevant help text for - * users. This can be used for citations, for example. + * @brief Indicate if your method only treats a subset of elements + * @return an element mask corresponding to the defined subset */ - virtual std::string description() const = 0; + virtual Core::Molecule::ElementMask elements() const = 0; /** * Set the dielectric constant for the model. @@ -83,7 +80,7 @@ class AVOGADROCALC_EXPORT ChargeModel */ virtual float dielectric() const { return m_dielectric; } - virtual Core::Array partialCharges( + virtual const MatrixX& partialCharges( const Core::Molecule& mol) const = 0; /** @@ -96,6 +93,18 @@ class AVOGADROCALC_EXPORT ChargeModel virtual double potential(const Core::Molecule& mol, const Vector3& point) const; + /** + * @brief Calculate the electrostatic potential at multiple points + * @param mol The molecule to calculate the potential for. + * @param array The points in space to calculate the potential at. + * @return The electrostatic potential at the points in an array. + * + * This method is used for batch calculation and defaults to simply + * calculating each point at a time. Some methods work faster in batches. + */ + virtual Core::Array& potentials(const Core::Molecule& mol, + const Core::Array& points) const; + protected: /** * @brief Append an error to the error string for the model. diff --git a/avogadro/calc/defaultmodel.cpp b/avogadro/calc/defaultmodel.cpp new file mode 100644 index 0000000000..1564ce141e --- /dev/null +++ b/avogadro/calc/defaultmodel.cpp @@ -0,0 +1,36 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "defaultmodel.h" + +#include +#include + +namespace Avogadro { + +using Core::Array; +using Core::Molecule; + +namespace Calc { + +DefaultModel::DefaultModel(const std::string id) + : m_identifier(id), ChargeModel() +{ + // we don't know which elements are in the molecule + // but we can just say all of them are okay + // (because this method should work for any molecule) + m_elements.set(); +} + +DefaultModel::~DefaultModel() {} + +const MatrixX& DefaultModel::partialCharges( + const Core::Molecule& mol) const +{ + return mol.partialCharges(m_identifier); +} + +} // namespace Calc +} // namespace Avogadro diff --git a/avogadro/calc/defaultmodel.h b/avogadro/calc/defaultmodel.h new file mode 100644 index 0000000000..59df25e7f0 --- /dev/null +++ b/avogadro/calc/defaultmodel.h @@ -0,0 +1,88 @@ +/****************************************************************************** + 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_DEFAULTMODEL_H +#define AVOGADRO_CALC_DEFAULTMODEL_H + +#include "avogadrocalcexport.h" + +#include + +namespace Avogadro { + +namespace Core { +class Molecule; +} + +namespace Calc { + +/** + * @class DefaultModel defaultmodel.h + * @brief Default charge model for file-provided atomic charges + * @author Geoffrey R. Hutchison + * + * This is a default model for using atomic partial charges from a + * file (e.g., quantum chemistry packages often provide Mulliken charges) + * + * The class + */ + +class AVOGADROCALC_EXPORT DefaultModel : public ChargeModel +{ +public: + DefaultModel(const std::string identifier = ""); + virtual ~DefaultModel(); + + /** + * Create a new instance of the file format class. Ownership passes to the + * caller. + */ + virtual DefaultModel* newInstance() const override + { + return new DefaultModel; + } + + /** + * @brief A unique identifier defined by the file + */ + virtual std::string identifier() const override { return m_identifier; } + + /** + * @brief Set the identifier + */ + virtual void setIdentifier(const std::string identifier) + { + m_identifier = identifier; + } + + /** + * @brief We don't have any other name beyond the identifier in the file + */ + virtual std::string name() const override { return m_identifier; } + + /** + * @brief This default method is defined for whatever is in a molecule + * @return all elements - technically not true, but we don't have the mol + */ + virtual Core::Molecule::ElementMask elements() const override + { + return (m_elements); + } + + /** + * @brief Retrieve the relevant charges from the molecule for our defined type + */ + virtual const MatrixX& partialCharges( + const Core::Molecule& mol) const override; + +protected: + std::string m_identifier; + Core::Molecule::ElementMask m_elements; +}; + +} // namespace Calc +} // namespace Avogadro + +#endif // AVOGADRO_CALC_CHARGEMODEL_H diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index 17e7badb25..0f5cc9f23c 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -251,7 +251,7 @@ void Molecule::setPartialCharges(const std::string& type, const MatrixX& value) m_charges[type] = value; } -MatrixX Molecule::partialCharges(const std::string& type) const +const MatrixX& Molecule::partialCharges(const std::string& type) const { auto search = m_charges.find(type); if (search != m_charges.end()) { diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index a40313b112..f781fb345e 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -93,7 +93,7 @@ class AVOGADROCORE_EXPORT Molecule void setPartialCharges(const std::string& type, const MatrixX& value); /** @return the atomic partial charges of type @p type */ - MatrixX partialCharges(const std::string& type) const; + const MatrixX& partialCharges(const std::string& type) const; /** @return the types of partial charges available stored with this molecule. */ std::vector partialChargeTypes() const; From ce7348544ae55b87b25d531090b51f0d6e357bf0 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 20 Jun 2022 13:41:40 -0400 Subject: [PATCH 4/8] Allow access through the ChargeManager class Includes element filtering / masking (e.g., this method only supports some elements) Signed-off-by: Geoff Hutchison --- avogadro/calc/chargemanager.cpp | 136 +++++++++++++----- avogadro/calc/chargemanager.h | 44 +++++- avogadro/core/molecule.cpp | 6 +- avogadro/core/molecule.h | 7 +- .../propertytables/propertymodel.cpp | 5 +- 5 files changed, 149 insertions(+), 49 deletions(-) diff --git a/avogadro/calc/chargemanager.cpp b/avogadro/calc/chargemanager.cpp index ae4f55cefa..e9f5011913 100644 --- a/avogadro/calc/chargemanager.cpp +++ b/avogadro/calc/chargemanager.cpp @@ -5,6 +5,7 @@ #include "chargemanager.h" #include "chargemodel.h" +#include "defaultmodel.h" #include #include @@ -37,11 +38,12 @@ bool ChargeManager::unregisterModel(const std::string& identifier) bool ChargeManager::addModel(ChargeModel* model) { - if (!model) { + if (model == nullptr) { appendError("Supplied model was null."); return false; } - if (m_identifiers.count(model->identifier()) > 0) { + + if (m_identifiers.find(model->identifier()) != m_identifiers.end()) { appendError("Model " + model->identifier() + " already loaded."); return false; } @@ -56,50 +58,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()].push_back(index); + m_identifiers[model->identifier()] = index; return true; } -namespace { -// Lookup each key from "keys" in "map", and remove "val" from the Map's -// data value (which is a vector of ValueType) -template -void removeFromMap(Map& map, const VectorOfKeys& keys, const ValueType& val) -{ - typedef typename VectorOfKeys::const_iterator KeysIter; - for (KeysIter key = keys.begin(), keyEnd = keys.end(); key != keyEnd; ++key) { - typename Map::iterator mapMatch = map.find(*key); - if (mapMatch == map.end()) - continue; - typename Map::mapped_type& vec = mapMatch->second; - if (vec.size() <= 1) { - map.erase(*key); - } else { - typename Map::mapped_type::iterator newEnd = - std::remove(vec.begin(), vec.end(), val); - vec.resize(newEnd - vec.begin()); - } - } -} -} // namespace - bool ChargeManager::removeModel(const std::string& identifier) { - ChargeIdVector ids = m_identifiers[identifier]; + auto ids = m_identifiers[identifier]; m_identifiers.erase(identifier); - if (ids.empty()) - return false; - - for (ChargeIdVector::const_iterator it = ids.begin(), itEnd = ids.end(); - it != itEnd; ++it) { - ChargeModel* model = m_models[*it]; - - if (model == nullptr) - continue; + ChargeModel* model = m_models[ids]; - m_models[*it] = nullptr; + if (model != nullptr) { + m_models[ids] = nullptr; delete model; } @@ -108,7 +80,7 @@ bool ChargeManager::removeModel(const std::string& identifier) ChargeManager::ChargeManager() { - // add any default models (EEM maybe?) + // add any default models here (EEM maybe?) } ChargeManager::~ChargeManager() @@ -121,5 +93,93 @@ ChargeManager::~ChargeManager() m_models.clear(); } +std::set ChargeManager::identifiersForMolecule( + const Core::Molecule& molecule) const +{ + // start with the types already in the molecule + std::set identifiers = molecule.partialChargeTypes(); + + // check our models for compatibility + for (std::vector::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& ChargeManager::potentials( + const std::string& identifier, const Core::Molecule& molecule, + const Core::Array& 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 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 diff --git a/avogadro/calc/chargemanager.h b/avogadro/calc/chargemanager.h index e2021cb6b4..f880884008 100644 --- a/avogadro/calc/chargemanager.h +++ b/avogadro/calc/chargemanager.h @@ -8,8 +8,13 @@ #include "avogadrocalcexport.h" +#include +#include +#include + #include #include +#include #include namespace Avogadro { @@ -88,7 +93,38 @@ class AVOGADROCALC_EXPORT ChargeManager /** * Get a list of all loaded identifiers */ - std::vector identifiers() const; + std::set 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 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& potentials(const std::string& identifier, + const Core::Molecule& mol, + const Core::Array& points) const; /** * Get any errors that have been logged when loading models. @@ -96,8 +132,7 @@ class AVOGADROCALC_EXPORT ChargeManager std::string error() const; private: - typedef std::vector ChargeIdVector; - typedef std::map ChargeIdMap; + typedef std::map ChargeIdMap; ChargeManager(); ~ChargeManager(); @@ -112,8 +147,7 @@ class AVOGADROCALC_EXPORT ChargeManager void appendError(const std::string& errorMessage); std::vector m_models; - - ChargeIdMap m_identifiers; + mutable ChargeIdMap m_identifiers; std::string m_error; }; diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index 0f5cc9f23c..614d3e1884 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -262,11 +262,11 @@ const MatrixX& Molecule::partialCharges(const std::string& type) const } } -std::vector Molecule::partialChargeTypes() const +std::set Molecule::partialChargeTypes() const { - std::vector types; + std::set types; for (auto& it : m_charges) - types.push_back(it.first); + types.insert(it.first); return types; } diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index f781fb345e..3d9e772b6b 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -96,7 +96,7 @@ class AVOGADROCORE_EXPORT Molecule const MatrixX& partialCharges(const std::string& type) const; /** @return the types of partial charges available stored with this molecule. */ - std::vector partialChargeTypes() const; + std::set partialChargeTypes() const; /** @return a vector of hybridizations for the atoms in the molecule. */ Array& hybridizations(); @@ -800,6 +800,11 @@ inline bool Molecule::setFormalCharge(Index atomId, signed char charge) return false; } +inline const Molecule::ElementMask& Molecule::elements() const +{ + return m_elements; +} + inline Vector3ub Molecule::color(Index atomId) const { if (atomId >= atomCount()) diff --git a/avogadro/qtplugins/propertytables/propertymodel.cpp b/avogadro/qtplugins/propertytables/propertymodel.cpp index 671486fe75..c5544aade5 100644 --- a/avogadro/qtplugins/propertytables/propertymodel.cpp +++ b/avogadro/qtplugins/propertytables/propertymodel.cpp @@ -120,9 +120,10 @@ QString partialCharge(Molecule* molecule, int atom) { // TODO: we need to track type and/or calling the charge calculator float charge = 0.0; - std::vector types = molecule->partialChargeTypes(); + std::set types = molecule->partialChargeTypes(); if (types.size() > 0) { - MatrixX charges = molecule->partialCharges(types[0]); + auto first = types.cbegin(); + MatrixX charges = molecule->partialCharges((*first)); charge = charges(atom, 0); } return QString("%L1").arg(charge, 0, 'f', 3); From e4bfa1f790fe1957b2707bd4cc46697be28f070c Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 20 Jun 2022 13:45:10 -0400 Subject: [PATCH 5/8] Include bitset Signed-off-by: Geoff Hutchison --- avogadro/core/molecule.h | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index 3d9e772b6b..5016176243 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -16,6 +16,7 @@ #include "variantmap.h" #include "vector.h" +#include #include #include #include @@ -48,7 +49,8 @@ class AVOGADROCORE_EXPORT Molecule /** Type for custom element map. */ typedef std::map CustomElementMap; - /** Type for element masks (e.g., does this molecule contain certain elements) */ + /** Type for element masks (e.g., does this molecule contain certain elements) + */ typedef std::bitset ElementMask; /** Creates a new, empty molecule. */ @@ -95,7 +97,8 @@ class AVOGADROCORE_EXPORT Molecule /** @return the atomic partial charges of type @p type */ const MatrixX& partialCharges(const std::string& type) const; - /** @return the types of partial charges available stored with this molecule. */ + /** @return the types of partial charges available stored with this molecule. + */ std::set partialChargeTypes() const; /** @return a vector of hybridizations for the atoms in the molecule. */ @@ -693,7 +696,8 @@ class AVOGADROCORE_EXPORT Molecule VariantMap m_data; std::map m_charges; //!< Sets of atomic partial charges CustomElementMap m_customElementMap; - ElementMask m_elements; //!< Which elements this molecule contains (e.g., for force fields) + ElementMask m_elements; //!< Which elements this molecule contains (e.g., for + //!< force fields) Array m_positions2d; Array m_positions3d; Array m_label; @@ -723,7 +727,7 @@ class AVOGADROCORE_EXPORT Molecule unsigned short m_hallNumber = 0; private: - mutable Graph m_graph; // A transformation of the molecule to a graph. + mutable Graph m_graph; // A transformation of the molecule to a graph. // edge information Array m_bondOrders; // vertex information From 0a006300bbe31ce34d115ebfa5ea10f5ed122c28 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 20 Jun 2022 14:31:33 -0400 Subject: [PATCH 6/8] Fix return by local reference Signed-off-by: Geoff Hutchison --- avogadro/calc/chargemanager.cpp | 4 ++-- avogadro/calc/chargemanager.h | 12 ++++++------ avogadro/calc/chargemodel.cpp | 2 +- avogadro/calc/chargemodel.h | 4 ++-- avogadro/calc/defaultmodel.cpp | 4 ++-- avogadro/calc/defaultmodel.h | 6 +++--- avogadro/core/molecule.h | 4 ++-- 7 files changed, 18 insertions(+), 18 deletions(-) diff --git a/avogadro/calc/chargemanager.cpp b/avogadro/calc/chargemanager.cpp index e9f5011913..345011cae0 100644 --- a/avogadro/calc/chargemanager.cpp +++ b/avogadro/calc/chargemanager.cpp @@ -113,7 +113,7 @@ std::set ChargeManager::identifiersForMolecule( return identifiers; } -const MatrixX& ChargeManager::partialCharges( +const MatrixX ChargeManager::partialCharges( const std::string& identifier, const Core::Molecule& molecule) const { // first check if the type is found in the molecule @@ -159,7 +159,7 @@ double ChargeManager::potential(const std::string& identifier, return model->potential(molecule, point); } -Core::Array& ChargeManager::potentials( +Core::Array ChargeManager::potentials( const std::string& identifier, const Core::Molecule& molecule, const Core::Array& points) const { diff --git a/avogadro/calc/chargemanager.h b/avogadro/calc/chargemanager.h index f880884008..f98c58b555 100644 --- a/avogadro/calc/chargemanager.h +++ b/avogadro/calc/chargemanager.h @@ -13,8 +13,8 @@ #include #include -#include #include +#include #include namespace Avogadro { @@ -108,8 +108,8 @@ class AVOGADROCALC_EXPORT ChargeManager * 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; + 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 @@ -122,9 +122,9 @@ class AVOGADROCALC_EXPORT ChargeManager * @return the potentials at the point for the molecule, or an array of 0.0 if * the model is not available */ - Core::Array& potentials(const std::string& identifier, - const Core::Molecule& mol, - const Core::Array& points) const; + Core::Array potentials(const std::string& identifier, + const Core::Molecule& mol, + const Core::Array& points) const; /** * Get any errors that have been logged when loading models. diff --git a/avogadro/calc/chargemodel.cpp b/avogadro/calc/chargemodel.cpp index 98261751d7..a1e1e7906f 100644 --- a/avogadro/calc/chargemodel.cpp +++ b/avogadro/calc/chargemodel.cpp @@ -46,7 +46,7 @@ double ChargeModel::potential(const Molecule& mol, const Vector3& point) const return potential / m_dielectric; } -Array& ChargeModel::potentials(const Core::Molecule& mol, +Array ChargeModel::potentials(const Core::Molecule& mol, const Array& points) const { // This is naive and slow, but can be re-implemented by methods diff --git a/avogadro/calc/chargemodel.h b/avogadro/calc/chargemodel.h index 088c7f7315..3c0a5d28e0 100644 --- a/avogadro/calc/chargemodel.h +++ b/avogadro/calc/chargemodel.h @@ -80,7 +80,7 @@ class AVOGADROCALC_EXPORT ChargeModel */ virtual float dielectric() const { return m_dielectric; } - virtual const MatrixX& partialCharges( + virtual const MatrixX partialCharges( const Core::Molecule& mol) const = 0; /** @@ -102,7 +102,7 @@ class AVOGADROCALC_EXPORT ChargeModel * This method is used for batch calculation and defaults to simply * calculating each point at a time. Some methods work faster in batches. */ - virtual Core::Array& potentials(const Core::Molecule& mol, + virtual Core::Array potentials(const Core::Molecule& mol, const Core::Array& points) const; protected: diff --git a/avogadro/calc/defaultmodel.cpp b/avogadro/calc/defaultmodel.cpp index 1564ce141e..d9662c9408 100644 --- a/avogadro/calc/defaultmodel.cpp +++ b/avogadro/calc/defaultmodel.cpp @@ -15,7 +15,7 @@ using Core::Molecule; namespace Calc { -DefaultModel::DefaultModel(const std::string id) +DefaultModel::DefaultModel(const std::string& id) : m_identifier(id), ChargeModel() { // we don't know which elements are in the molecule @@ -26,7 +26,7 @@ DefaultModel::DefaultModel(const std::string id) DefaultModel::~DefaultModel() {} -const MatrixX& DefaultModel::partialCharges( +const MatrixX DefaultModel::partialCharges( const Core::Molecule& mol) const { return mol.partialCharges(m_identifier); diff --git a/avogadro/calc/defaultmodel.h b/avogadro/calc/defaultmodel.h index 59df25e7f0..a423ca2e36 100644 --- a/avogadro/calc/defaultmodel.h +++ b/avogadro/calc/defaultmodel.h @@ -32,7 +32,7 @@ namespace Calc { class AVOGADROCALC_EXPORT DefaultModel : public ChargeModel { public: - DefaultModel(const std::string identifier = ""); + DefaultModel(const std::string& identifier = ""); virtual ~DefaultModel(); /** @@ -52,7 +52,7 @@ class AVOGADROCALC_EXPORT DefaultModel : public ChargeModel /** * @brief Set the identifier */ - virtual void setIdentifier(const std::string identifier) + virtual void setIdentifier(const std::string& identifier) { m_identifier = identifier; } @@ -74,7 +74,7 @@ class AVOGADROCALC_EXPORT DefaultModel : public ChargeModel /** * @brief Retrieve the relevant charges from the molecule for our defined type */ - virtual const MatrixX& partialCharges( + virtual const MatrixX partialCharges( const Core::Molecule& mol) const override; protected: diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index 5016176243..17fc63b9f6 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -282,7 +282,7 @@ class AVOGADROCORE_EXPORT Molecule /** @} */ /** @return the elements currently in this molecule */ - const ElementMask& elements() const; + const ElementMask elements() const; /** Adds an atom to the molecule. */ virtual AtomType addAtom(unsigned char atomicNumber); @@ -804,7 +804,7 @@ inline bool Molecule::setFormalCharge(Index atomId, signed char charge) return false; } -inline const Molecule::ElementMask& Molecule::elements() const +inline const Molecule::ElementMask Molecule::elements() const { return m_elements; } From 47f47e37be3239504f92ab3264d5f4be8c530021 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 20 Jun 2022 15:39:20 -0400 Subject: [PATCH 7/8] Fixed a pass local reference Signed-off-by: Geoff Hutchison --- avogadro/core/molecule.cpp | 2 +- avogadro/core/molecule.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index 614d3e1884..0747dfc510 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -251,7 +251,7 @@ void Molecule::setPartialCharges(const std::string& type, const MatrixX& value) m_charges[type] = value; } -const MatrixX& Molecule::partialCharges(const std::string& type) const +const MatrixX Molecule::partialCharges(const std::string& type) const { auto search = m_charges.find(type); if (search != m_charges.end()) { diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index 17fc63b9f6..50aa2d686b 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -95,7 +95,7 @@ class AVOGADROCORE_EXPORT Molecule void setPartialCharges(const std::string& type, const MatrixX& value); /** @return the atomic partial charges of type @p type */ - const MatrixX& partialCharges(const std::string& type) const; + const MatrixX partialCharges(const std::string& type) const; /** @return the types of partial charges available stored with this molecule. */ From 38f4fbf74e531836f7f09a157b16d0d7d0a979ae Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 20 Jun 2022 19:49:15 -0400 Subject: [PATCH 8/8] Fix suggestions and test failures Signed-off-by: Geoff Hutchison --- avogadro/calc/chargemanager.cpp | 13 +++------ avogadro/core/molecule.cpp | 47 ++++++++++++++++++++++----------- avogadro/core/molecule.h | 2 +- 3 files changed, 35 insertions(+), 27 deletions(-) diff --git a/avogadro/calc/chargemanager.cpp b/avogadro/calc/chargemanager.cpp index 345011cae0..cb2b536368 100644 --- a/avogadro/calc/chargemanager.cpp +++ b/avogadro/calc/chargemanager.cpp @@ -47,13 +47,6 @@ bool ChargeManager::addModel(ChargeModel* model) appendError("Model " + model->identifier() + " already loaded."); return false; } - for (std::vector::const_iterator it = m_models.begin(); - it != m_models.end(); ++it) { - if (*it == model) { - appendError("The model object was already loaded."); - return false; - } - } // If we got here then the format is unique enough to be added. size_t index = m_models.size(); @@ -125,7 +118,7 @@ const MatrixX ChargeManager::partialCharges( } // otherwise go through our list - if (m_identifiers.find(identifier) != m_identifiers.end()) { + if (m_identifiers.find(identifier) == m_identifiers.end()) { MatrixX charges(molecule.atomCount(), 1); // we have to return something, so zeros return charges; @@ -150,7 +143,7 @@ double ChargeManager::potential(const std::string& identifier, } // otherwise go through our list - if (m_identifiers.find(identifier) != m_identifiers.end()) { + if (m_identifiers.find(identifier) == m_identifiers.end()) { return 0.0; } @@ -171,7 +164,7 @@ Core::Array ChargeManager::potentials( return model.potentials(molecule, points); } - if (m_identifiers.find(identifier) != m_identifiers.end()) { + if (m_identifiers.find(identifier) == m_identifiers.end()) { Core::Array potentials(points.size(), 0.0); return potentials; } diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index 0747dfc510..fc583faacd 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -33,7 +33,7 @@ Molecule::Molecule() } Molecule::Molecule(const Molecule& other) - : m_data(other.m_data), m_charges(other.m_charges), + : m_data(other.m_data), m_partialCharges(other.m_partialCharges), m_customElementMap(other.m_customElementMap), m_elements(other.m_elements), m_positions2d(other.m_positions2d), m_positions3d(other.m_positions3d), m_label(other.m_label), m_coordinates3d(other.m_coordinates3d), @@ -74,7 +74,8 @@ Molecule::Molecule(const Molecule& other) } Molecule::Molecule(Molecule&& other) noexcept - : m_data(std::move(other.m_data)), m_charges(std::move(other.m_charges)), + : m_data(std::move(other.m_data)), + m_partialCharges(std::move(other.m_partialCharges)), m_customElementMap(std::move(other.m_customElementMap)), m_elements(other.m_elements), m_positions2d(std::move(other.m_positions2d)), m_positions3d(std::move(other.m_positions3d)), @@ -116,7 +117,7 @@ Molecule& Molecule::operator=(const Molecule& other) { if (this != &other) { m_data = other.m_data; - m_charges = other.m_charges; + m_partialCharges = other.m_partialCharges; m_customElementMap = other.m_customElementMap; m_elements = other.m_elements; m_positions2d = other.m_positions2d; @@ -176,7 +177,7 @@ Molecule& Molecule::operator=(Molecule&& other) noexcept { if (this != &other) { m_data = std::move(other.m_data); - m_charges = std::move(other.m_charges); + m_partialCharges = std::move(other.m_partialCharges); m_customElementMap = std::move(other.m_customElementMap); m_elements = std::move(other.m_elements); m_positions2d = std::move(other.m_positions2d); @@ -248,13 +249,13 @@ void Molecule::setPartialCharges(const std::string& type, const MatrixX& value) if (value.size() != atomCount()) return; - m_charges[type] = value; + m_partialCharges[type] = value; } const MatrixX Molecule::partialCharges(const std::string& type) const { - auto search = m_charges.find(type); - if (search != m_charges.end()) { + auto search = m_partialCharges.find(type); + if (search != m_partialCharges.end()) { return search->second; // value from the map } else { MatrixX charges(atomCount(), 1); @@ -265,7 +266,7 @@ const MatrixX Molecule::partialCharges(const std::string& type) const std::set Molecule::partialChargeTypes() const { std::set types; - for (auto& it : m_charges) + for (auto& it : m_partialCharges) types.insert(it.first); return types; } @@ -364,9 +365,14 @@ Molecule::AtomType Molecule::addAtom(unsigned char number) { m_graph.addVertex(); m_atomicNumbers.push_back(number); - m_elements.set(number); + // we're not going to easily handle custom elements + if (number <= element_count) + m_elements.set(number); + else + m_elements.set(element_count - 1); // custom element + m_layers.addAtomToActiveLayer(atomCount() - 1); - m_charges.clear(); + m_partialCharges.clear(); return AtomType(this, static_cast(atomCount() - 1)); } @@ -425,7 +431,7 @@ bool Molecule::removeAtom(Index index) m_selectedAtoms.pop_back(); } - m_charges.clear(); + m_partialCharges.clear(); removeBonds(index); // before we remove, check if there's any other atom of this element @@ -468,7 +474,7 @@ void Molecule::clearAtoms() m_atomicNumbers.clear(); m_bondOrders.clear(); m_graph.clear(); - m_charges.clear(); + m_partialCharges.clear(); m_elements.reset(); } @@ -492,7 +498,7 @@ Molecule::BondType Molecule::addBond(Index atom1, Index atom2, m_bondOrders[index] = order; } // any existing charges are invalidated - m_charges.clear(); + m_partialCharges.clear(); return BondType(this, index); } @@ -522,7 +528,7 @@ bool Molecule::removeBond(Index index) return false; m_graph.removeEdge(index); m_bondOrders.swapAndPop(index); - m_charges.clear(); + m_partialCharges.clear(); return true; } @@ -546,7 +552,7 @@ void Molecule::clearBonds() m_bondOrders.clear(); m_graph.removeEdges(); m_graph.setSize(atomCount()); - m_charges.clear(); + m_partialCharges.clear(); } Molecule::BondType Molecule::bond(Index index) const @@ -1035,10 +1041,13 @@ bool Molecule::setAtomicNumbers(const Core::Array& nums) if (nums.size() == atomCount()) { m_atomicNumbers = nums; + // update element mask + m_elements.reset(); // update colors too if (nums.size() == m_colors.size()) { for (Index i = 0; i < nums.size(); ++i) { - m_colors[i] = Vector3ub(Elements::color(atomicNumber(i))); + m_colors[i] = Vector3ub(Elements::color(m_atomicNumbers[i])); + m_elements.set(m_atomicNumbers[i]); } } @@ -1052,6 +1061,12 @@ bool Molecule::setAtomicNumber(Index atomId, unsigned char number) if (atomId < atomCount()) { m_atomicNumbers[atomId] = number; + // recalculate the element mask + m_elements.reset(); + for (Index i = 0; i < m_atomicNumbers.size(); ++i) { + m_elements.set(m_atomicNumbers[i]); + } + // update colors too if (atomId < m_colors.size()) m_colors[atomId] = Vector3ub(Elements::color(number)); diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index 50aa2d686b..39eef29801 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -694,7 +694,7 @@ class AVOGADROCORE_EXPORT Molecule protected: VariantMap m_data; - std::map m_charges; //!< Sets of atomic partial charges + std::map m_partialCharges; //!< Sets of atomic partial charges CustomElementMap m_customElementMap; ElementMask m_elements; //!< Which elements this molecule contains (e.g., for //!< force fields)