From 78696ba7563419a159daa0abfe163e663fe338bd Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sun, 22 Oct 2023 15:58:04 -0400 Subject: [PATCH] Add initial obmm support Signed-off-by: Geoff Hutchison --- avogadro/calc/energycalculator.cpp | 9 +- avogadro/qtgui/pythonscript.cpp | 35 +- avogadro/qtgui/pythonscript.h | 23 +- avogadro/qtplugins/forcefield/CMakeLists.txt | 3 +- avogadro/qtplugins/forcefield/forcefield.cpp | 97 +++-- avogadro/qtplugins/forcefield/forcefield.h | 2 +- avogadro/qtplugins/forcefield/obmmenergy.cpp | 366 ++++++++++++++++++ avogadro/qtplugins/forcefield/obmmenergy.h | 92 +++++ avogadro/qtplugins/forcefield/obmmprocess.cpp | 164 -------- avogadro/qtplugins/forcefield/obmmprocess.h | 163 -------- .../qtplugins/forcefield/scriptenergy.cpp | 65 +++- .../qtplugins/forcefield/scripts/gfnff.py | 26 +- avogadro/qtplugins/forcefield/scripts/uff.py | 25 +- 13 files changed, 654 insertions(+), 416 deletions(-) create mode 100644 avogadro/qtplugins/forcefield/obmmenergy.cpp create mode 100644 avogadro/qtplugins/forcefield/obmmenergy.h delete mode 100644 avogadro/qtplugins/forcefield/obmmprocess.cpp delete mode 100644 avogadro/qtplugins/forcefield/obmmprocess.h diff --git a/avogadro/calc/energycalculator.cpp b/avogadro/calc/energycalculator.cpp index 88eb83cec8..89222b569d 100644 --- a/avogadro/calc/energycalculator.cpp +++ b/avogadro/calc/energycalculator.cpp @@ -20,13 +20,18 @@ void EnergyCalculator::cleanGradients(TVector& grad) unsigned int size = grad.rows(); // check for overflows -- in case of divide by zero, etc. for (unsigned int i = 0; i < size; ++i) { - if (!std::isfinite(grad[i])) { + if (!std::isfinite(grad[i]) || std::isnan(grad[i])) { grad[i] = 0.0; } } // freeze any masked atoms or coordinates - grad = grad.cwiseProduct(m_mask); + /* + if (m_mask.rows() == size) + grad = grad.cwiseProduct(m_mask); + else + std::cerr << "Error: mask size " << m_mask.rows() << " " << grad.rows() << std::endl; + */ } } // namespace Avogadro diff --git a/avogadro/qtgui/pythonscript.cpp b/avogadro/qtgui/pythonscript.cpp index e5e37ab9b4..8b567ceb43 100644 --- a/avogadro/qtgui/pythonscript.cpp +++ b/avogadro/qtgui/pythonscript.cpp @@ -221,37 +221,32 @@ void PythonScript::processFinished(int, QProcess::ExitStatus) emit finished(); } -QByteArray PythonScript::asyncWriteAndResponse(QByteArray input, const int expectedLines) +void PythonScript::asyncTerminate() +{ + if (m_process != nullptr) { + m_process->terminate(); + disconnect(m_process, SIGNAL(finished()), this, SLOT(processFinished())); + m_process->deleteLater(); + m_process = nullptr; + } +} + +QByteArray PythonScript::asyncWriteAndResponse(QByteArray input) { if (m_process == nullptr) { return QByteArray(); // wait } - qDebug() << " writing to process: " << input << " " << expectedLines; - m_process->write(input); QByteArray buffer; - if (expectedLines == -1) { - m_process->waitForFinished(); - buffer = m_process->readAll(); - } else { - // only read a certain number of lines - // (e.g., energy and gradients) - int remainingLines = expectedLines; - bool ready = m_process->waitForReadyRead(); - if (ready) { - while (m_process->canReadLine() && remainingLines > 0) { - buffer += m_process->readLine(); - remainingLines--; - } - // clear the rest of the buffer - qDebug() << " remaining " << m_process->readAll(); + bool ready = m_process->waitForReadyRead(); + if (ready) { + while (m_process->canReadLine()) { + buffer += m_process->readLine(); } } - qDebug() << " read from process: " << buffer; - return buffer; } diff --git a/avogadro/qtgui/pythonscript.h b/avogadro/qtgui/pythonscript.h index b49491076e..4b19d5cdc8 100644 --- a/avogadro/qtgui/pythonscript.h +++ b/avogadro/qtgui/pythonscript.h @@ -12,9 +12,9 @@ #include #include +#include #include #include -#include namespace Avogadro { namespace QtGui { @@ -92,23 +92,26 @@ class AVOGADROQTGUI_EXPORT PythonScript : public QObject * Start a new process to execute asynchronously * " [args ...]", * optionally passing scriptStdin to the processes standard input. - * + * * Will send asyncFinished() signal when finished */ void asyncExecute(const QStringList& args, - const QByteArray& scriptStdin = QByteArray()); + const QByteArray& scriptStdin = QByteArray()); /** * Write input to the asynchronous process' standard input and return the * standard output when ready. Does not wait for the process to terminate * before returning (e.g. "server mode"). - * + * * @param input The input to write to the process' standard input - * @param expectedLines The number of lines to expect in the output. If -1, - * the output is read until the process terminates. * @return The standard output of the process + */ + QByteArray asyncWriteAndResponse(QByteArray input); + + /** + * Terminate the asynchronous process. */ - QByteArray asyncWriteAndResponse(QByteArray input, const int expectedLines = -1); + void asyncTerminate(); /** * Returns the standard output of the asynchronous process when finished. @@ -116,9 +119,9 @@ class AVOGADROQTGUI_EXPORT PythonScript : public QObject QByteArray asyncResponse(); signals: -/** - * The asynchronous execution is finished or timed out - */ + /** + * The asynchronous execution is finished or timed out + */ void finished(); public slots: diff --git a/avogadro/qtplugins/forcefield/CMakeLists.txt b/avogadro/qtplugins/forcefield/CMakeLists.txt index 115725018d..238d46ba9f 100644 --- a/avogadro/qtplugins/forcefield/CMakeLists.txt +++ b/avogadro/qtplugins/forcefield/CMakeLists.txt @@ -1,5 +1,6 @@ set(forcefield_srcs forcefield.cpp + obmmenergy.cpp scriptenergy.cpp ) @@ -18,8 +19,6 @@ set(forcefields scripts/gfn1.py scripts/gfn2.py scripts/gfnff.py - scripts/mmff94.py - scripts/uff.py ) install(PROGRAMS ${forcefields} diff --git a/avogadro/qtplugins/forcefield/forcefield.cpp b/avogadro/qtplugins/forcefield/forcefield.cpp index ba86fa6194..22ed9658c4 100644 --- a/avogadro/qtplugins/forcefield/forcefield.cpp +++ b/avogadro/qtplugins/forcefield/forcefield.cpp @@ -5,6 +5,7 @@ ******************************************************************************/ #include "forcefield.h" +#include "obmmenergy.h" #include "scriptenergy.h" #include @@ -50,7 +51,10 @@ const int constraintAction = 5; Forcefield::Forcefield(QObject* parent_) : ExtensionPlugin(parent_) { - refreshScripts(); + //refreshScripts(); + Calc::EnergyManager::registerModel(new OBMMEnergy("MMFF94")); + Calc::EnergyManager::registerModel(new OBMMEnergy("UFF")); + Calc::EnergyManager::registerModel(new OBMMEnergy("GAFF")); QAction* action = new QAction(this); action->setEnabled(true); @@ -119,7 +123,8 @@ void Forcefield::optimize() bool isInteractive = m_molecule->undoMolecule()->isInteractive(); m_molecule->undoMolecule()->setInteractive(true); - cppoptlib::LbfgsSolver solver; + //cppoptlib::LbfgsSolver solver; + cppoptlib::ConjugatedGradientDescentSolver solver; int n = m_molecule->atomCount(); // we have to cast the current 3d positions into a VectorXd @@ -127,6 +132,7 @@ void Forcefield::optimize() double* p = pos[0].data(); Eigen::Map map(p, 3 * n); Eigen::VectorXd positions = map; + Eigen::VectorXd lastPositions = positions; Eigen::VectorXd gradient = Eigen::VectorXd::Zero(3 * n); // just to get the right size / shape @@ -137,14 +143,14 @@ void Forcefield::optimize() cppoptlib::Criteria crit = cppoptlib::Criteria::defaults(); // e.g., every 5 steps, update coordinates - crit.iterations = m_nSteps; + crit.iterations = 5; // we don't set function or gradient criteria // .. these seem to be broken in the solver code // .. so we handle ourselves solver.setStopCriteria(crit); // set the method - std::string recommended = recommendedForceField(); + std::string recommended = "UFF"; qDebug() << "Energy method: " << recommended.c_str(); if (m_method == nullptr) { @@ -152,39 +158,78 @@ void Forcefield::optimize() m_method = Calc::EnergyManager::instance().model(recommended); } m_method->setMolecule(m_molecule); - m_method->setMask(m_molecule->frozenAtomMask()); + auto mask = m_molecule->frozenAtomMask(); + if (mask.rows() != m_molecule->atomCount()) { + mask = Eigen::VectorXd::Zero(3 * m_molecule->atomCount()); + // set to 1.0 + for (Index i = 0; i < 3 * m_molecule->atomCount(); ++i) { + mask[i] = 1.0; + } + } + m_method->setMask(mask); Real energy = m_method->value(positions); + m_method->gradient(positions, gradient); + qDebug() << " initial " << energy + << " gradNorm: " << gradient.norm() + << " posNorm: " << positions.norm(); + + Real currentEnergy = 0.0; for (unsigned int i = 0; i < m_maxSteps / crit.iterations; ++i) { solver.minimize(*m_method, positions); - Real currentEnergy = m_method->value(positions); + currentEnergy = m_method->value(positions); // get the current gradient for force visualization m_method->gradient(positions, gradient); + qDebug() << " optimize " << i << currentEnergy + << " gradNorm: " << gradient.norm() + << " posNorm: " << positions.norm(); // update coordinates - const double* d = positions.data(); - // casting back would be lovely... - for (size_t i = 0; i < n; ++i) { - pos[i] = Vector3(*(d), *(d + 1), *(d + 2)); - d += 3; - - forces[i] = -1.0 * Vector3(gradient[3 * i], gradient[3 * i + 1], - gradient[3 * i + 2]); + bool isFinite = std::isfinite(currentEnergy); + if (isFinite) { + const double* d = positions.data(); + bool isFinite = true; + // casting back would be lovely... + for (size_t i = 0; i < n; ++i) { + if (!std::isfinite(*d) || !std::isfinite(*(d + 1)) || + !std::isfinite(*(d + 2))) { + isFinite = false; + break; + } + + pos[i] = Vector3(*(d), *(d + 1), *(d + 2)); + d += 3; + + forces[i] = -1.0 * 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); - - // check for convergence - if (fabs(gradient.maxCoeff()) < m_gradientTolerance) - break; - if (fabs(currentEnergy - energy) < m_tolerance) - break; - energy = currentEnergy; + if (isFinite) { + qDebug() << " finite! "; + m_molecule->undoMolecule()->setAtomPositions3d(pos, + tr("Optimize Geometry")); + m_molecule->setForceVectors(forces); + Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified; + m_molecule->emitChanged(changes); + 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); + } } m_molecule->undoMolecule()->setInteractive(isInteractive); diff --git a/avogadro/qtplugins/forcefield/forcefield.h b/avogadro/qtplugins/forcefield/forcefield.h index 2388d24f95..3cf7f8fd14 100644 --- a/avogadro/qtplugins/forcefield/forcefield.h +++ b/avogadro/qtplugins/forcefield/forcefield.h @@ -80,7 +80,7 @@ private slots: // defaults Minimizer m_minimizer = LBFGS; - unsigned int m_maxSteps = 250; + unsigned int m_maxSteps = 100; unsigned int m_nSteps = 5; double m_tolerance = 1.0e-6; double m_gradientTolerance = 1.0e-4; diff --git a/avogadro/qtplugins/forcefield/obmmenergy.cpp b/avogadro/qtplugins/forcefield/obmmenergy.cpp new file mode 100644 index 0000000000..42af19a97e --- /dev/null +++ b/avogadro/qtplugins/forcefield/obmmenergy.cpp @@ -0,0 +1,366 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "obmmenergy.h" + +#include + +#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) + { + connect(m_process, SIGNAL(readyRead()), SLOT(readyRead())); + } + + bool waitForOutput(QByteArray output, int msTimeout = 2000) + { + if (!wait(msTimeout)) + return false; + + // success! + output = m_output; + 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), +#if defined(_WIN32) + m_executable("obmm.exe") +#else + m_executable("obmm") +#endif +{ + // eventually CJSON might be nice + m_inputFormat = new Io::CmlFormat; + + if (method == "UFF") { + m_description = tr("Universal Force Field"); + m_elements.reset(); + for (unsigned int i = 1; i < 102; ++i) + m_elements.set(i); + } else if (method == "GAFF") { + m_description = tr("Generalized Amber Force Field"); + + // H, C, N, O, F, P, S, Cl, Br, and I + m_elements.set(1); + m_elements.set(6); + m_elements.set(7); + m_elements.set(8); + m_elements.set(9); + m_elements.set(15); + m_elements.set(16); + m_elements.set(17); + m_elements.set(35); + m_elements.set(53); + } else if (method == "MMFF94") { + m_description = tr("Merck Molecular Force Field 94"); + m_elements.reset(); + + // H, C, N, O, F, Si, P, S, Cl, Br, and I + m_elements.set(1); + m_elements.set(6); + m_elements.set(7); + m_elements.set(8); + m_elements.set(9); + m_elements.set(14); + m_elements.set(15); + m_elements.set(16); + m_elements.set(17); + m_elements.set(35); + m_elements.set(53); + } +} + +OBMMEnergy::~OBMMEnergy() +{ + delete m_inputFormat; + delete m_process; +} + +void OBMMEnergy::setupProcess() +{ + if (m_process != nullptr) { + m_process->kill(); + delete m_process; + } + + m_process = new QProcess(); + + // Read the AVO_OBMM_EXECUTABLE env var to optionally override the + // executable used. + QByteArray obmmExec = qgetenv("AVO_OBMM_EXECUTABLE"); + if (!obmmExec.isEmpty()) { + m_executable = obmmExec; + } else { + // If not overridden, look for an obmm next to the executable. + QDir baseDir(QCoreApplication::applicationDirPath()); + if (!baseDir.absolutePath().startsWith("/usr/") && + QFileInfo(baseDir.absolutePath() + '/' + m_executable).exists()) { + m_executable = baseDir.absolutePath() + '/' + m_executable; + QProcessEnvironment env = QProcessEnvironment::systemEnvironment(); +#if defined(_WIN32) + env.insert("BABEL_DATADIR", + QCoreApplication::applicationDirPath() + "/data"); +#else + QDir dir(QCoreApplication::applicationDirPath() + "/../share/openbabel"); + QStringList filters; + filters << "3.*"; + QStringList dirs = dir.entryList(filters); + if (dirs.size() == 1) { + env.insert("BABEL_DATADIR", QCoreApplication::applicationDirPath() + + "/../share/openbabel/" + dirs[0]); + } else { + qDebug() << "Error, Open Babel data directory not found."; + } + dir.setPath(QCoreApplication::applicationDirPath() + "/../lib/openbabel"); + dirs = dir.entryList(filters); + if (dirs.size() == 1) { + env.insert("BABEL_LIBDIR", QCoreApplication::applicationDirPath() + + "/../lib/openbabel/" + dirs[0]); + } else { + qDebug() << "Error, Open Babel plugins directory not found."; + } +#endif + m_process->setProcessEnvironment(env); + } + } +} + +Calc::EnergyCalculator* OBMMEnergy::newInstance() const +{ + return new OBMMEnergy(m_name); +} + +void OBMMEnergy::setMolecule(Core::Molecule* mol) +{ + m_molecule = mol; + + // should check if the molecule is valid for this script + // .. this should never happen, but let's be defensive + if (mol == nullptr) { + return; // nothing to do + } + + setupProcess(); + + // start the process + // we need a tempory file to write the molecule + // get a temporary filename + QString tempPath = QDir::tempPath(); + if (!tempPath.endsWith(QDir::separator())) + tempPath += QDir::separator(); + QString tempPattern = tempPath + "avogadroOBMMXXXXXX.cml"; + m_tempFile.setFileTemplate(tempPattern); + if (!m_tempFile.open()) { + // appendError("Error creating temporary file."); + return; + } + + // write the molecule + m_inputFormat->writeFile(m_tempFile.fileName().toStdString(), *mol); + m_tempFile.close(); + + // start the process + m_process->start(m_executable, QStringList() << m_tempFile.fileName()); + if (!m_process->waitForStarted()) { + // appendError("Error starting process."); + return; + } + ProcessListener listener(m_process); + QByteArray result; + if (!listener.waitForOutput(result)) { + // appendError("Error running process."); + return; + } + + // 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); + if (!listener.waitForOutput(result)) { + // appendError("Error running process."); + return; + } + qDebug() << "OBMM energy: " << result; +} + +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; + + // write the new coordinates and read the energy + QByteArray 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->waitForBytesWritten(); + m_process->waitForReadyRead(); + result.clear(); + while (!result.contains("command >")) { + result += m_process->readLine(); + } + + qDebug() << " asking energy " << result << m_process->state(); + + // now ask for the energy + input = "energy\n\n"; + result.clear(); + m_process->write(input); + m_process->waitForBytesWritten(); + m_process->waitForReadyRead(); + while (!result.contains("command >")) { + result += m_process->readLine(); + } + + qDebug() << "OBMM: " << result; + + // 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 +} + +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"; + 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; + + // 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; + + // go through lines in result until we see "gradient " + QStringList lines = QString(result).remove('\r').split('\n'); + bool readingGradient = false; + unsigned int i = 0; + for (auto line : lines) { + if (line.contains("gradient")) { + readingGradient = true; + continue; + } + if (readingGradient) { + QStringList items = line.split(" ", Qt::SkipEmptyParts); + if (items.size() == 3) { + grad[3 * i] = items[0].toDouble(); + grad[3 * i + 1] = items[1].toDouble(); + grad[3 * i + 2] = items[2].toDouble(); + ++i; + } + } + } + + cleanGradients(grad); +} + +} // namespace Avogadro::QtPlugins + +#include "obmmenergy.moc" diff --git a/avogadro/qtplugins/forcefield/obmmenergy.h b/avogadro/qtplugins/forcefield/obmmenergy.h new file mode 100644 index 0000000000..1c6bf2e902 --- /dev/null +++ b/avogadro/qtplugins/forcefield/obmmenergy.h @@ -0,0 +1,92 @@ +/****************************************************************************** + 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_OBMMENERGY_H +#define AVOGADRO_QTPLUGINS_OBMMENERGY_H + +#include + +#include + +#include +#include +#include +#include + +class QJsonObject; + +namespace Avogadro { + +namespace Io { +class FileFormat; +} + +namespace QtGui { +class PythonScript; +} + +namespace QtPlugins { + +class OBMMEnergy : public Avogadro::Calc::EnergyCalculator +{ + Q_DECLARE_TR_FUNCTIONS(OBMMEnergy) + +public: + /** Formats that may be written to the script's input. */ + enum Format + { + NotUsed, + Cjson, + Cml, + Mdl, // sdf + Pdb, + Xyz + }; + + OBMMEnergy(const std::string& method = ""); + ~OBMMEnergy() override; + + std::string method() const { return m_identifier; } + void setupProcess(); + + Calc::EnergyCalculator* newInstance() const override; + + std::string identifier() const override { return m_identifier; } + std::string name() const override { return m_name; } + std::string description() const override { return m_description.toStdString(); } + + Core::Molecule::ElementMask elements() const override { return (m_elements); } + + // This will check if the molecule is valid for this script + // and then start the external process + void setMolecule(Core::Molecule* mol) override; + // energy + Real value(const Eigen::VectorXd& x) override; + // gradient (which may be unsupported and fall back to numeric) + void gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) override; + + /** + * @brief Synchronous use of the QProcess. + */ + class ProcessListener; + +private: + Core::Molecule* m_molecule; + Io::FileFormat* m_inputFormat; + QProcess* m_process; + QString m_executable; + + Core::Molecule::ElementMask m_elements; + std::string m_identifier; + std::string m_name; + QString m_description; + + QTemporaryFile m_tempFile; +}; + +} // namespace QtPlugins +} // namespace Avogadro + +#endif // AVOGADRO_QTPLUGINS_OBMMENERGY_H diff --git a/avogadro/qtplugins/forcefield/obmmprocess.cpp b/avogadro/qtplugins/forcefield/obmmprocess.cpp deleted file mode 100644 index 5fae2b990f..0000000000 --- a/avogadro/qtplugins/forcefield/obmmprocess.cpp +++ /dev/null @@ -1,164 +0,0 @@ -/****************************************************************************** - - This source file is part of the Avogadro project. - - Copyright 2019 Geoffrey R. Hutchison - - This source code is released under the New BSD License, (the "License"). - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. - -******************************************************************************/ - -#include "OBMMProcess.h" - -#include -#include -#include -#include -#include -#include - -namespace Avogadro { -namespace QtPlugins { - -OBMMProcess::OBMMProcess(QObject* parent_) - : QObject(parent_), m_processLocked(false), m_aborted(false), - m_process(new QProcess(this)), -#if defined(_WIN32) - m_obmmExecutable("obmm.exe") -#else - m_obmmExecutable("obmm") -#endif -{ - // Read the AVO_OBABEL_EXECUTABLE env var to optionally override the - // executable used. - QByteArray obmmExec = qgetenv("AVO_OBMM_EXECUTABLE"); - if (!obabelExec.isEmpty()) { - m_obmmExecutable = obmmExec; - } else { - // If not overridden, look for an obabel next to the executable. - QDir baseDir(QCoreApplication::applicationDirPath()); - if (!baseDir.absolutePath().startsWith("/usr/") && - QFileInfo(baseDir.absolutePath() + '/' + m_obabelExecutable).exists()) { - m_obabelExecutable = baseDir.absolutePath() + '/' + m_obmmExecutable; - QProcessEnvironment env = QProcessEnvironment::systemEnvironment(); -#if defined(_WIN32) - env.insert("BABEL_DATADIR", - QCoreApplication::applicationDirPath() + "/data"); -#else - QDir dir(QCoreApplication::applicationDirPath() + "/../share/openbabel"); - QStringList filters; - filters << "2.*"; - QStringList dirs = dir.entryList(filters); - if (dirs.size() == 1) { - env.insert("BABEL_DATADIR", QCoreApplication::applicationDirPath() + - "/../share/openbabel/" + dirs[0]); - } else { - qDebug() << "Error, Open Babel data directory not found."; - } - dir.setPath(QCoreApplication::applicationDirPath() + "/../lib/openbabel"); - dirs = dir.entryList(filters); - if (dirs.size() == 1) { - env.insert("BABEL_LIBDIR", QCoreApplication::applicationDirPath() + - "/../lib/openbabel/" + dirs[0]); - } else { - qDebug() << "Error, Open Babel plugins directory not found."; - } -#endif - m_process->setProcessEnvironment(env); - } - } -} - -void OBMMProcess::abort() -{ - m_aborted = true; - emit aborted(); -} - -void OBMMProcess::obError() -{ - qDebug() << "Process encountered an error, and did not execute correctly."; - if (m_process) { - qDebug() << "\tExit code:" << m_process->exitCode(); - qDebug() << "\tExit status:" << m_process->exitStatus(); - qDebug() << "\tExit output:" << m_process->readAll(); - } -} - -bool OBMMProcess::queryForceFields() -{ - if (!tryLockProcess()) { - qWarning() << "OBMMProcess::queryForceFields(): process already in use."; - return false; - } - - QStringList options; - options << "-L" - << "forcefields"; - - executeObabel(options, this, SLOT(queryForceFieldsPrepare())); - return true; -} - -void OBMMProcess::queryForceFieldsPrepare() -{ - if (m_aborted) { - releaseProcess(); - return; - } - - QMap result; - - QString output = QString::fromLatin1(m_process->readAllStandardOutput()); - - QRegExp parser("([^\\s]+)\\s+(\\S[^\\n]*[^\\n\\.]+)\\.?\\n"); - int pos = 0; - while ((pos = parser.indexIn(output, pos)) != -1) { - QString key = parser.cap(1); - QString desc = parser.cap(2); - result.insertMulti(key, desc); - pos += parser.matchedLength(); - } - - releaseProcess(); - emit queryForceFieldsFinished(result); -} - -void OBMMProcess::executeObabel(const QStringList& options, QObject* receiver, - const char* slot, const QByteArray& obabelStdin) -{ - // Setup exit handler - if (receiver) { - connect(m_process, SIGNAL(finished(int)), receiver, slot); - connect(m_process, SIGNAL(error(QProcess::ProcessError)), receiver, slot); - connect(m_process, SIGNAL(error(QProcess::ProcessError)), this, - SLOT(obError())); - } - - // Start process - qDebug() << "OBMMProcess::executeObabel: " - "Running" - << m_obabelExecutable << options.join(" "); - m_process->start(m_obabelExecutable, options); - if (!obabelStdin.isNull()) { - m_process->write(obabelStdin); - m_process->closeWriteChannel(); - } -} - -void OBMMProcess::resetState() -{ - m_aborted = false; - m_process->disconnect(this); - disconnect(m_process); - connect(this, SIGNAL(aborted()), m_process, SLOT(kill())); -} - -} // namespace QtPlugins -} // namespace Avogadro diff --git a/avogadro/qtplugins/forcefield/obmmprocess.h b/avogadro/qtplugins/forcefield/obmmprocess.h deleted file mode 100644 index 87e026d5f3..0000000000 --- a/avogadro/qtplugins/forcefield/obmmprocess.h +++ /dev/null @@ -1,163 +0,0 @@ -/****************************************************************************** - - This source file is part of the Avogadro project. - - Copyright 2019 Geoffrey R. Hutchison - - This source code is released under the New BSD License, (the "License"). - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. - -******************************************************************************/ - -#ifndef AVOGADRO_QTPLUGINS_OBMMPROCESS_H -#define AVOGADRO_QTPLUGINS_OBMMPROCESS_H - -#include -#include -#include - -class QProcess; - -namespace Avogadro { -namespace QtPlugins { - -/** - * @brief The OBMMProcess class provides an interface to the `obmm` executable, - * which is run in a separate process. - * - * The `obmm` executable used by this class can be overridden by setting the - * AVO_OBABEL_EXECUTABLE environment variable. - */ -class OBMMProcess : public QObject -{ - Q_OBJECT -public: - explicit OBMMProcess(QObject* parent_ = 0); - - /** - * @name Process Management - * Methods, slots, and signals used to interact with the OpenBabel process. - * @{ - */ -public: - /** - * The `obmm` executable used by the process. - */ - QString obmmExecutable() const { return m_obmmExecutable; } - - /** - * @return True if the process is in use, false otherwise. - */ - bool inUse() const { return m_processLocked; } - -public slots: - /** - * Abort any currently running processes. - * - * This will cause aborted() to be emitted, but not any of the - * operation-specific "finished" signals. - */ - void abort(); - - /** - * Called when an error in the process occurs. - */ - void obError(); - -signals: - /** - * Emitted when the abort() method has been called. - */ - void aborted(); - - // end Process Management doxygen group - /**@}*/ - -public slots: - /** - * Request a list of all supported force fields from obabel. - * - * After calling this method, the queryForceFieldsFinished signal will be - * emitted. This method executes - * - * `obabel -L forcefields` - * - * and parses the output. - * - * If an error occurs, queryReadFormatsFinished will be emitted with an empty - * argument. - * - * @return True if the process started successfully, false otherwise. - */ - bool queryForceFields(); - -signals: - /** - * Triggered when the process started by queryForceFields() completes. - * @param forceFields The force fields supported by OpenBabel. Keys - * are unique identifiers for the force fields, and the values are - * non-translated (english), human-readable descriptions. - * - * If an error occurs, forceFields will be empty. - */ - void queryForceFieldsFinished(const QMap& forceFields); - -private slots: - void queryForceFieldsPrepare(); - - -private: - /** - * Internal method for launching the obmm executable. - * @param options List of options to pass to QProcess::start - * @param receiver A QObject subclass instance that has @a slot as a member. - * @param slot The slot to call when completed. Must have no arguments. - * @param obmmStdin Standard input for the obmm process (optional). - * - * Call this method like so: -@code -QStringList options; - -executeobmm(options, this, SLOT(mySlot())); -@endcode - * - * @a slot will be connected to QProcess::finished(int) and - * QProcess::error(QProcess::ProcessError) with @a receiver as receiver and - * @a m_process as sender. @a m_process is then started using - * m_obmmExecutable and options as arguments. If provided, the obmmStdin - * data will be written to the obmm stdin channel. - */ - void executeobmm(const QStringList& options, QObject* receiver = nullptr, - const char* slot = nullptr, - const QByteArray& obmmStdin = QByteArray()); - - void resetState(); - - // Not thread safe -- just uses a bool. - bool tryLockProcess() - { - if (m_processLocked) - return false; - m_processLocked = true; - resetState(); - return true; - } - - void releaseProcess() { m_processLocked = false; } - - bool m_processLocked; - bool m_aborted; - QProcess* m_process; - QString m_obmmExecutable; - -}; - -} // namespace QtPlugins -} // namespace Avogadro - -#endif // AVOGADRO_QTPLUGINS_OBMMProcess_H diff --git a/avogadro/qtplugins/forcefield/scriptenergy.cpp b/avogadro/qtplugins/forcefield/scriptenergy.cpp index d5e09dbfda..939d6ca1b8 100644 --- a/avogadro/qtplugins/forcefield/scriptenergy.cpp +++ b/avogadro/qtplugins/forcefield/scriptenergy.cpp @@ -100,9 +100,9 @@ void ScriptEnergy::setMolecule(Core::Molecule* mol) QStringList options; options << "-f" << m_tempFile.fileName(); + // if there was a previous process, kill it + m_interpreter->asyncTerminate(); // start the interpreter - //@ todo - check if there was a previous process running - // .. if so, kill it m_interpreter->asyncExecute(options); } @@ -118,15 +118,66 @@ Real ScriptEnergy::value(const Eigen::VectorXd& x) input += QString::number(x[i]) + " " + QString::number(x[i + 1]) + " " + QString::number(x[i + 2]) + "\n"; } - QByteArray result = m_interpreter->asyncWriteAndResponse(input, 1); + QByteArray result = m_interpreter->asyncWriteAndResponse(input); + + // go through lines in result until we see "AvogadroEnergy: " + QStringList lines = QString(result).remove('\r').split('\n'); + double energy = 0.0; + for (auto line : lines) { + if (line.startsWith("AvogadroEnergy:")) { + QStringList items = line.split(" "); + if (items.size() > 1) + energy = items[1].toDouble(); + } + } - return result.toDouble(); // if conversion fails, returns 0.0 + return energy; // if conversion fails, returns 0.0 } void ScriptEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) { - //@ todo read the gradient if available from the process - EnergyCalculator::gradient(x, grad); + if (!m_gradients) { + EnergyCalculator::gradient(x, grad); + return; + } + + // Get the gradient from the script + // write the new coordinates and read the energy + QByteArray input; + 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"; + } + QByteArray result = m_interpreter->asyncWriteAndResponse(input); + + // parse the result + // first split on newlines + QStringList lines = QString(result).remove('\r').split('\n'); + double energy = 0.0; + unsigned int i = 0; + bool readingGrad = false; + for (auto line : lines) { + if (line.startsWith("AvogadroGradient:")) { + readingGrad = true; + continue; // next line + } + + if (readingGrad) { + QStringList items = line.split(" "); + if (items.size() == 3) { + grad[i] = items[0].toDouble(); + grad[i + 1] = items[1].toDouble(); + grad[i + 2] = items[2].toDouble(); + i += 3; + } + + if (i > x.size()) + break; + } + } + + cleanGradients(grad); } ScriptEnergy::Format ScriptEnergy::stringToFormat(const std::string& str) @@ -293,8 +344,6 @@ void ScriptEnergy::readMetaData() // get the element mask // (if it doesn't exist, the default is no elements anyway) m_valid = parseElements(metaData); - - qDebug() << " finished parsing metadata "; } bool ScriptEnergy::parseString(const QJsonObject& ob, const QString& key, diff --git a/avogadro/qtplugins/forcefield/scripts/gfnff.py b/avogadro/qtplugins/forcefield/scripts/gfnff.py index c415ec5eb3..24b2b3876f 100644 --- a/avogadro/qtplugins/forcefield/scripts/gfnff.py +++ b/avogadro/qtplugins/forcefield/scripts/gfnff.py @@ -67,18 +67,6 @@ def run(filename): # we loop forever - Avogadro will kill our process when done while True: - # first print the energy of these coordinates - print(res.get_energy()) # in Hartree - - with open("/Users/ghutchis/gfnff.log", "a") as f: - f.write(str(res.get_energy()) + "\n") - - # now print the gradient - # .. we don't want the "[]" in the output - output = np.array2string(res.get_gradient()) - output = output.replace("[", "").replace("]", "") - print(output) - # read new coordinates from stdin for i in range(len(atoms)): coordinates[i] = np.fromstring(input(), sep=" ") @@ -89,6 +77,20 @@ def run(filename): calc.update(coordinates) calc.singlepoint(res) + # first print the energy of these coordinates + print("AvogadroEnergy:", res.get_energy()) # in Hartree + + with open("/Users/ghutchis/gfnff.log", "a") as f: + f.write(str(res.get_energy()) + "\n") + + # now print the gradient + # .. we don't want the "[]" in the output + print("AvogadroGradient:") + grad = res.get_gradient() * 4961.475 + output = np.array2string(grad) + output = output.replace("[", "").replace("]", "") + print(output) + if __name__ == "__main__": parser = argparse.ArgumentParser("GFN2 calculator") diff --git a/avogadro/qtplugins/forcefield/scripts/uff.py b/avogadro/qtplugins/forcefield/scripts/uff.py index 7d35226262..a2d7945902 100644 --- a/avogadro/qtplugins/forcefield/scripts/uff.py +++ b/avogadro/qtplugins/forcefield/scripts/uff.py @@ -47,23 +47,32 @@ def run(filename): # we loop forever - Avogadro will kill the process when done num_atoms = len(mol.atoms) while True: - # first print the energy of these coordinates - print(ff.Energy(True)) # in Hartree - - # now print the gradient on each atom - for atom in mol.atoms: - grad = ff.GetGradient(atom.OBAtom) - print(grad.GetX(), grad.GetY(), grad.GetZ()) - # read new coordinates from stdin for i in range(num_atoms): coordinates = np.fromstring(input(), sep=" ") atom = mol.atoms[i] atom.OBAtom.SetVector(coordinates[0], coordinates[1], coordinates[2]) + with open("/Users/ghutchis/uff.log", "a") as f: + f.write("Coordinates:\n") + for atom in mol.atoms: + f.write(str(atom.OBAtom.GetX()) + " " + str(atom.OBAtom.GetY()) + " " + str(atom.OBAtom.GetZ()) + "\n") + # update the molecule geometry for the next energy ff.SetCoordinates(mol.OBMol) + # first print the energy of these coordinates + energy = ff.Energy(True) # in Hartree + print("AvogadroEnergy:", energy) # in Hartree + with open("/Users/ghutchis/uff.log", "a") as f: + f.write(str(energy) + "\n") + + # now print the gradient on each atom + print("AvogadroGradient:") + for atom in mol.atoms: + grad = ff.GetGradient(atom.OBAtom) + print(grad.GetX(), grad.GetY(), grad.GetZ()) + if __name__ == "__main__": parser = argparse.ArgumentParser("UFF calculator")