Skip to content

Commit

Permalink
Add support for unit cell optimizing through minimum distances
Browse files Browse the repository at this point in the history 
Signed-off-by: Geoff Hutchison <[email protected]>
  • Loading branch information
@ghutchis
ghutchis committed Oct 12, 2023
1 parent f33a553 commit 05bd338
Show file tree
Hide file tree
Showing 4 changed files with 110 additions and 45 deletions.
124 changes: 89 additions & 35 deletions avogadro/calc/lennardjones.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,14 @@

#include <avogadro/core/elements.h>
#include <avogadro/core/molecule.h>
#include <avogadro/core/unitcell.h>

namespace Avogadro::Calc {

LennardJones::LennardJones() : m_vdw(true), m_depth(100.0), m_exponent(6) {}
LennardJones::LennardJones()
: m_vdw(true), m_depth(100.0), m_exponent(6), m_cell(nullptr)
{
}

LennardJones::~LennardJones() {}

Expand All @@ -21,11 +25,16 @@ void LennardJones::setMolecule(Core::Molecule* mol)
if (mol == nullptr) {
return; // nothing to do
}

m_cell = mol->unitCell(); // could be nullptr
int numAtoms = mol->atomCount();

// track atomic radii for this molecule
m_radii.setZero();
Eigen::MatrixXd radii(numAtoms, numAtoms);
Eigen::MatrixXd mask(numAtoms*3, 1);
mask.setOnes();
m_mask = mask;

for (Index i = 0; i < numAtoms; ++i) {
Core::Atom atom1 = mol->atom(i);
Expand All @@ -50,7 +59,6 @@ void LennardJones::setMolecule(Core::Molecule* mol)
}

m_radii = radii;
//@todo store unit cell if available
}

Real LennardJones::value(const Eigen::VectorXd& x)
Expand All @@ -63,16 +71,34 @@ Real LennardJones::value(const Eigen::VectorXd& x)
int numAtoms = m_molecule->atomCount();

