Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Calculate and render dipole moments #1801

Merged
merged 16 commits into from
Nov 20, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 48 additions & 10 deletions avogadro/calc/chargemanager.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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];

Expand All @@ -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;
}
Expand Down Expand Up @@ -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
Expand Down
17 changes: 17 additions & 0 deletions avogadro/calc/chargemanager.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
25 changes: 25 additions & 0 deletions avogadro/calc/chargemodel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
#include <avogadro/core/array.h>
#include <avogadro/core/molecule.h>

#include <iostream>

namespace Avogadro {

using Core::Array;
Expand All @@ -19,6 +21,27 @@ namespace Calc {
constexpr double M_PI = 3.14159265358979323846;
#endif

Vector3 ChargeModel::dipoleMoment(const Molecule& mol) const
{
// default is to get the set of partial atomic charges
// (some models might do something more sophisticated)
const MatrixX charges = partialCharges(mol);
// also get the positions of the atoms
const Array<Vector3> positions = mol.atomPositions3d();

Vector3 dipole(0.0, 0.0, 0.0);
if (charges.rows() != positions.size())
std::cout << "Error: charges " << charges.rows() << " != positions "
<< positions.size() << std::endl;

for (unsigned int i = 0; i < charges.size(); ++i)
dipole += charges(i, 0) * positions[i];

dipole *= 4.80320471257; // convert to Debye from electron-Angstrom

return dipole;
}

double ChargeModel::potential(Molecule& mol, const Vector3& point) const
{
// default is to get the set of partial atomic charges
Expand All @@ -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();
Expand Down
11 changes: 11 additions & 0 deletions avogadro/calc/chargemodel.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
5 changes: 5 additions & 0 deletions avogadro/calc/defaultmodel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 2 additions & 0 deletions avogadro/calc/defaultmodel.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
88 changes: 88 additions & 0 deletions avogadro/core/variant-inline.h
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,19 @@

#include "variant.h"

#include <iostream>
#include <sstream>

namespace Avogadro::Core {

inline Variant::Variant() : m_type(Null) {}

inline Variant::Variant(double x, double y, double z) : m_type(Vector)
{
Vector3* v = new Vector3(x, y, z);
m_value.vector = v;
}

template <typename T>
inline Variant::Variant(T v) : m_type(Null)
{
Expand All @@ -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;
}
Expand All @@ -59,11 +82,31 @@ inline bool Variant::isNull() const
return m_type == Null;
}

inline bool Variant::setValue(double x, double y, double z)
{
clear();

m_type = Vector;
m_value.vector = new Vector3(x, y, z);

return true;
}

template <typename T>
inline bool Variant::setValue(T v)
{
AVO_UNUSED(v);

#ifndef NDEBUG
#if defined(_MSC_VER)
std::cerr << " Variant::setValue() not implemented for " << __FUNCSIG__
<< std::endl;
#else
std::cerr << " Variant::setValue() not implemented for "
<< __PRETTY_FUNCTION__ << std::endl;
#endif
#endif

clear();

return false;
Expand Down Expand Up @@ -185,6 +228,28 @@ inline bool Variant::setValue(MatrixX matrix)
return true;
}

template <>
inline bool Variant::setValue(Vector3 vector)
{
clear();

m_type = Vector;
m_value.vector = new Vector3(vector);

return true;
}

template <>
inline bool Variant::setValue(Vector3f vector)
{
clear();

m_type = Vector;
m_value.vector = new Vector3(vector.x(), vector.y(), vector.z());

return true;
}

template <typename T>
inline T Variant::value() const
{
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -424,6 +507,11 @@ inline const MatrixX& Variant::toMatrixRef() const
return value<const MatrixX&>();
}

inline Vector3 Variant::toVector3() const
{
return value<Vector3>();
}

// --- Operators ----------------------------------------------------------- //
inline Variant& Variant::operator=(const Variant& variant)
{
Expand Down
Loading
Loading