diff --git a/avogadro/calc/lennardjones.cpp b/avogadro/calc/lennardjones.cpp index b4d0c1f8ab..63e90e0c58 100644 --- a/avogadro/calc/lennardjones.cpp +++ b/avogadro/calc/lennardjones.cpp @@ -7,10 +7,14 @@ #include #include +#include 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() {} @@ -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); @@ -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) @@ -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)); + } } } @@ -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]; + } } } } diff --git a/avogadro/calc/lennardjones.h b/avogadro/calc/lennardjones.h index d54f9151f6..3364ff3f3e 100644 --- a/avogadro/calc/lennardjones.h +++ b/avogadro/calc/lennardjones.h @@ -13,6 +13,7 @@ namespace Avogadro { namespace Core { class Molecule; +class UnitCell; } namespace Calc { @@ -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; diff --git a/avogadro/core/unitcell.h b/avogadro/core/unitcell.h index 99ef4de1db..8565f95a2a 100644 --- a/avogadro/core/unitcell.h +++ b/avogadro/core/unitcell.h @@ -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 diff --git a/avogadro/qtplugins/forcefield/forcefield.cpp b/avogadro/qtplugins/forcefield/forcefield.cpp index a6ab813d91..ca20d539cf 100644 --- a/avogadro/qtplugins/forcefield/forcefield.cpp +++ b/avogadro/qtplugins/forcefield/forcefield.cpp @@ -60,6 +60,7 @@ 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); @@ -67,6 +68,7 @@ Forcefield::Forcefield(QObject* parent_) 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); @@ -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; } @@ -122,20 +124,23 @@ void Forcefield::optimize() m_molecule->undoMolecule()->setInteractive(true); //@todo check m_minimizer for method to use - cppoptlib::LbfgsSolver solver; + //cppoptlib::LbfgsSolver solver; + cppoptlib::ConjugatedGradientDescentSolver solver; int n = m_molecule->atomCount(); Core::Array pos = m_molecule->atomPositions3d(); + // just to get the right size / shape + Core::Array forces = m_molecule->atomPositions3d(); double* p = pos[0].data(); Eigen::Map 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 crit = cppoptlib::Criteria::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; @@ -148,8 +153,11 @@ 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(); @@ -157,16 +165,17 @@ void Forcefield::optimize() 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); }