Real energy = 0.0;
for (Index i = 0; i < numAtoms; ++i) {
Vector3 ipos(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
for (Index j = i + 1; j < numAtoms; ++j) {
Vector3 jpos(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Real r = (ipos - jpos).norm();
if (r < 0.1)
r = 0.1; // ensure we don't divide by zero

Real ratio = pow((m_radii(i, j) / r), m_exponent);
energy += m_depth * (ratio * ratio - 2.0 * (ratio));
// we put the conditional here outside the double loop
if (m_cell == nullptr) {
// regular molecule
for (Index i = 0; i < numAtoms; ++i) {
Vector3 ipos(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
for (Index j = i + 1; j < numAtoms; ++j) {
Vector3 jpos(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Real r = (ipos - jpos).norm();
if (r < 0.1)
r = 0.1; // ensure we don't divide by zero

Real ratio = pow((m_radii(i, j) / r), m_exponent);
energy += m_depth * (ratio * ratio - 2.0 * (ratio));
}
}
} else {
// use the unit cell to get minimum distances
for (Index i = 0; i < numAtoms; ++i) {
Vector3 ipos(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
for (Index j = i + 1; j < numAtoms; ++j) {
Vector3 jpos(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Real r = m_cell->distance(ipos, jpos);
if (r < 0.1)
r = 0.1; // ensure we don't divide by zero

Real ratio = pow((m_radii(i, j) / r), m_exponent);
energy += m_depth * (ratio * ratio - 2.0 * (ratio));
}
}
}

Expand All @@ -88,30 +114,58 @@ void LennardJones::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad)
// clear the gradients
grad.setZero();

//@todo handle unit cells and minimum distances
int numAtoms = m_molecule->atomCount();

for (Index i = 0; i < numAtoms; ++i) {
Vector3 ipos(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
for (Index j = i + 1; j < numAtoms; ++j) {
Vector3 jpos(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3 force = ipos - jpos;

Real r = force.norm();
if (r < 0.1)
r = 0.1; // ensure we don't divide by zero

Real rad = pow(m_radii(i, j), m_exponent);
Real term1 = -2 * (m_exponent)*rad * rad * pow(r, -2 * m_exponent - 1);
Real term2 = 2 * (m_exponent)*rad * pow(r, -1 * m_exponent - 1);
Real dE = m_depth * (term1 + term2);

force = (dE / r) * force;

// update gradients
for (unsigned int c = 0; c < 3; ++c) {
grad[3 * i + c] += force[c];
grad[3 * j + c] -= force[c];
// we put the conditional here outside the double loop
if (m_cell == nullptr) {
// regular molecule
for (Index i = 0; i < numAtoms; ++i) {
Vector3 ipos(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
for (Index j = i + 1; j < numAtoms; ++j) {
Vector3 jpos(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3 force = ipos - jpos;

Real r = force.norm();
if (r < 0.1)
r = 0.1; // ensure we don't divide by zero

Real rad = pow(m_radii(i, j), m_exponent);
Real term1 = -2 * (m_exponent)*rad * rad * pow(r, -2 * m_exponent - 1);
Real term2 = 2 * (m_exponent)*rad * pow(r, -1 * m_exponent - 1);
Real dE = m_depth * (term1 + term2);

force = (dE / r) * force;

// update gradients
for (unsigned int c = 0; c < 3; ++c) {
grad[3 * i + c] += force[c];
grad[3 * j + c] -= force[c];
}
}
}
} else {
// unit cell
for (Index i = 0; i < numAtoms; ++i) {
Vector3 ipos(x[3 * i], x[3 * i + 1], x[3 * i + 2]);
for (Index j = i + 1; j < numAtoms; ++j) {
Vector3 jpos(x[3 * j], x[3 * j + 1], x[3 * j + 2]);
Vector3 force = m_cell->minimumImage(ipos - jpos);

Real r = force.norm();
if (r < 0.1)
r = 0.1; // ensure we don't divide by zero

Real rad = pow(m_radii(i, j), m_exponent);
Real term1 = -2 * (m_exponent)*rad * rad * pow(r, -2 * m_exponent - 1);
Real term2 = 2 * (m_exponent)*rad * pow(r, -1 * m_exponent - 1);
Real dE = m_depth * (term1 + term2);

force = (dE / r) * force;

// update gradients
for (unsigned int c = 0; c < 3; ++c) {
grad[3 * i + c] += force[c];
grad[3 * j + c] -= force[c];
}
}
}
}
Expand Down
2 changes: 2 additions & 0 deletions avogadro/calc/lennardjones.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
namespace Avogadro {
namespace Core {
class Molecule;
class UnitCell;
}

namespace Calc {
Expand Down Expand Up @@ -43,6 +44,7 @@ class AVOGADROCALC_EXPORT LennardJones : public EnergyCalculator

protected:
Core::Molecule* m_molecule;
Core::UnitCell* m_cell;
Eigen::MatrixXd m_radii;
bool m_vdw;
Real m_depth;
Expand Down
2 changes: 1 addition & 1 deletion avogadro/core/unitcell.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,7 +342,7 @@ inline Vector3 UnitCell::minimumImage(const Vector3& v) const

inline Real UnitCell::distance(const Vector3& v1, const Vector3& v2) const
{
return std::fabs(minimumImage(v1 - v2).norm());
return minimumImage(v1 - v2).norm();
}

} // namespace Core
Expand Down
27 changes: 18 additions & 9 deletions avogadro/qtplugins/forcefield/forcefield.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,13 +60,15 @@ Forcefield::Forcefield(QObject* parent_)
action->setEnabled(true);
action->setText(tr("Optimize"));
action->setData(optimizeAction);
action->setProperty("menu priority", 920);
connect(action, SIGNAL(triggered()), SLOT(optimize()));
m_actions.push_back(action);

action = new QAction(this);
action->setEnabled(true);
action->setText(tr("Energy")); // calculate energy
action->setData(energyAction);
action->setProperty("menu priority", 910);
connect(action, SIGNAL(triggered()), SLOT(energy()));
m_actions.push_back(action);

Expand All @@ -93,7 +95,7 @@ QStringList Forcefield::menuPath(QAction* action) const
path << tr("&Extensions");
return path;
}
path << tr("&Extensions") << tr("Calculate");
path << tr("&Extensions") << tr("&Calculate");
return path;
}

Expand Down Expand Up @@ -122,20 +124,23 @@ void Forcefield::optimize()
m_molecule->undoMolecule()->setInteractive(true);

//@todo check m_minimizer for method to use
cppoptlib::LbfgsSolver<EnergyCalculator> solver;
//cppoptlib::LbfgsSolver<EnergyCalculator> solver;
cppoptlib::ConjugatedGradientDescentSolver<EnergyCalculator> solver;

int n = m_molecule->atomCount();
Core::Array<Vector3> pos = m_molecule->atomPositions3d();
// just to get the right size / shape
Core::Array<Vector3> forces = m_molecule->atomPositions3d();
double* p = pos[0].data();
Eigen::Map<Eigen::VectorXd> map(p, 3 * n);
Eigen::VectorXd positions = map;
Eigen::VectorXd gradient = Eigen::VectorXd::Zero(3 * n);

// Create a Criteria class to adjust stopping criteria
cppoptlib::Criteria<Real> crit = cppoptlib::Criteria<Real>::defaults();
// @todo allow criteria to be set
crit.iterations = 5;
crit.xDelta = 1.0e-4; // positions converged to 1.0e-4
crit.fDelta = 1.0e-4; // energy converged to 1.0e-4
crit.fDelta = 1.0e-6;
solver.setStopCriteria(crit); // every 5 steps, update coordinates
cppoptlib::Status status = cppoptlib::Status::NotStarted;

Expand All @@ -148,25 +153,29 @@ void Forcefield::optimize()
solver.minimize(lj, positions);

cppoptlib::Status currentStatus = solver.status();
if (currentStatus != status || currentStatus == cppoptlib::Status::Continue) {
// status has changed or minimizer hasn't converged
// qDebug() << " status: " << (int)currentStatus;
// qDebug() << " energy: " << lj.value(positions);
// check the gradient norm
lj.gradient(positions, gradient);
// qDebug() << " gradient: " << gradient.norm();

// update coordinates
const double* d = positions.data();
// casting would be lovely...
for (size_t i = 0; i < n; ++i) {
pos[i] = Vector3(*(d), *(d + 1), *(d + 2));
d += 3;

forces[i] = Vector3(gradient[3 * i], gradient[3 * i + 1],
gradient[3 * i + 2]);
}
// todo - merge these into one undo step
m_molecule->undoMolecule()->setAtomPositions3d(pos, tr("Optimize Geometry"));
m_molecule->setForceVectors(forces);
Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified;
m_molecule->emitChanged(changes);
}
}

qDebug() << " energy: " << lj.value(positions);

m_molecule->undoMolecule()->setInteractive(isInteractive);
}

Expand Down

0 comments on commit 05bd338

Please sign in to comment.