From 389321120f8c6f83066f7a4d4e0496b474944b58 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Wed, 20 Nov 2024 13:04:24 -0500 Subject: [PATCH] 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 --------- Signed-off-by: Geoff Hutchison --- avogadro/calc/chargemanager.cpp | 58 +++++- avogadro/calc/chargemanager.h | 17 ++ avogadro/calc/chargemodel.cpp | 25 +++ avogadro/calc/chargemodel.h | 11 ++ avogadro/calc/defaultmodel.cpp | 5 + avogadro/calc/defaultmodel.h | 2 + avogadro/core/variant-inline.h | 88 +++++++++ avogadro/core/variant.h | 50 +++-- avogadro/io/cjsonformat.cpp | 13 +- avogadro/qtgui/layermodel.cpp | 15 +- avogadro/qtgui/meshgenerator.h | 175 +++++++++--------- avogadro/qtgui/molecule.cpp | 26 ++- avogadro/qtgui/molecule.h | 15 +- avogadro/qtplugins/CMakeLists.txt | 1 + avogadro/qtplugins/dipole/CMakeLists.txt | 9 + avogadro/qtplugins/dipole/dipole.cpp | 68 +++++++ avogadro/qtplugins/dipole/dipole.h | 55 ++++++ .../molecularproperties/molecularmodel.cpp | 9 +- avogadro/qtplugins/openbabel/obcharges.cpp | 26 ++- avogadro/qtplugins/openbabel/obcharges.h | 1 + avogadro/qtplugins/openbabel/obprocess.cpp | 2 + .../scriptcharges/scriptchargemodel.cpp | 11 +- .../scriptcharges/scriptchargemodel.h | 1 + avogadro/quantumio/gaussianfchk.cpp | 5 +- avogadro/quantumio/orca.h | 2 +- avogadro/rendering/arrowgeometry.cpp | 22 ++- avogadro/rendering/arrowgeometry.h | 4 + 27 files changed, 559 insertions(+), 157 deletions(-) create mode 100644 avogadro/qtplugins/dipole/CMakeLists.txt create mode 100644 avogadro/qtplugins/dipole/dipole.cpp create mode 100644 avogadro/qtplugins/dipole/dipole.h diff --git a/avogadro/calc/chargemanager.cpp b/avogadro/calc/chargemanager.cpp index 14ed5ef85e..e4a397b869 100644 --- a/avogadro/calc/chargemanager.cpp +++ b/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 diff --git a/avogadro/calc/chargemanager.h b/avogadro/calc/chargemanager.h index 2e7785f7f4..350a785853 100644 --- a/avogadro/calc/chargemanager.h +++ b/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 diff --git a/avogadro/calc/chargemodel.cpp b/avogadro/calc/chargemodel.cpp index 526a2382c7..995aea758c 100644 --- a/avogadro/calc/chargemodel.cpp +++ b/avogadro/calc/chargemodel.cpp @@ -8,6 +8,8 @@ #include #include +#include + 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 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(); diff --git a/avogadro/calc/chargemodel.h b/avogadro/calc/chargemodel.h index f275b2692f..758a0e8ace 100644 --- a/avogadro/calc/chargemodel.h +++ b/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. diff --git a/avogadro/calc/defaultmodel.cpp b/avogadro/calc/defaultmodel.cpp index 313eaf1f91..ff98b29b68 100644 --- a/avogadro/calc/defaultmodel.cpp +++ b/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 diff --git a/avogadro/calc/defaultmodel.h b/avogadro/calc/defaultmodel.h index 3e25d3cf7b..eb1b48292e 100644 --- a/avogadro/calc/defaultmodel.h +++ b/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; diff --git a/avogadro/core/variant-inline.h b/avogadro/core/variant-inline.h index 3ce6945d62..77379df89c 100644 --- a/avogadro/core/variant-inline.h +++ b/avogadro/core/variant-inline.h @@ -8,12 +8,19 @@ #include "variant.h" +#include #include 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 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 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 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(); } +inline Vector3 Variant::toVector3() const +{ + return value(); +} + // --- Operators ----------------------------------------------------------- // inline Variant& Variant::operator=(const Variant& variant) { diff --git a/avogadro/core/variant.h b/avogadro/core/variant.h index 09599981a0..c7b1d09f53 100644 --- a/avogadro/core/variant.h +++ b/avogadro/core/variant.h @@ -10,6 +10,7 @@ #include "avogadrocore.h" #include "matrix.h" +#include "vector.h" #include @@ -37,6 +38,7 @@ class AVOGADROCORE_EXPORT Variant Double, Pointer, String, + Vector, Matrix }; @@ -50,73 +52,82 @@ class AVOGADROCORE_EXPORT Variant /** Creates a new copy of @p variant. */ inline Variant(const Variant& variant); + /** Creates a variant to store a 3D vector */ + Variant(double x, double y, double z); + /** Destroys the variant object. */ inline ~Variant(); - /** Returns variant's type. */ + /** @return variant's type. */ inline Type type() const; - /** Returns \c true if the variant is null. */ + /** @return \c true if the variant is null. */ inline bool isNull() const; /** Sets the value of the variant to @p value. */ template bool setValue(T value); - /** Returns the value of the variant in the type given by \c T. */ + /** Sets the value of the variant to a 3D vector */ + bool setValue(double x, double y, double z); + + /** @return the value of the variant in the type given by \c T. */ template T value() const; /** Clears the variant's data and sets the variant to null. */ inline void clear(); - /** Returns the value of the variant as a \c bool. */ + /** @return the value of the variant as a \c bool. */ inline bool toBool() const; - /** Returns the value of the variant as a \c char. */ + /** @return the value of the variant as a \c char. */ inline char toChar() const; - /** Returns the value of the variant as an \c unsigned \c char. */ + /** @return the value of the variant as an \c unsigned \c char. */ inline unsigned char toUChar() const; - /** Returns the value of the variant as a \c short. */ + /** @return the value of the variant as a \c short. */ inline short toShort() const; - /** Returns the value of the variant as an \c unsigned \c short. */ + /** @return the value of the variant as an \c unsigned \c short. */ inline unsigned short toUShort() const; - /** Returns the value of the variant as an \c int. */ + /** @return the value of the variant as an \c int. */ inline int toInt() const; - /** Returns the value of the variant as an \c unsigned \c int. */ + /** @return the value of the variant as an \c unsigned \c int. */ inline unsigned int toUInt() const; - /** Returns the value of the variant as a \c long. */ + /** @return the value of the variant as a \c long. */ inline long toLong() const; - /** Returns the value of the variant as an \c unsigned \c long. */ + /** @return the value of the variant as an \c unsigned \c long. */ inline unsigned long toULong() const; - /** Returns the value of the variant as a \c float. */ + /** @return the value of the variant as a \c float. */ inline float toFloat() const; - /** Returns the value of the variant as a \c double. */ + /** @return the value of the variant as a \c double. */ inline double toDouble() const; - /** Returns the value of the variant as a \c Real. */ + /** @return the value of the variant as a \c Real. */ inline Real toReal() const; - /** Returns the value of the variant as a pointer. */ + /** @return the value of the variant as a pointer. */ inline void* toPointer() const; - /** Returns the value of the variant as a string. */ + /** @return the value of the variant as a string. */ inline std::string toString() const; - /** Returns the value of the variant as a MatrixX. */ + /** @return the value of the variant as a MatrixX. */ inline MatrixX toMatrix() const; + /** @return the value of the variant as a Vector3 */ + inline Vector3 toVector3() const; + /** - * Returns a reference to the value of the variant as a MatrixX. + * @return a reference to the value of the variant as a MatrixX. * This method will not perform any casting -- if type() is not exactly * MatrixX, the function will fail and return a reference to an empty MatrixX. */ @@ -141,6 +152,7 @@ class AVOGADROCORE_EXPORT Variant double _double; void* pointer; std::string* string; + Vector3* vector; MatrixX* matrix; } m_value; }; diff --git a/avogadro/io/cjsonformat.cpp b/avogadro/io/cjsonformat.cpp index f9257258e5..8bde3f0a88 100644 --- a/avogadro/io/cjsonformat.cpp +++ b/avogadro/io/cjsonformat.cpp @@ -609,15 +609,22 @@ bool CjsonFormat::deserialize(std::istream& file, Molecule& molecule, if (properties.find("totalCharge") != properties.end()) { molecule.setData("totalCharge", static_cast(properties["totalCharge"])); - } - if (properties.find("totalSpinMultiplicity") != properties.end()) { + } else if (properties.find("totalSpinMultiplicity") != properties.end()) { molecule.setData("totalSpinMultiplicity", static_cast(properties["totalSpinMultiplicity"])); + } else if (properties.find("dipoleMoment") != properties.end()) { + // read the numeric array + json dipole = properties["dipoleMoment"]; + if (isNumericArray(dipole) && dipole.size() == 3) { + Core::Variant dipoleMoment(dipole[0], dipole[1], dipole[2]); + molecule.setData("dipoleMoment", dipoleMoment); + } } // iterate through everything else for (auto& element : properties.items()) { if (element.key() == "totalCharge" || - element.key() == "totalSpinMultiplicity") { + element.key() == "totalSpinMultiplicity" || + element.key() == "dipoleMoment") { continue; } if (element.value().type() == json::value_t::array) { diff --git a/avogadro/qtgui/layermodel.cpp b/avogadro/qtgui/layermodel.cpp index 17ca8a13bf..c34c31e61f 100644 --- a/avogadro/qtgui/layermodel.cpp +++ b/avogadro/qtgui/layermodel.cpp @@ -20,7 +20,8 @@ namespace { const int QTTY_COLUMNS = 6; } -LayerModel::LayerModel(QObject* p) : QAbstractItemModel(p), m_item(0) { +LayerModel::LayerModel(QObject* p) : QAbstractItemModel(p), m_item(0) +{ m_plusIcon = QIcon(":/icons/fallback/32x32/plus.png"); m_dotsIcon = QIcon(":/icons/fallback/32x32/dots.png"); m_previewIcon = QIcon(":/icons/fallback/32x32/preview.png"); @@ -53,8 +54,7 @@ Qt::ItemFlags LayerModel::flags(const QModelIndex&) const return Qt::ItemIsEnabled; } -bool LayerModel::setData(const QModelIndex&, const QVariant&, - int) +bool LayerModel::setData(const QModelIndex&, const QVariant&, int) { return false; } @@ -82,7 +82,8 @@ QVariant LayerModel::data(const QModelIndex& idx, int role) const if (idx.column() == ColumnType::Name) { switch (role) { case Qt::DisplayRole: { - return QString(tr("Layer %1")).arg(layer + 1); // count starts at 0 internally + return QString(tr("Layer %1")) + .arg(layer + 1); // count starts at 0 internally } case Qt::ForegroundRole: if (layer == getMoleculeLayer().activeLayer()) @@ -139,6 +140,8 @@ QString LayerModel::getTranslatedName(const std::string& name) const return tr("Close Contacts", "rendering of non-covalent close contacts"); else if (name == "Crystal Lattice") return tr("Crystal Lattice"); + else if (name == "Dipole Moment") + return tr("Dipole Moment"); else if (name == "Force") return tr("Force"); else if (name == "Labels") @@ -155,8 +158,6 @@ QString LayerModel::getTranslatedName(const std::string& name) const return tr("Symmetry Elements"); else if (name == "Van der Waals") return tr("Van der Waals"); - else if (name == "Van der Waals (AO)") - return tr("Van der Waals (AO)", "ambient occlusion"); else if (name == "Wireframe") return tr("Wireframe"); @@ -247,4 +248,4 @@ size_t LayerModel::layerCount() const return LayerManager::layerCount(); } -} // namespace Avogadro +} // namespace Avogadro::QtGui diff --git a/avogadro/qtgui/meshgenerator.h b/avogadro/qtgui/meshgenerator.h index 95a9a70a3a..99dd7298fc 100644 --- a/avogadro/qtgui/meshgenerator.h +++ b/avogadro/qtgui/meshgenerator.h @@ -18,7 +18,7 @@ namespace Avogadro { namespace Core { class Cube; class Mesh; -} +} // namespace Core namespace QtGui { @@ -58,7 +58,8 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread * @return True if the MeshGenerator was successfully initialized. */ MeshGenerator(const Core::Cube* cube, Core::Mesh* mesh, float iso, - int passes = 6, bool reverse = false, QObject* parent = nullptr); + int passes = 6, bool reverse = false, + QObject* parent = nullptr); /** * Destructor. @@ -83,74 +84,69 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread void run() override; /** - * It holds the range and starting offsets of isosurface-intersected - * edges along the x, y, and z axes for each grid cell. - */ - struct gridEdge - { - gridEdge() - : xl(0), - xr(0), - xstart(0), - ystart(0), - zstart(0) - {} - - // trim values - // set on pass 1 - int xl; - int xr; - - // modified on pass 2 - // set on pass 3 - int xstart; - int ystart; - int zstart; - }; - - /** - * Handles duplicate vertices (Not implemented). Placeholder for future functionality. - */ - unsigned long duplicate(const Vector3i& c, const Vector3f& pos); + * It holds the range and starting offsets of isosurface-intersected + * edges along the x, y, and z axes for each grid cell. + */ + struct gridEdge + { + gridEdge() : xl(0), xr(0), xstart(0), ystart(0), zstart(0) {} + + // trim values + // set on pass 1 + int xl; + int xr; + + // modified on pass 2 + // set on pass 3 + int xstart; + int ystart; + int zstart; + }; + /** + * Handles duplicate vertices (Not implemented). Placeholder for future + * functionality. + */ + unsigned long duplicate(const Vector3i& c, const Vector3f& pos); - /** + /** * @name Flying Edges - * Methods to implement the "flying edges" method for isosurface mesh generation. - * Flying edges: A high-performance scalable isocontouring algorithm - * Schroeder; Maynard; Geveci; - * 2015 IEEE 5th Symposium on Large Data Analysis and Visualization (LDAV) + * Methods to implement the "flying edges" method for isosurface mesh + * generation. Flying edges: A high-performance scalable isocontouring + * algorithm Schroeder; Maynard; Geveci; 2015 IEEE 5th Symposium on Large Data + * Analysis and Visualization (LDAV) * [10.1109/LDAV.2015.7348069](https://doi.org/10.1109/LDAV.2015.7348069) * Alternate (non-VTK) implementation at * https://github.com/sandialabs/miniIsosurface/blob/master/flyingEdges/ * @{ + */ /** - * Pass 1 for flying edges. Pass1 detects and records - * where the isosurface intersects each row of grid edges - * along the x-axis. - */ + * Pass 1 for flying edges. Pass1 detects and records + * where the isosurface intersects each row of grid edges + * along the x-axis. + */ void FlyingEdgesAlgorithmPass1(); /** - * Pass2 assigns case identifiers to each grid cell based on - * intersected edges and tallies the number of triangles needed - * for mesh construction. - */ + * Pass2 assigns case identifiers to each grid cell based on + * intersected edges and tallies the number of triangles needed + * for mesh construction. + */ void FlyingEdgesAlgorithmPass2(); /** - * Pass3 computes cumulative offsets for triangles - * and vertices and allocates memory for the mesh structures. - */ + * Pass3 computes cumulative offsets for triangles + * and vertices and allocates memory for the mesh structures. + */ void FlyingEdgesAlgorithmPass3(); /** - * Calculates normals, triangles and vertices. - */ + * Calculates normals, triangles and vertices. + */ void FlyingEdgesAlgorithmPass4(); - /**@}*/ + /**@}*/ /** * @return The Cube being used by the class. @@ -158,18 +154,19 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread const Core::Cube* cube() const { return m_cube; } /** - * Determines the x-range (xl to xr) where isosurface intersections - * occur, optimizing calculations within this range. - */ - inline void calcTrimValues( - int& xl, int& xr, int const& j, int const& k) const; + * Determines the x-range (xl to xr) where isosurface intersections + * occur, optimizing calculations within this range. + */ + inline void calcTrimValues(int& xl, int& xr, int const& j, + int const& k) const; /** - * Indicates which edges intersects the isosurface. - */ - inline unsigned char - calcCubeCase(unsigned char const& ec0, unsigned char const& ec1, - unsigned char const& ec2, unsigned char const& ec3) const; + * Indicates which edges intersects the isosurface. + */ + inline unsigned char calcCubeCase(unsigned char const& ec0, + unsigned char const& ec1, + unsigned char const& ec2, + unsigned char const& ec3) const; /** * @return The Mesh being generated by the class. @@ -199,35 +196,34 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread void progressValueChanged(int); protected: - /** - * isCutEdge checks whether the grid edge at position (i, j, k) is - * intersected by the isosurface based on edge case conditions. - * @return Boolean if it's intersected or not. - */ + * isCutEdge checks whether the grid edge at position (i, j, k) is + * intersected by the isosurface based on edge case conditions. + * @return Boolean if it's intersected or not. + */ bool isCutEdge(int const& i, int const& j, int const& k) const; /** - * It computes the 3D intersection point on a cube edge via interpolation. - */ + * It computes the 3D intersection point on a cube edge via interpolation. + */ inline std::array interpolateOnCube( std::array, 8> const& pts, - std::array const& isovals, - unsigned char const& edge) const; - + std::array const& isovals, unsigned char const& edge) const; + /** - * It linearly interpolates between two 3D points, a and b, using - * the given weight to determine the intermediate position. - */ - inline std::array interpolate( - std::array const& a, - std::array const& b, - float const& weight) const; + * It linearly interpolates between two 3D points, a and b, using + * the given weight to determine the intermediate position. + */ + inline std::array interpolate(std::array const& a, + std::array const& b, + float const& weight) const; /** - * calcCaseEdge determines an edge case code (0–3) based on two boolean edge comparisons. - */ - inline unsigned char calcCaseEdge(bool const& prevEdge, bool const& currEdge) const; + * calcCaseEdge determines an edge case code (0–3) based on two boolean edge + * comparisons. + */ + inline unsigned char calcCaseEdge(bool const& prevEdge, + bool const& currEdge) const; float m_iso; /** The value of the isosurface. */ int m_passes; /** Number of smoothing passes to perform. */ @@ -240,21 +236,24 @@ class AVOGADROQTGUI_EXPORT MeshGenerator : public QThread Core::Array m_normals, m_vertices; std::vector gridEdges; // size (m_dim.y() * m_dim.z()) - std::vector cubeCases; //size ((m_dim.x() - 1) * (m_dim.y() - 1) * (m_dim.z() - 1)) - std::vector triCounter; // size ((m_dim.y() - 1) * (m_dim.z() - 1)) - std::vector edgeCases; // size ((m_dim.x() - 1) * (m_dim.y()) * (m_dim.z())) - Core::Array m_triangles; // triangles of a mesh + std::vector + cubeCases; // size ((m_dim.x() - 1) * (m_dim.y() - 1) * (m_dim.z() - 1)) + std::vector triCounter; // size ((m_dim.y() - 1) * (m_dim.z() - 1)) + std::vector + edgeCases; // size ((m_dim.x() - 1) * (m_dim.y()) * (m_dim.z())) + Core::Array m_triangles; // triangles of a mesh int m_progmin; int m_progmax; /** * These are the lookup tables for flying edges. - * Reference : https://github.com/sandialabs/miniIsosurface/blob/master/flyingEdges/util/MarchingCubesTables.h + * Reference : + * https://github.com/sandialabs/miniIsosurface/blob/master/flyingEdges/util/MarchingCubesTables.h */ - static const unsigned char m_numTris[256]; - static const bool m_isCut[256][12]; + static const unsigned char m_numTris[256]; + static const bool m_isCut[256][12]; static const char m_caseTriangles[256][16]; - static const unsigned char m_edgeVertices[12][2]; + static const unsigned char m_edgeVertices[12][2]; }; } // End namespace QtGui diff --git a/avogadro/qtgui/molecule.cpp b/avogadro/qtgui/molecule.cpp index eb91db2f10..93102d97f7 100644 --- a/avogadro/qtgui/molecule.cpp +++ b/avogadro/qtgui/molecule.cpp @@ -13,17 +13,16 @@ namespace Avogadro::QtGui { using std::swap; Molecule::Molecule(QObject* p) - : QObject(p), Core::Molecule(), - m_undoMolecule(new RWMolecule(*this, this)) + : QObject(p), Core::Molecule(), m_undoMolecule(new RWMolecule(*this, this)) { - m_undoMolecule->setInteractive(true); + m_undoMolecule->setInteractive(false); } Molecule::Molecule(const Molecule& other) : QObject(), Core::Molecule(other), m_undoMolecule(new RWMolecule(*this, this)) { - m_undoMolecule->setInteractive(true); + m_undoMolecule->setInteractive(false); // Now assign the unique ids for (Index i = 0; i < atomCount(); i++) m_atomUniqueIds.push_back(i); @@ -204,8 +203,8 @@ void Molecule::swapAtom(Index a, Index b) Core::Molecule::swapAtom(a, b); } -Molecule::BondType Molecule::addBond(Index a, Index b, - unsigned char order, Index uniqueId) +Molecule::BondType Molecule::addBond(Index a, Index b, unsigned char order, + Index uniqueId) { if (uniqueId >= static_cast(m_bondUniqueIds.size()) || m_bondUniqueIds[uniqueId] != MaxIndex) { @@ -290,6 +289,11 @@ void Molecule::emitChanged(unsigned int change) emit changed(change); } +void Molecule::emitUpdate() const +{ + emit update(); +} + Index Molecule::findAtomUniqueId(Index index) const { for (Index i = 0; i < static_cast(m_atomUniqueIds.size()); ++i) @@ -311,4 +315,12 @@ RWMolecule* Molecule::undoMolecule() return m_undoMolecule; } -} // namespace Avogadro +bool Molecule::isInteractive() const +{ + if (m_undoMolecule == nullptr) + return false; + + return m_undoMolecule->isInteractive(); +} + +} // namespace Avogadro::QtGui diff --git a/avogadro/qtgui/molecule.h b/avogadro/qtgui/molecule.h index 816e77c538..20c2c98d45 100644 --- a/avogadro/qtgui/molecule.h +++ b/avogadro/qtgui/molecule.h @@ -234,6 +234,8 @@ class AVOGADROQTGUI_EXPORT Molecule : public QObject, public Core::Molecule RWMolecule* undoMolecule(); + bool isInteractive() const; + void swapBond(Index a, Index b); void swapAtom(Index a, Index b); @@ -244,6 +246,11 @@ public slots: */ void emitChanged(unsigned int change); + /** + * @brief Request an update through the update() signal + */ + void emitUpdate() const; + signals: /** * @brief Indicates that the molecule has changed. @@ -253,7 +260,13 @@ public slots: * change & Atoms == true then atoms were changed in some way, and if * change & Removed == true then one or more atoms were removed. */ - void changed(unsigned int change); + void changed(unsigned int change) const; + + /** + * @brief Request an update of the molecule. + * (e.g., re-compute properties) + */ + void update() const; private: Core::Array m_atomUniqueIds; diff --git a/avogadro/qtplugins/CMakeLists.txt b/avogadro/qtplugins/CMakeLists.txt index c837fb61df..1926f7f749 100644 --- a/avogadro/qtplugins/CMakeLists.txt +++ b/avogadro/qtplugins/CMakeLists.txt @@ -106,6 +106,7 @@ add_subdirectory(configurepython) add_subdirectory(copypaste) add_subdirectory(crystal) add_subdirectory(customelements) +add_subdirectory(dipole) add_subdirectory(editor) add_subdirectory(fetchpdb) add_subdirectory(focus) diff --git a/avogadro/qtplugins/dipole/CMakeLists.txt b/avogadro/qtplugins/dipole/CMakeLists.txt new file mode 100644 index 0000000000..3dc71f20b0 --- /dev/null +++ b/avogadro/qtplugins/dipole/CMakeLists.txt @@ -0,0 +1,9 @@ +avogadro_plugin(Dipole + "Dipole rendering scheme" + ScenePlugin + dipole.h + Dipole + dipole.cpp + "") + +target_link_libraries(Dipole PRIVATE Avogadro::Calc Avogadro::Rendering) diff --git a/avogadro/qtplugins/dipole/dipole.cpp b/avogadro/qtplugins/dipole/dipole.cpp new file mode 100644 index 0000000000..7fb8014f94 --- /dev/null +++ b/avogadro/qtplugins/dipole/dipole.cpp @@ -0,0 +1,68 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "dipole.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Avogadro::QtPlugins { + +using QtGui::Molecule; +using Rendering::ArrowGeometry; +using Rendering::GeometryNode; +using Rendering::GroupNode; + +Dipole::Dipole(QObject* p) : ScenePlugin(p) +{ + m_layerManager = QtGui::PluginLayerManager(m_name); +} + +Dipole::~Dipole() {} + +void Dipole::process(const QtGui::Molecule& molecule, + Rendering::GroupNode& node) +{ + auto* geometry = new GeometryNode; + node.addChild(geometry); + + auto* arrow = new ArrowGeometry; + arrow->identifier().molecule = &molecule; + + arrow->setColor(Vector3ub(255, 0, 0)); + geometry->addDrawable(arrow); + + Vector3f origin = Vector3f::Zero(); + + // check if the molecule has the dipole set + if (!m_customDipole) { + if (molecule.hasData("dipoleMoment")) { + m_dipoleVector = molecule.data("dipoleMoment").toVector3(); + } else { + if (!molecule.isInteractive()) { + m_dipoleVector = + Calc::ChargeManager::instance().dipoleMoment(m_type, molecule); + } + } + } + + arrow->addSingleArrow(m_dipoleVector.cast(), origin); +} + +QWidget* Dipole::setupWidget() +{ + return nullptr; +} + +} // namespace Avogadro::QtPlugins diff --git a/avogadro/qtplugins/dipole/dipole.h b/avogadro/qtplugins/dipole/dipole.h new file mode 100644 index 0000000000..0abae0feb2 --- /dev/null +++ b/avogadro/qtplugins/dipole/dipole.h @@ -0,0 +1,55 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#ifndef AVOGADRO_QTPLUGINS_DIPOLE_H +#define AVOGADRO_QTPLUGINS_DIPOLE_H + +#include +#include + +namespace Avogadro { +namespace QtPlugins { + +/** + * @brief Render a molecule dipole moment arrow + */ +class Dipole : public QtGui::ScenePlugin +{ + Q_OBJECT + +public: + explicit Dipole(QObject* parent = nullptr); + ~Dipole() override; + + void process(const QtGui::Molecule& molecule, + Rendering::GroupNode& node) override; + + QString name() const override { return tr("Dipole Moment"); } + + QString description() const override + { + return tr("Render the dipole moment of the molecule."); + } + + DefaultBehavior defaultBehavior() const override + { + return DefaultBehavior::False; + } + + QWidget* setupWidget() override; + bool hasSetupWidget() const override { return true; } + +private: + std::string m_name = "Dipole Moment"; + std::string m_type = "MMFF94"; + std::vector m_types; + Vector3 m_dipoleVector; + bool m_customDipole = false; // Custom dipole moment set +}; + +} // end namespace QtPlugins +} // end namespace Avogadro + +#endif // AVOGADRO_QTPLUGINS_FORCE_H diff --git a/avogadro/qtplugins/molecularproperties/molecularmodel.cpp b/avogadro/qtplugins/molecularproperties/molecularmodel.cpp index 73866a44d3..7577fc8caa 100644 --- a/avogadro/qtplugins/molecularproperties/molecularmodel.cpp +++ b/avogadro/qtplugins/molecularproperties/molecularmodel.cpp @@ -53,7 +53,8 @@ void MolecularModel::setMolecule(QtGui::Molecule* molecule) } // make sure we know if the molecule changed - connect(m_molecule, SIGNAL(changed(uint)), SLOT(updateTable(uint))); + connect(m_molecule, &QtGui::Molecule::changed, this, + &MolecularModel::updateTable); updateTable(QtGui::Molecule::Added); } @@ -291,6 +292,8 @@ QVariant MolecularModel::headerData(int section, Qt::Orientation orientation, return tr("Net Charge"); else if (it->first == " 10totalSpinMultiplicity") return tr("Net Spin Multiplicity"); + else if (it->first == "dipoleMoment") + return tr("Dipole Moment (Debye)"); else if (it->first == "homoEnergy") return tr("HOMO Energy (eV)", "highest occupied molecular orbital"); else if (it->first == "lumoEnergy") @@ -389,6 +392,10 @@ void MolecularModel::updateTable(unsigned int flags) if (m_molecule->totalSpinMultiplicity() != 1) m_propertiesCache.setValue(" 10totalSpinMultiplicity", m_molecule->totalSpinMultiplicity()); + if (m_molecule->hasData("dipoleMoment")) { + auto dipole = m_molecule->data("dipoleMoment").toVector3(); + m_propertiesCache.setValue("dipoleMoment", dipole.norm()); + } // TODO check for homo, lumo, or somo energies // m_propertiesCache.setValue("homoEnergy", energy); diff --git a/avogadro/qtplugins/openbabel/obcharges.cpp b/avogadro/qtplugins/openbabel/obcharges.cpp index b2c2c9d405..166310586e 100644 --- a/avogadro/qtplugins/openbabel/obcharges.cpp +++ b/avogadro/qtplugins/openbabel/obcharges.cpp @@ -156,9 +156,8 @@ std::string OBCharges::name() const return ""; } -MatrixX OBCharges::partialCharges(Core::Molecule& molecule) const +MatrixX OBCharges::partialCharges(const Core::Molecule& molecule) const { - // get the charges from obabel MatrixX charges(molecule.atomCount(), 1); if (m_identifier.empty()) { @@ -168,10 +167,13 @@ MatrixX OBCharges::partialCharges(Core::Molecule& molecule) const // check to see if we already have them in the molecule charges = molecule.partialCharges(m_identifier); - // if there's a non-zero charge, then we're done - for (unsigned int i = 0; i < charges.rows(); ++i) { - if (abs(charges(i, 0)) > 0.00001) - return charges; + // if the number of charges matches the number of atoms + // and there's a non-zero charge, then we're done + if (charges.rows() == molecule.atomCount()) { + for (unsigned int i = 0; i < charges.rows(); ++i) { + if (abs(charges(i, 0)) > 0.00001) + return charges; + } } // otherwise, we're going to run obprocess to get the charges @@ -202,6 +204,18 @@ MatrixX OBCharges::partialCharges(Core::Molecule& molecule) const if (abs(charges(0, 0)) < 0.00001) charges(0, 0) = 0.0001; + // check the size + if (output.size() != molecule.atomCount()) { + qDebug() << "Charges size mismatch."; + return charges; + } + return charges; +} + +MatrixX OBCharges::partialCharges(Core::Molecule& molecule) const +{ + MatrixX charges = + partialCharges(static_cast(molecule)); // cache the charges and allow them to show up in output molecule.setPartialCharges(m_identifier, charges); return charges; diff --git a/avogadro/qtplugins/openbabel/obcharges.h b/avogadro/qtplugins/openbabel/obcharges.h index 23e3b5f052..4727b2c72a 100644 --- a/avogadro/qtplugins/openbabel/obcharges.h +++ b/avogadro/qtplugins/openbabel/obcharges.h @@ -55,6 +55,7 @@ class OBCharges : public Avogadro::Calc::ChargeModel * @brief Retrieve the relevant charges from the molecule for our defined type */ virtual MatrixX partialCharges(Core::Molecule& mol) const override; + virtual MatrixX partialCharges(const Core::Molecule& mol) const override; /** * @brief Synchronous use of the OBProcess. diff --git a/avogadro/qtplugins/openbabel/obprocess.cpp b/avogadro/qtplugins/openbabel/obprocess.cpp index 108828f4e8..d2d924fc78 100644 --- a/avogadro/qtplugins/openbabel/obprocess.cpp +++ b/avogadro/qtplugins/openbabel/obprocess.cpp @@ -570,9 +570,11 @@ void OBProcess::executeObabel(const QStringList& options, QObject* receiver, } // Start process +#ifndef NDEBUG qDebug() << "OBProcess::executeObabel: " "Running" << m_obabelExecutable << options.join(" "); +#endif m_process->start(m_obabelExecutable, options); if (!obabelStdin.isNull()) { m_process->write(obabelStdin); diff --git a/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp b/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp index 92ab999e2f..6802ef8e65 100644 --- a/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp +++ b/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp @@ -49,7 +49,7 @@ Calc::ChargeModel* ScriptChargeModel::newInstance() const return new ScriptChargeModel(m_interpreter->scriptFilePath()); } -MatrixX ScriptChargeModel::partialCharges(Core::Molecule& mol) const +MatrixX ScriptChargeModel::partialCharges(const Core::Molecule& mol) const { MatrixX charges(mol.atomCount(), 1); @@ -108,10 +108,15 @@ MatrixX ScriptChargeModel::partialCharges(Core::Molecule& mol) const charges(atom, 0) = charge; ++atom; } + return charges; +} - // cache the charges +MatrixX ScriptChargeModel::partialCharges(Core::Molecule& mol) const +{ + // just create a copy of the const version + MatrixX charges = partialCharges(static_cast(mol)); + // cache them mol.setPartialCharges(m_identifier, charges); - return charges; } diff --git a/avogadro/qtplugins/scriptcharges/scriptchargemodel.h b/avogadro/qtplugins/scriptcharges/scriptchargemodel.h index 290c2f10bb..218be67e57 100644 --- a/avogadro/qtplugins/scriptcharges/scriptchargemodel.h +++ b/avogadro/qtplugins/scriptcharges/scriptchargemodel.h @@ -59,6 +59,7 @@ class ScriptChargeModel : public Avogadro::Calc::ChargeModel Core::Molecule::ElementMask elements() const override { return m_elements; } MatrixX partialCharges(Core::Molecule& mol) const override; + MatrixX partialCharges(const Core::Molecule& mol) const override; double potential(Core::Molecule& mol, const Vector3& point) const override; diff --git a/avogadro/quantumio/gaussianfchk.cpp b/avogadro/quantumio/gaussianfchk.cpp index 496c90cd72..3a85ee162b 100644 --- a/avogadro/quantumio/gaussianfchk.cpp +++ b/avogadro/quantumio/gaussianfchk.cpp @@ -73,7 +73,10 @@ bool GaussianFchk::read(std::istream& in, Core::Molecule& molecule) // set the spin multiplicity molecule.setData("totalSpinMultiplicity", m_spin); // dipole moment - molecule.setData("dipoleMoment", m_dipoleMoment); + // TODO: This should be a Vector3d + Core::Variant dipole(m_dipoleMoment.x(), m_dipoleMoment.y(), + m_dipoleMoment.z()); + molecule.setData("dipoleMoment", dipole); // Do simple bond perception. molecule.perceiveBondsSimple(); diff --git a/avogadro/quantumio/orca.h b/avogadro/quantumio/orca.h index caeaff9501..006872a0b9 100644 --- a/avogadro/quantumio/orca.h +++ b/avogadro/quantumio/orca.h @@ -62,7 +62,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat std::vector m_atomNums; std::vector m_atomPos; - Eigen::Vector3d m_dipoleMoment; + Vector3 m_dipoleMoment; std::vector> m_bondOrders; diff --git a/avogadro/rendering/arrowgeometry.cpp b/avogadro/rendering/arrowgeometry.cpp index df99b967d1..b1aac510fb 100644 --- a/avogadro/rendering/arrowgeometry.cpp +++ b/avogadro/rendering/arrowgeometry.cpp @@ -41,12 +41,17 @@ class ArrowGeometry::Private ShaderProgram program; }; -ArrowGeometry::ArrowGeometry() : m_dirty(false), d(new Private) {} +ArrowGeometry::ArrowGeometry() + : m_dirty(false), d(new Private), m_color(0, 255, 0) +{ +} ArrowGeometry::ArrowGeometry(const ArrowGeometry& other) : Drawable(other), m_vertices(other.m_vertices), - m_lineStarts(other.m_lineStarts), m_dirty(true), d(new Private) -{} + m_lineStarts(other.m_lineStarts), m_color(other.m_color), m_dirty(true), + d(new Private) +{ +} ArrowGeometry::~ArrowGeometry() { @@ -95,10 +100,8 @@ void ArrowGeometry::render(const Camera& camera) } // Render the arrows using the shader. - for (auto & m_vertice : m_vertices) { - Vector3f v3 = - m_vertice.first + - 0.8 * (m_vertice.second - m_vertice.first); + for (auto& m_vertice : m_vertices) { + Vector3f v3 = m_vertice.first + 0.8 * (m_vertice.second - m_vertice.first); drawLine(m_vertice.first, v3, 2); drawCone(v3, m_vertice.second, 0.05, 1.0); } @@ -117,7 +120,6 @@ void ArrowGeometry::drawLine(const Vector3f& start, const Vector3f& end, double lineWidth) { // Draw a line between two points of the specified thickness - glPushAttrib(GL_LIGHTING_BIT); glDisable(GL_LIGHTING); @@ -158,7 +160,7 @@ void ArrowGeometry::drawCone(const Vector3f& base, const Vector3f& cap, Vector3f n = (cap - v).cross(v - vPrec).normalized(); Vector3f nNext = (cap - vNext).cross(vNext - v).normalized(); glBegin(GL_TRIANGLES); - glColor3ub(0, 255, 0); + glColor3ub(m_color[0], m_color[1], m_color[2]); glNormal3fv((n + nNext).normalized().data()); glVertex3fv(cap.data()); glNormal3fv(nNext.data()); @@ -188,4 +190,4 @@ void ArrowGeometry::addSingleArrow(const Vector3f& pos1, const Vector3f& pos2) m_dirty = true; } -} // End namespace Avogadro +} // namespace Avogadro::Rendering diff --git a/avogadro/rendering/arrowgeometry.h b/avogadro/rendering/arrowgeometry.h index 35807ffc3f..10e72de731 100644 --- a/avogadro/rendering/arrowgeometry.h +++ b/avogadro/rendering/arrowgeometry.h @@ -66,6 +66,9 @@ class AVOGADRORENDERING_EXPORT ArrowGeometry : public Drawable return m_vertices; } + /** Set the color of the arrow */ + void setColor(const Vector3ub& c) { m_color = c; } + private: /** * @brief Update the shaders ready for rendering. @@ -74,6 +77,7 @@ class AVOGADRORENDERING_EXPORT ArrowGeometry : public Drawable Core::Array> m_vertices; Core::Array m_lineStarts; + Vector3ub m_color; bool m_dirty;