diff --git a/avogadro/qtplugins/forcefield/forcefield.cpp b/avogadro/qtplugins/forcefield/forcefield.cpp index 31c7095ee1..285a9e2920 100644 --- a/avogadro/qtplugins/forcefield/forcefield.cpp +++ b/avogadro/qtplugins/forcefield/forcefield.cpp @@ -31,9 +31,6 @@ #include #include -#include -#include -#include #include namespace Avogadro { @@ -49,6 +46,7 @@ const int configureAction = 2; const int freezeAction = 3; const int unfreezeAction = 4; const int constraintAction = 5; +const int forcesAction = 6; Forcefield::Forcefield(QObject* parent_) : ExtensionPlugin(parent_), m_method(nullptr) @@ -63,12 +61,12 @@ Forcefield::Forcefield(QObject* parent_) m_gradientTolerance = settings.value("gradientTolerance", 1.0e-4).toDouble(); settings.endGroup(); - refreshScripts(); - /* @todo - finish OBMM interface + // add the openbabel calculators in case they don't exist Calc::EnergyManager::registerModel(new OBMMEnergy("MMFF94")); Calc::EnergyManager::registerModel(new OBMMEnergy("UFF")); Calc::EnergyManager::registerModel(new OBMMEnergy("GAFF")); - */ + + refreshScripts(); QAction* action = new QAction(this); action->setEnabled(true); @@ -86,6 +84,14 @@ Forcefield::Forcefield(QObject* parent_) connect(action, SIGNAL(triggered()), SLOT(energy())); m_actions.push_back(action); + action = new QAction(this); + action->setEnabled(true); + action->setText(tr("Forces")); // calculate gradients + action->setData(forcesAction); + action->setProperty("menu priority", 910); + connect(action, SIGNAL(triggered()), SLOT(forces())); + m_actions.push_back(action); + action = new QAction(this); action->setEnabled(true); action->setText(tr("Configure…")); @@ -123,7 +129,11 @@ QList Forcefield::actions() const QStringList Forcefield::menuPath(QAction* action) const { QStringList path; - path << tr("&Extensions") << tr("&Calculate"); + if (action->data().toInt() == optimizeAction) + path << tr("&Extensions"); + else + path << tr("&Extensions") << tr("&Calculate"); + return path; } @@ -220,7 +230,6 @@ void Forcefield::optimize() m_molecule->undoMolecule()->setInteractive(true); cppoptlib::LbfgsSolver solver; - // cppoptlib::ConjugatedGradientDescentSolver solver; int n = m_molecule->atomCount(); @@ -252,7 +261,7 @@ void Forcefield::optimize() cppoptlib::Criteria crit = cppoptlib::Criteria::defaults(); // e.g., every N steps, update coordinates - crit.iterations = 5; + crit.iterations = 2; // we don't set function or gradient criteria // .. these seem to be broken in the solver code // .. so we handle ourselves @@ -295,6 +304,11 @@ void Forcefield::optimize() forces[i] = -0.1 * Vector3(gradient[3 * i], gradient[3 * i + 1], gradient[3 * i + 2]); } + } else { + // reset to last positions + positions = lastPositions; + gradient = Eigen::VectorXd::Zero(3 * n); + break; } // todo - merge these into one undo step @@ -307,18 +321,12 @@ void Forcefield::optimize() lastPositions = positions; // check for convergence - /* if (fabs(gradient.maxCoeff()) < m_gradientTolerance) break; if (fabs(currentEnergy - energy) < m_tolerance) break; - */ energy = currentEnergy; - } else { - // reset to last positions - positions = lastPositions; - gradient = Eigen::VectorXd::Zero(3 * n); } } @@ -345,6 +353,52 @@ void Forcefield::energy() QMessageBox::information(nullptr, tr("Avogadro"), msg); } +void Forcefield::forces() +{ + if (m_molecule == nullptr || m_method == nullptr) + return; + + int n = m_molecule->atomCount(); + + // double-check the mask + auto mask = m_molecule->frozenAtomMask(); + if (mask.rows() != 3 * n) { + mask = Eigen::VectorXd::Zero(3 * n); + // set to 1.0 + for (Index i = 0; i < 3 * n; ++i) { + mask[i] = 1.0; + } + } + m_method->setMolecule(m_molecule); + m_method->setMask(mask); + + // we have to cast the current 3d positions into a VectorXd + Core::Array pos = 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); + // just to get the right size / shape + // we'll use this to draw the force arrows + Core::Array forces = m_molecule->atomPositions3d(); + + m_method->gradient(positions, gradient); + + for (size_t i = 0; i < n; ++i) { + forces[i] = + -0.1 * Vector3(gradient[3 * i], gradient[3 * i + 1], gradient[3 * i + 2]); + } + + m_molecule->setForceVectors(forces); + Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified; + m_molecule->emitChanged(changes); + + QString msg( + tr("%1 Force Norm = %L2").arg(m_methodName.c_str()).arg(gradient.norm())); + QMessageBox::information(nullptr, tr("Avogadro"), msg); +} + std::string Forcefield::recommendedForceField() const { // if we have a unit cell, we need to use the LJ calculator diff --git a/avogadro/qtplugins/forcefield/forcefield.h b/avogadro/qtplugins/forcefield/forcefield.h index 05170b66cf..a08e75675c 100644 --- a/avogadro/qtplugins/forcefield/forcefield.h +++ b/avogadro/qtplugins/forcefield/forcefield.h @@ -71,6 +71,7 @@ public slots: private slots: void energy(); + void forces(); void optimize(); void freezeSelected(); void unfreezeSelected(); diff --git a/avogadro/qtplugins/forcefield/obmmenergy.cpp b/avogadro/qtplugins/forcefield/obmmenergy.cpp index e58d798243..594dba9f80 100644 --- a/avogadro/qtplugins/forcefield/obmmenergy.cpp +++ b/avogadro/qtplugins/forcefield/obmmenergy.cpp @@ -12,58 +12,14 @@ #include #include #include +#include #include namespace Avogadro::QtPlugins { -class OBMMEnergy::ProcessListener : public QObject -{ - Q_OBJECT -public: - ProcessListener(QProcess* proc) - : QObject(), m_finished(false), m_process(proc) - { - } - - bool waitForOutput(QByteArray output, int msTimeout = 2000) - { - connect(m_process, SIGNAL(readyRead()), SLOT(readyRead())); - if (!wait(msTimeout)) - return false; - - // success! - output = m_output; - disconnect(m_process, nullptr, nullptr, nullptr); - m_finished = false; - return true; - } - -public slots: - void readyRead() - { - m_finished = true; - m_output = m_process->readAllStandardOutput(); - } - -private: - bool wait(int msTimeout) - { - QTimer timer; - timer.start(msTimeout); - - while (timer.isActive() && !m_finished) - qApp->processEvents(QEventLoop::AllEvents, 500); - - return m_finished; - } - - QProcess* m_process; - bool m_finished; - QByteArray m_output; -}; - OBMMEnergy::OBMMEnergy(const std::string& method) : m_identifier(method), m_name(method), m_process(nullptr), + m_molecule(nullptr), #if defined(_WIN32) m_executable("obmm.exe") #else @@ -117,6 +73,33 @@ OBMMEnergy::~OBMMEnergy() delete m_process; } +QByteArray OBMMEnergy::writeAndRead(const QByteArray& input) +{ + if (m_process == nullptr) + return QByteArray(); + + QByteArray result, line; + m_process->write(input + "\n"); + QThread::msleep(1); + m_process->waitForReadyRead(5); + while (m_process->canReadLine() && !line.startsWith("command >")) { + line = m_process->readLine(); + result += line; + } + // check if we've really flushed the output + if (!result.contains("invalid command\n command >")) { + m_process->write(" \n"); + QThread::msleep(1); + m_process->waitForReadyRead(5); + while (m_process->canReadLine()) { + line = m_process->readLine(); + result += line; + } + } + result += m_process->readAllStandardOutput(); + return result; +} + void OBMMEnergy::setupProcess() { if (m_process != nullptr) { @@ -158,6 +141,8 @@ void OBMMEnergy::setupProcess() env.insert("BABEL_LIBDIR", QCoreApplication::applicationDirPath() + "/../lib/openbabel/" + dirs[0]); } else { + env.insert("BABEL_LIBDIR", QCoreApplication::applicationDirPath() + + "/../lib/openbabel/"); qDebug() << "Error, Open Babel plugins directory not found."; } #endif @@ -177,7 +162,7 @@ void OBMMEnergy::setMolecule(Core::Molecule* mol) // should check if the molecule is valid for this script // .. this should never happen, but let's be defensive - if (mol == nullptr) { + if (mol == nullptr || mol->atomCount() == 0) { return; // nothing to do } @@ -204,44 +189,28 @@ void OBMMEnergy::setMolecule(Core::Molecule* mol) m_process->start(m_executable, QStringList() << m_tempFile.fileName()); if (!m_process->waitForStarted()) { // appendError("Error starting process."); + qDebug() << "OBMM: Error starting process."; return; } - ProcessListener listener(m_process); - QByteArray result; - if (!listener.waitForOutput(result)) { - // appendError("Error running process."); - return; - } - qDebug() << "OBMM start: " << result; - - // okay, we need to write "load " to the interpreter - // and then read the response - QByteArray input = "load " + m_tempFile.fileName().toLocal8Bit() + "\n"; - m_process->write(input); - if (!listener.waitForOutput(result)) { - // appendError("Error running process."); - return; - } - qDebug() << "OBMM: " << result; - - // set the method m_identifier.c_str() + - input = QByteArray("ff MMFF94\n"); - m_process->write(input); - if (!listener.waitForOutput(result)) { - // appendError("Error running process."); - return; - } - qDebug() << "OBMM ff: " << result; - // check for an energy - input = QByteArray("energy\n"); - result.clear(); - m_process->write(input); + QByteArray input, line, result; + m_process->waitForReadyRead(); result.clear(); while (!result.contains("command >")) { result += m_process->readLine(); + if (!m_process->canReadLine()) + break; } - qDebug() << "OBMM energy: " << result; + result += m_process->readAllStandardOutput(); + + // set the method m_identifier.c_str() + + input = QByteArray("ff ") + m_identifier.c_str(); + result = writeAndRead(input); + + // okay, we need to write "load " to the interpreter + // and then read the response + input = "load " + m_tempFile.fileName().toLocal8Bit(); + result = writeAndRead(input); } Real OBMMEnergy::value(const Eigen::VectorXd& x) @@ -249,56 +218,33 @@ Real OBMMEnergy::value(const Eigen::VectorXd& x) if (m_molecule == nullptr || m_process == nullptr) return 0.0; // nothing to do - m_process->waitForReadyRead(); - QByteArray result; - while (!result.contains("command >")) { - result += m_process->readLine(); - } - qDebug() << " starting " << result; + QByteArray input, result; // write the new coordinates and read the energy - QByteArray input = "coord\n"; + input = "coord\n"; for (Index i = 0; i < x.size(); i += 3) { // write as x y z (space separated) input += QString::number(x[i]) + " " + QString::number(x[i + 1]) + " " + QString::number(x[i + 2]) + "\n"; } - input += "\n"; - m_process->write(input); - m_process->waitForReadyRead(); - result.clear(); - while (!result.contains("command >")) { - result += m_process->readLine(); - } - - qDebug() << " asking energy " << result << m_process->state(); + result = writeAndRead(input); // now ask for the energy - input = "energy\n\n"; - result.clear(); - m_process->write(input); - m_process->waitForReadyRead(); - while (!result.contains("command >")) { - result += m_process->readLine(); - } - - qDebug() << "OBMM: " << result; + input = "energy\n"; + result = writeAndRead(input); // go through lines in result until we see "total energy" QStringList lines = QString(result).remove('\r').split('\n'); double energy = 0.0; for (auto line : lines) { if (line.contains("total energy =")) { - qDebug() << " OBMM: " << line; QStringList items = line.split(" ", Qt::SkipEmptyParts); if (items.size() > 4) energy = items[3].toDouble(); } } - qDebug() << " OBMM: " << energy << " done"; - return energy; // if conversion fails, returns 0.0 } @@ -307,38 +253,19 @@ void OBMMEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) if (m_molecule == nullptr || m_process == nullptr) return; - EnergyCalculator::gradient(x, grad); - return; - - qDebug() << "OBMM: gradient"; - // write the new coordinates and read the energy - QByteArray input = "coord\n"; + QByteArray result, input = "coord\n"; for (Index i = 0; i < x.size(); i += 3) { // write as x y z (space separated) input += QString::number(x[i]) + " " + QString::number(x[i + 1]) + " " + QString::number(x[i + 2]) + "\n"; } - m_process->write(input); - qDebug() << "OBMM Grad wrote coords"; - m_process->waitForReadyRead(); - QByteArray result; - while (m_process->canReadLine()) { - result += m_process->readLine(); - } - - qDebug() << "OBMM: " << result; + result = writeAndRead(input); // now ask for the energy - input = "grad\n"; - m_process->write(input); - m_process->waitForReadyRead(); - while (m_process->canReadLine()) { - result += m_process->readLine(); - } - - qDebug() << "OBMM: " << result; + input = "grad"; + result = writeAndRead(input); // go through lines in result until we see "gradient " QStringList lines = QString(result).remove('\r').split('\n'); @@ -352,7 +279,7 @@ void OBMMEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) if (readingGradient) { QStringList items = line.split(" ", Qt::SkipEmptyParts); if (items.size() == 3) { - grad[3 * i] = -1.0 * items[0].toDouble(); + grad[3 * i] = items[0].toDouble(); grad[3 * i + 1] = items[1].toDouble(); grad[3 * i + 2] = items[2].toDouble(); ++i; @@ -360,6 +287,8 @@ void OBMMEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) } } + grad *= -1; // OpenBabel outputs forces, not grads + cleanGradients(grad); } diff --git a/avogadro/qtplugins/forcefield/obmmenergy.h b/avogadro/qtplugins/forcefield/obmmenergy.h index f12e7ba1fd..63d0d9d2b3 100644 --- a/avogadro/qtplugins/forcefield/obmmenergy.h +++ b/avogadro/qtplugins/forcefield/obmmenergy.h @@ -23,10 +23,6 @@ namespace Io { class FileFormat; } -namespace QtGui { -class PythonScript; -} - namespace QtPlugins { class OBMMEnergy : public Avogadro::Calc::EnergyCalculator @@ -73,7 +69,7 @@ class OBMMEnergy : public Avogadro::Calc::EnergyCalculator /** * @brief Synchronous use of the QProcess. */ - class ProcessListener; + QByteArray writeAndRead(const QByteArray& input); private: Core::Molecule* m_molecule;