diff --git a/avogadro/calc/energycalculator.cpp b/avogadro/calc/energycalculator.cpp index 9eac17dd5e..25fbcda18a 100644 --- a/avogadro/calc/energycalculator.cpp +++ b/avogadro/calc/energycalculator.cpp @@ -27,22 +27,4 @@ void EnergyCalculator::cleanGradients(TVector& grad) grad = grad.cwiseProduct(m_mask); } -void EnergyCalculator::freezeAtom(Index atomId) -{ - if (atomId * 3 <= m_mask.rows() - 3) { - m_mask[atomId*3] = 0.0; - m_mask[atomId*3+1] = 0.0; - m_mask[atomId*3+2] = 0.0; - } -} - -void EnergyCalculator::unfreezeAtom(Index atomId) -{ - if (atomId * 3 <= m_mask.rows() - 3) { - m_mask[atomId*3] = 1.0; - m_mask[atomId*3+1] = 1.0; - m_mask[atomId*3+2] = 1.0; - } -} - } // namespace Avogadro diff --git a/avogadro/calc/energycalculator.h b/avogadro/calc/energycalculator.h index f59c73546b..b9784e4055 100644 --- a/avogadro/calc/energycalculator.h +++ b/avogadro/calc/energycalculator.h @@ -27,6 +27,12 @@ class AVOGADROCALC_EXPORT EnergyCalculator : public cppoptlib::Problem EnergyCalculator() {} ~EnergyCalculator() {} + /** + * Create a new instance of the model. Ownership passes to the + * caller. + */ + virtual EnergyCalculator* newInstance() const = 0; + /** * @return a unique identifier for this calculator. */ @@ -96,16 +102,23 @@ class AVOGADROCALC_EXPORT EnergyCalculator : public cppoptlib::Problem */ TVector mask() const { return m_mask; } - void freezeAtom(Index atomId); - void unfreezeAtom(Index atomId); - /** * Called when the current molecule changes. */ virtual void setMolecule(Core::Molecule* mol) = 0; protected: + /** + * @brief Append an error to the error string for the model. + * @param errorString The error to be added. + * @param newLine Add a new line after the error string? + */ + void appendError(const std::string& errorString, bool newLine = true) const; + TVector m_mask; // optimize or frozen atom mask + +private: + mutable std::string m_error; }; } // end namespace Calc diff --git a/avogadro/calc/energymanager.cpp b/avogadro/calc/energymanager.cpp index fdd23cd6dd..76e93d78f0 100644 --- a/avogadro/calc/energymanager.cpp +++ b/avogadro/calc/energymanager.cpp @@ -81,7 +81,11 @@ std::string EnergyManager::nameForModel(const std::string& identifier) const EnergyManager::EnergyManager() { - // add any default models here (EEM maybe?) + // add any default models here + + // LJ is the fallback, since it can handle anything + // (maybe not well, but it can handle it) + addModel(new LennardJones); } EnergyManager::~EnergyManager() diff --git a/avogadro/calc/lennardjones.h b/avogadro/calc/lennardjones.h index 5fb7466644..3334ffe48d 100644 --- a/avogadro/calc/lennardjones.h +++ b/avogadro/calc/lennardjones.h @@ -14,7 +14,7 @@ namespace Avogadro { namespace Core { class Molecule; class UnitCell; -} +} // namespace Core namespace Calc { @@ -24,25 +24,23 @@ class AVOGADROCALC_EXPORT LennardJones : public EnergyCalculator LennardJones(); ~LennardJones(); - std::string identifier() const override - { return "LJ"; } + LennardJones* newInstance() const override { return new LennardJones; } - std::string name() const override - { return "Lennard-Jones"; } + std::string identifier() const override { return "LJ"; } + + std::string name() const override { return "Lennard-Jones"; } std::string description() const override - { return "Universal Lennard-Jones potential"; } + { + return "Universal Lennard-Jones potential"; + } bool acceptsUnitCell() const override { return true; } - Core::Molecule::ElementMask elements() const override - { - return (m_elements); - } + Core::Molecule::ElementMask elements() const override { return (m_elements); } Real value(const Eigen::VectorXd& x) override; - void gradient(const Eigen::VectorXd& x, - Eigen::VectorXd& grad) override; + void gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) override; /** * Called when the current molecule changes. diff --git a/avogadro/core/molecule.cpp b/avogadro/core/molecule.cpp index e04046c925..db6567f0a9 100644 --- a/avogadro/core/molecule.cpp +++ b/avogadro/core/molecule.cpp @@ -49,6 +49,7 @@ Molecule::Molecule(const Molecule& other) m_residues(other.m_residues), m_hallNumber(other.m_hallNumber), m_graph(other.m_graph), m_bondOrders(other.m_bondOrders), m_atomicNumbers(other.m_atomicNumbers), + m_frozenAtomMask(other.m_frozenAtomMask), m_layers(LayerManager::getMoleculeLayer(this)) { // Copy over any meshes @@ -90,6 +91,7 @@ Molecule::Molecule(Molecule&& other) noexcept m_residues(other.m_residues), m_hallNumber(other.m_hallNumber), m_graph(other.m_graph), m_bondOrders(other.m_bondOrders), m_atomicNumbers(other.m_atomicNumbers), + m_frozenAtomMask(other.m_frozenAtomMask), m_layers(LayerManager::getMoleculeLayer(this)) { m_basisSet = other.m_basisSet; @@ -133,6 +135,7 @@ Molecule& Molecule::operator=(const Molecule& other) m_bondOrders = other.m_bondOrders; m_atomicNumbers = other.m_atomicNumbers; m_hallNumber = other.m_hallNumber; + m_frozenAtomMask = other.m_frozenAtomMask; clearMeshes(); @@ -193,6 +196,7 @@ Molecule& Molecule::operator=(Molecule&& other) noexcept m_bondOrders = other.m_bondOrders; m_atomicNumbers = other.m_atomicNumbers; m_hallNumber = other.m_hallNumber; + m_frozenAtomMask = other.m_frozenAtomMask; clearMeshes(); m_meshes = std::move(other.m_meshes); @@ -267,6 +271,59 @@ std::set Molecule::partialChargeTypes() const return types; } +void Molecule::setFrozenAtom(Index atomId, bool frozen) +{ + if (atomId >= m_atomicNumbers.size()) + return; + + // check if we need to resize + unsigned int size = m_frozenAtomMask.rows(); + if (m_frozenAtomMask.rows() != 3*m_atomicNumbers.size()) + m_frozenAtomMask.conservativeResize(3*m_atomicNumbers.size()); + + // do we need to initialize new values? + if (m_frozenAtomMask.rows() > size) + for (unsigned int i = size; i < m_frozenAtomMask.rows(); ++i) + m_frozenAtomMask[i] = 1.0; + + float value = frozen ? 0.0 : 1.0; + if (atomId * 3 <= m_frozenAtomMask.rows() - 3) { + m_frozenAtomMask[atomId*3] = value; + m_frozenAtomMask[atomId*3+1] = value; + m_frozenAtomMask[atomId*3+2] = value; + } +} + +bool Molecule::frozenAtom(Index atomId) const +{ + bool frozen = false; + if (atomId * 3 <= m_frozenAtomMask.rows() - 3) { + frozen = (m_frozenAtomMask[atomId*3] == 0.0 && + m_frozenAtomMask[atomId*3+1] == 0.0 && + m_frozenAtomMask[atomId*3+2] == 0.0); + } + return frozen; +} + +void Molecule::setFrozenAtomAxis(Index atomId, int axis, bool frozen) +{ + // check if we need to resize + unsigned int size = m_frozenAtomMask.rows(); + if (m_frozenAtomMask.rows() != 3*m_atomicNumbers.size()) + m_frozenAtomMask.conservativeResize(3*m_atomicNumbers.size()); + + // do we need to initialize new values? + if (m_frozenAtomMask.rows() > size) + for (unsigned int i = size; i < m_frozenAtomMask.rows(); ++i) + m_frozenAtomMask[i] = 1.0; + + float value = frozen ? 0.0 : 1.0; + if (atomId * 3 <= m_frozenAtomMask.rows() - 3) { + m_frozenAtomMask[atomId*3 + axis] = value; + } +} + + void Molecule::setData(const std::string& name, const Variant& value) { m_data.setValue(name, value); diff --git a/avogadro/core/molecule.h b/avogadro/core/molecule.h index f2b921ebea..b4c61297d3 100644 --- a/avogadro/core/molecule.h +++ b/avogadro/core/molecule.h @@ -693,6 +693,26 @@ class AVOGADROCORE_EXPORT Molecule */ bool setAtomicNumber(Index atomId, unsigned char atomicNumber); + /** + * Freeze or unfreeze an atom for optimization + */ + void setFrozenAtom(Index atomId, bool frozen); + + /** + * Get the frozen status of an atom + */ + bool frozenAtom(Index atomId) const; + + /** + * Freeze or unfreeze X, Y, or Z coordinate of an atom for optimization + * @param atomId The index of the atom to modify. + * @param axis The axis to freeze (0, 1, or 2 for X, Y, or Z) + * @param frozen True to freeze, false to unfreeze + */ + void setFrozenAtomAxis(Index atomId, int axis, bool frozen); + + Eigen::VectorXd frozenAtomMask() const { return m_frozenAtomMask; } + /** * @return a map of components and count. */ @@ -763,6 +783,8 @@ class AVOGADROCORE_EXPORT Molecule // This will be stored from the last space group operation unsigned short m_hallNumber = 0; + Eigen::VectorXd m_frozenAtomMask; + private: mutable Graph m_graph; // A transformation of the molecule to a graph. // edge information diff --git a/avogadro/qtgui/pythonscript.cpp b/avogadro/qtgui/pythonscript.cpp index b075b91b75..04650e1704 100644 --- a/avogadro/qtgui/pythonscript.cpp +++ b/avogadro/qtgui/pythonscript.cpp @@ -221,6 +221,34 @@ void PythonScript::processFinished(int, QProcess::ExitStatus) emit finished(); } +QByteArray PythonScript::asyncWriteAndResponse(QByteArray input, const int expectedLines) +{ + if (m_process == nullptr) { + return QByteArray(); // wait + } + + 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--; + } + } + } + + return buffer; +} + QByteArray PythonScript::asyncResponse() { if (m_process == nullptr || m_process->state() == QProcess::Running) { diff --git a/avogadro/qtgui/pythonscript.h b/avogadro/qtgui/pythonscript.h index 28a259ce3e..b49491076e 100644 --- a/avogadro/qtgui/pythonscript.h +++ b/avogadro/qtgui/pythonscript.h @@ -98,6 +98,18 @@ class AVOGADROQTGUI_EXPORT PythonScript : public QObject void asyncExecute(const QStringList& args, 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, const int expectedLines = -1); + /** * Returns the standard output of the asynchronous process when finished. */ diff --git a/avogadro/qtplugins/forcefield/CMakeLists.txt b/avogadro/qtplugins/forcefield/CMakeLists.txt index 74fd65392d..115725018d 100644 --- a/avogadro/qtplugins/forcefield/CMakeLists.txt +++ b/avogadro/qtplugins/forcefield/CMakeLists.txt @@ -1,5 +1,6 @@ set(forcefield_srcs forcefield.cpp + scriptenergy.cpp ) avogadro_plugin(Forcefield @@ -22,4 +23,4 @@ set(forcefields ) install(PROGRAMS ${forcefields} -DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/forcefields/") +DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/energy/") diff --git a/avogadro/qtplugins/forcefield/forcefield.cpp b/avogadro/qtplugins/forcefield/forcefield.cpp index ca20d539cf..ab63b9a831 100644 --- a/avogadro/qtplugins/forcefield/forcefield.cpp +++ b/avogadro/qtplugins/forcefield/forcefield.cpp @@ -5,25 +5,25 @@ ******************************************************************************/ #include "forcefield.h" +#include "scriptenergy.h" #include #include #include -#include -#include #include #include +#include +#include #include -#include -#include -#include -#include #include #include #include +#include + +#include #include #include @@ -37,22 +37,16 @@ namespace Avogadro { namespace QtPlugins { +using Avogadro::Calc::EnergyCalculator; using Avogadro::QtGui::Molecule; using Avogadro::QtGui::RWMolecule; -using Avogadro::Calc::EnergyCalculator; const int energyAction = 0; const int optimizeAction = 1; const int configureAction = 2; const int freezeAction = 3; -Forcefield::Forcefield(QObject* parent_) - : ExtensionPlugin(parent_) - , m_molecule(nullptr) - , m_minimizer(LBFGS) - , m_method(0) // just Lennard-Jones for now - , m_maxSteps(100) - , m_outputFormat(nullptr) +Forcefield::Forcefield(QObject* parent_) : ExtensionPlugin(parent_) { refreshScripts(); @@ -74,7 +68,7 @@ Forcefield::Forcefield(QObject* parent_) action = new QAction(this); action->setEnabled(true); - action->setText(tr("Freeze Atoms")); // calculate energy + action->setText(tr("Freeze Atoms")); action->setData(freezeAction); connect(action, SIGNAL(triggered()), SLOT(freezeSelected())); m_actions.push_back(action); @@ -106,12 +100,7 @@ void Forcefield::setMolecule(QtGui::Molecule* mol) m_molecule = mol; - // @todo set molecule for a calculator -} - -void Forcefield::refreshScripts() -{ - // call the script loader + // @todo set the available method list } void Forcefield::optimize() @@ -123,9 +112,7 @@ void Forcefield::optimize() bool isInteractive = m_molecule->undoMolecule()->isInteractive(); m_molecule->undoMolecule()->setInteractive(true); - //@todo check m_minimizer for method to use - //cppoptlib::LbfgsSolver solver; - cppoptlib::ConjugatedGradientDescentSolver solver; + cppoptlib::LbfgsSolver solver; int n = m_molecule->atomCount(); Core::Array pos = m_molecule->atomPositions3d(); @@ -138,42 +125,49 @@ void Forcefield::optimize() // Create a Criteria class to adjust stopping criteria cppoptlib::Criteria crit = cppoptlib::Criteria::defaults(); - // @todo allow criteria to be set - crit.iterations = 5; - crit.fDelta = 1.0e-6; - solver.setStopCriteria(crit); // every 5 steps, update coordinates - cppoptlib::Status status = cppoptlib::Status::NotStarted; + + // e.g., every 5 steps, update coordinates + crit.iterations = m_nSteps; + // 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 //@todo check m_method for a particular calculator Calc::LennardJones lj; lj.setMolecule(m_molecule); + double energy = lj.value(positions); for (unsigned int i = 0; i < m_maxSteps / crit.iterations; ++i) { solver.minimize(lj, positions); - cppoptlib::Status currentStatus = solver.status(); - // qDebug() << " status: " << (int)currentStatus; - // qDebug() << " energy: " << lj.value(positions); - // check the gradient norm + double currentEnergy = lj.value(positions); 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); + + // 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] = -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; } m_molecule->undoMolecule()->setInteractive(isInteractive); @@ -212,5 +206,47 @@ void Forcefield::freezeSelected() } } -} // end QtPlugins +void Forcefield::refreshScripts() +{ + unregisterScripts(); + qDeleteAll(m_scripts); + m_scripts.clear(); + + QMap scriptPaths = + QtGui::ScriptLoader::scriptList("energy"); + foreach (const QString& filePath, scriptPaths) { + auto* model = new ScriptEnergy(filePath); + if (model->isValid()) + m_scripts.push_back(model); + else + delete model; + } + + registerScripts(); +} + +void Forcefield::unregisterScripts() +{ + for (QList::const_iterator + it = m_scripts.constBegin(), + itEnd = m_scripts.constEnd(); + it != itEnd; ++it) { + Calc::EnergyManager::unregisterModel((*it)->identifier()); + } } + +void Forcefield::registerScripts() +{ + for (QList::const_iterator + it = m_scripts.constBegin(), + itEnd = m_scripts.constEnd(); + it != itEnd; ++it) { + if (!Calc::EnergyManager::registerModel((*it)->newInstance())) { + qDebug() << "Could not register model" << (*it)->identifier().c_str() + << "due to name conflict."; + } + } +} + +} // namespace QtPlugins +} // namespace Avogadro diff --git a/avogadro/qtplugins/forcefield/forcefield.h b/avogadro/qtplugins/forcefield/forcefield.h index a96fffaa99..2c51afe531 100644 --- a/avogadro/qtplugins/forcefield/forcefield.h +++ b/avogadro/qtplugins/forcefield/forcefield.h @@ -16,6 +16,11 @@ class QAction; class QDialog; namespace Avogadro { + +namespace Calc { +class EnergyCalculator; +} + namespace QtPlugins { /** @@ -57,6 +62,8 @@ public slots: * Scan for new scripts in the Forcefield directories. */ void refreshScripts(); + void registerScripts(); + void unregisterScripts(); private slots: void energy(); @@ -65,18 +72,19 @@ private slots: private: QList m_actions; - QtGui::Molecule* m_molecule; + QtGui::Molecule* m_molecule = nullptr; - Minimizer m_minimizer; - unsigned int m_method; - unsigned int m_maxSteps; + // defaults + Minimizer m_minimizer = LBFGS; + unsigned int m_maxSteps = 250; + unsigned int m_nSteps = 5; + double m_tolerance = 1.0e-6; + double m_gradientTolerance = 1.0e-4; + Calc::EnergyCalculator *m_method = nullptr; - // maps program name --> script file path - QMap m_forcefieldScripts; - - const Io::FileFormat* m_outputFormat; - QString m_tempFileName; + QList m_scripts; }; + } } diff --git a/avogadro/qtplugins/forcefield/scriptenergy.cpp b/avogadro/qtplugins/forcefield/scriptenergy.cpp new file mode 100644 index 0000000000..13cc2c2117 --- /dev/null +++ b/avogadro/qtplugins/forcefield/scriptenergy.cpp @@ -0,0 +1,367 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "scriptenergy.h" + +#include +#include + +// formats supported in scripts +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace Avogadro::QtPlugins { + +ScriptEnergy::ScriptEnergy(const QString& scriptFileName_) + : m_interpreter(new QtGui::PythonScript(scriptFileName_)), m_valid(false), + m_inputFormat(NotUsed), m_gradients(false), m_ions(false), + m_radicals(false), m_unitCells(false) +{ + m_elements.reset(); + readMetaData(); +} + +ScriptEnergy::~ScriptEnergy() +{ + delete m_interpreter; +} + +QString ScriptEnergy::scriptFilePath() const +{ + return m_interpreter->scriptFilePath(); +} + +Calc::EnergyCalculator* ScriptEnergy::newInstance() const +{ + return new ScriptEnergy(m_interpreter->scriptFilePath()); +} + +void ScriptEnergy::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 || m_interpreter == nullptr) { + return; // nothing to do + } + + if (!m_unitCells && mol->unitCell()) { + // appendError("Unit cell not supported for this script."); + return; + } + if (!m_ions && mol->totalCharge() != 0) { + // appendError("Ionized molecules not supported for this script."); + return; + } + if (!m_radicals && mol->totalSpinMultiplicity() != 1) { + // appendError("Radical molecules not supported for this script."); + return; + } + + // start the process + // we need a tempory file to write the molecule + QScopedPointer format(createFileFormat(m_inputFormat)); + if (format.isNull()) { + // appendError("Invalid input format."); + return; + } + // get a temporary filename + QString tempPath = QDir::tempPath(); + if (!tempPath.endsWith(QDir::separator())) + tempPath += QDir::separator(); + QString tempPattern = + tempPath + "avogadroenergyXXXXXX." + format->fileExtensions()[0].c_str(); + m_tempFile.setFileTemplate(tempPattern); + if (!m_tempFile.open()) { + // appendError("Error creating temporary file."); + return; + } + + // write the molecule + format->writeFile(m_tempFile.fileName().toStdString(), *mol); + m_tempFile.close(); + + // construct the command line options + QStringList options; + options << "-f" << m_tempFile.fileName(); + + // start the interpreter + m_interpreter->asyncExecute(options); +} + +Real ScriptEnergy::value(const Eigen::VectorXd& x) +{ + if (m_molecule == nullptr || m_interpreter == nullptr) + return 0.0; // nothing to do + + // 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, 1); + + return result.toDouble(); // if conversion fails, returns 0.0 +} + +void ScriptEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) +{ + EnergyCalculator::gradient(x, grad); +} + +ScriptEnergy::Format ScriptEnergy::stringToFormat(const std::string& str) +{ + if (str == "cjson") + return Cjson; + else if (str == "cml") + return Cml; + else if (str == "mdl" || str == "mol" || str == "sdf" || str == "sd") + return Mdl; + else if (str == "pdb") + return Pdb; + else if (str == "xyz") + return Xyz; + return NotUsed; +} + +Io::FileFormat* ScriptEnergy::createFileFormat(ScriptEnergy::Format fmt) +{ + switch (fmt) { + case Cjson: + return new Io::CjsonFormat; + case Cml: + return new Io::CmlFormat; + case Mdl: + return new Io::MdlFormat; + case Pdb: + return new Io::PdbFormat; + case Xyz: + return new Io::XyzFormat; + default: + case NotUsed: + return nullptr; + } +} + +void ScriptEnergy::resetMetaData() +{ + m_valid = false; + m_gradients = false; + m_ions = false; + m_radicals = false; + m_unitCells = false; + m_inputFormat = NotUsed; + m_identifier.clear(); + m_name.clear(); + m_description.clear(); + m_formatString.clear(); + m_elements.reset(); +} + +void ScriptEnergy::readMetaData() +{ + resetMetaData(); + + QByteArray output(m_interpreter->execute(QStringList() << "--metadata")); + + if (m_interpreter->hasErrors()) { + qWarning() << tr("Error retrieving metadata for energy script: %1") + .arg(scriptFilePath()) + << "\n" + << m_interpreter->errorList(); + return; + } + + QJsonParseError parseError; + QJsonDocument doc(QJsonDocument::fromJson(output, &parseError)); + if (parseError.error != QJsonParseError::NoError) { + qWarning() << tr("Error parsing metadata for energy script: %1") + .arg(scriptFilePath()) + << "\n" + << parseError.errorString(); + return; + } + + if (!doc.isObject()) { + qWarning() << tr("Error parsing metadata for energy script: %1\n" + "Result is not a JSON object.\n") + .arg(scriptFilePath()); + return; + } + + const QJsonObject metaData(doc.object()); + + // Read required inputs first. + std::string identifierTmp; + if (!parseString(metaData, "identifier", identifierTmp)) { + qWarning() << "Error parsing metadata for energy script:" + << scriptFilePath() << "\n" + << "Error parsing required member 'identifier'" + << "\n" + << output; + return; + } + m_identifier = identifierTmp; + + std::string nameTmp; + if (!parseString(metaData, "name", nameTmp)) { + qWarning() << "Error parsing metadata for energy script:" + << scriptFilePath() << "\n" + << "Error parsing required member 'name'" + << "\n" + << output; + return; + } + m_name = nameTmp; + + std::string descriptionTmp; + parseString(metaData, "description", descriptionTmp); + m_description = descriptionTmp; // optional + + Format inputFormatTmp = NotUsed; + std::string inputFormatStrTmp; + if (!parseString(metaData, "inputFormat", inputFormatStrTmp)) { + qWarning() << "Error parsing metadata for energy script:" + << scriptFilePath() << "\n" + << "Member 'inputFormat' required for writable formats." + << "\n" + << output; + return; + } + m_formatString = inputFormatStrTmp.c_str(); // for the json key + + // Validate the input format + inputFormatTmp = stringToFormat(inputFormatStrTmp); + if (inputFormatTmp == NotUsed) { + qWarning() << "Error parsing metadata for energy script:" + << scriptFilePath() << "\n" + << "Member 'inputFormat' not recognized:" + << inputFormatStrTmp.c_str() + << "\nValid values are cjson, cml, mdl/sdf, pdb, or xyz.\n" + << output; + return; + } + m_inputFormat = inputFormatTmp; + + // check ions, radicals, unit cells + /* + "unitCell": False, + "gradients": True, + "ion": False, + "radical": False, + */ + if (!metaData["gradients"].isBool()) { + return; // not valid + } + m_gradients = metaData["gradients"].toBool(); + + if (!metaData["unitCell"].isBool()) { + return; // not valid + } + m_unitCells = metaData["unitCell"].toBool(); + + if (!metaData["ion"].isBool()) { + return; // not valid + } + m_ions = metaData["ion"].toBool(); + + if (!metaData["radical"].isBool()) { + return; // not valid + } + m_radicals = metaData["radical"].toBool(); + + // get the element mask + // (if it doesn't exist, the default is no elements anyway) + m_valid = parseElements(metaData); +} + +bool ScriptEnergy::parseString(const QJsonObject& ob, const QString& key, + std::string& str) +{ + if (!ob[key].isString()) + return false; + + str = ob[key].toString().toStdString(); + + return !str.empty(); +} + +void ScriptEnergy::processElementString(const QString& str) +{ + // parse the QString + // first turn any commas into whitespace + QString str2(str); + str2.replace(',', ' '); + // then split on whitespace + QStringList strList = str2.split(QRegExp("\\s+"), QString::SkipEmptyParts); + foreach (QString sstr, strList) { + // these should be numbers or ranges (e.g., 1-84) + if (sstr.contains('-')) { + // range, so split on the dash + QStringList strList2 = sstr.split('-'); + if (strList2.size() != 2) + return; + + // get the two numbers + bool ok; + int start = strList2[0].toInt(&ok); + if (!ok || start < 1 || start > 119) + return; + int end = strList2[1].toInt(&ok); + if (!ok || end < 1 || end > 119) + return; + for (int i = start; i <= end; ++i) + m_elements.set(i); + } + + bool ok; + int i = sstr.toInt(&ok); + if (!ok || i < 1 || i > 119) + return; + + m_elements.set(i); + } +} + +bool ScriptEnergy::parseElements(const QJsonObject& object) +{ + m_elements.reset(); + + // we could either get a string or an array (of numbers) + if (object["elements"].isString()) { + auto str = object["elements"].toString(); + processElementString(str); + + } else if (object["elements"].isArray()) { + QJsonArray arr = object["elements"].toArray(); + for (auto&& i : arr) { + if (i.isString()) { + processElementString(i.toString()); + } else if (i.isDouble()) { + int element = i.toInt(); + if (element >= 1 && element <= 119) // check the range + m_elements.set(element); + } + } + } + return true; +} + +} // namespace Avogadro::QtPlugins diff --git a/avogadro/qtplugins/forcefield/scriptenergy.h b/avogadro/qtplugins/forcefield/scriptenergy.h new file mode 100644 index 0000000000..32783d2ff5 --- /dev/null +++ b/avogadro/qtplugins/forcefield/scriptenergy.h @@ -0,0 +1,107 @@ +/****************************************************************************** + 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_SCRIPTENERGY_H +#define AVOGADRO_QTPLUGINS_SCRIPTENERGY_H + +#include + +#include + +#include +#include +#include + +class QJsonObject; + +namespace Avogadro { + +namespace Io { +class FileFormat; +} + +namespace QtGui { +class PythonScript; +} + +namespace QtPlugins { + +class ScriptEnergy : public Avogadro::Calc::EnergyCalculator +{ + Q_DECLARE_TR_FUNCTIONS(ScriptEnergy) + +public: + /** Formats that may be written to the script's input. */ + enum Format + { + NotUsed, + Cjson, + Cml, + Mdl, // sdf + Pdb, + Xyz + }; + + ScriptEnergy(const QString& scriptFileName = ""); + ~ScriptEnergy() override; + + QString scriptFilePath() const; + + Format inputFormat() const { return m_inputFormat; } + bool isValid() const { return m_valid; } + + 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; } + + Core::Molecule::ElementMask elements() const override { return m_elements; } + bool supportsGradients() const { return m_gradients; } + bool supportsIons() const { return m_ions; } + bool supportsRadicals() const { return m_radicals; } + bool supportsUnitCells() const { return m_unitCells; } + + // 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; + +private: + static Format stringToFormat(const std::string& str); + static Io::FileFormat* createFileFormat(Format fmt); + void resetMetaData(); + void readMetaData(); + bool parseString(const QJsonObject& ob, const QString& key, std::string& str); + void processElementString(const QString& str); + bool parseElements(const QJsonObject& ob); + +private: + QtGui::PythonScript* m_interpreter; + Format m_inputFormat; + Core::Molecule* m_molecule; + + // what's supported by this script + Core::Molecule::ElementMask m_elements; + bool m_valid; + bool m_gradients; + bool m_ions; + bool m_radicals; + bool m_unitCells; + + std::string m_identifier; + std::string m_name; + std::string m_description; + QString m_formatString; + QTemporaryFile m_tempFile; +}; + +} // namespace QtPlugins +} // namespace Avogadro + +#endif // AVOGADRO_QTPLUGINS_SCRIPTENERGY_H diff --git a/avogadro/qtplugins/forcefield/scripts/gfn1.py b/avogadro/qtplugins/forcefield/scripts/gfn1.py new file mode 100644 index 0000000000..3673b1a8ad --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/gfn1.py @@ -0,0 +1,102 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + import numpy as np + from xtb.libxtb import VERBOSITY_MUTED + from xtb.interface import Calculator, Param + imported = True +except ImportError: + imported = False + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "GFN1", + "identifier": "GFN1", + "description": "Calculate GFN1-xtb energies and gradients", + "inputFormat": "cjson", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": True, + "radical": True, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + with open(filename, "r") as f: + mol_cjson = json.load(f) + + # first setup the calculator + atoms = np.array(mol_cjson["atoms"]["elements"]["number"]) + coord_list = mol_cjson["atoms"]["coords"]["3d"] + coordinates = np.array(coord_list, dtype=float).reshape(-1, 3) + # .. we need to convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # check for total charge + # and spin multiplicity + charge = None # neutral + spin = None # singlet + if "properties" in mol_cjson: + if "totalCharge" in mol_cjson["properties"]: + charge = mol_cjson["properties"]["totalCharge"] + if "totalSpinMultiplicity" in mol_cjson["properties"]: + spin = mol_cjson["properties"]["totalSpinMultiplicity"] + + calc = Calculator(Param.GFN1xTB, atoms, coordinates, + charge=charge, uhf=spin) + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() + + # we loop forever - Avogadro will kill the process when done + while(True): + # first print the energy of these coordinates + print(res.get_energy()) # in Hartree + + # 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=" ") + # .. convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # update the calculator and run a new calculation + calc.update(coordinates) + calc.singlepoint(res) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("GFN1 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("xtb-python is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/gfn2.py b/avogadro/qtplugins/forcefield/scripts/gfn2.py new file mode 100644 index 0000000000..e7035d147b --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/gfn2.py @@ -0,0 +1,102 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + import numpy as np + from xtb.libxtb import VERBOSITY_MUTED + from xtb.interface import Calculator, Param + imported = True +except ImportError: + imported = False + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "GFN2", + "identifier": "GFN2", + "description": "Calculate GFN2-xtb energies and gradients", + "inputFormat": "cjson", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": True, + "radical": True, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + with open(filename, "r") as f: + mol_cjson = json.load(f) + + # first setup the calculator + atoms = np.array(mol_cjson["atoms"]["elements"]["number"]) + coord_list = mol_cjson["atoms"]["coords"]["3d"] + coordinates = np.array(coord_list, dtype=float).reshape(-1, 3) + # .. we need to convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # check for total charge + # and spin multiplicity + charge = None # neutral + spin = None # singlet + if "properties" in mol_cjson: + if "totalCharge" in mol_cjson["properties"]: + charge = mol_cjson["properties"]["totalCharge"] + if "totalSpinMultiplicity" in mol_cjson["properties"]: + spin = mol_cjson["properties"]["totalSpinMultiplicity"] + + calc = Calculator(Param.GFN2xTB, atoms, coordinates, + charge=charge, uhf=spin) + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() + + # we loop forever - Avogadro will kill the process when done + while(True): + # first print the energy of these coordinates + print(res.get_energy()) # in Hartree + + # 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=" ") + # .. convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # update the calculator and run a new calculation + calc.update(coordinates) + calc.singlepoint(res) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("GFN2 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("xtb-python is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/gfnff.py b/avogadro/qtplugins/forcefield/scripts/gfnff.py new file mode 100644 index 0000000000..c6e6223f9c --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/gfnff.py @@ -0,0 +1,104 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + import numpy as np + from xtb.libxtb import VERBOSITY_MUTED + from xtb.interface import Calculator, Param + + imported = True +except ImportError: + imported = False + + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "GFN-FF", + "identifier": "GFN-FF", + "description": "Calculate GFNFF-xtb energies and gradients", + "inputFormat": "cjson", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": True, + "radical": False, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + with open(filename, "r") as f: + mol_cjson = json.load(f) + + # first setup the calculator + atoms = np.array(mol_cjson["atoms"]["elements"]["number"]) + coord_list = mol_cjson["atoms"]["coords"]["3d"] + coordinates = np.array(coord_list, dtype=float).reshape(-1, 3) + # .. we need to convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # check for total charge + # and spin multiplicity + charge = None # neutral + spin = None # singlet + if "properties" in mol_cjson: + if "totalCharge" in mol_cjson["properties"]: + charge = mol_cjson["properties"]["totalCharge"] + if "totalSpinMultiplicity" in mol_cjson["properties"]: + spin = mol_cjson["properties"]["totalSpinMultiplicity"] + + calc = Calculator(Param.GFNFF, atoms, coordinates, + charge=charge, uhf=spin) + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() + + # 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 + + # 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=" ") + # .. convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # update the calculator and run a new calculation + calc.update(coordinates) + calc.singlepoint(res) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("GFN2 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("xtb-python is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/mmff94.py b/avogadro/qtplugins/forcefield/scripts/mmff94.py new file mode 100644 index 0000000000..e0f826e3da --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/mmff94.py @@ -0,0 +1,85 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + from openbabel import pybel + import numpy as np + + imported = True +except ImportError: + imported = False + + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "MMFF94", + "identifier": "MMFF94", + "description": "Calculate MMFF94 energies and gradients", + "inputFormat": "cml", + "elements": "1,6-9,14-17,35,53", + "unitCell": False, + "gradients": True, + "ion": False, + "radical": False, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + mol = next(pybel.readfile("cml", filename)) + + ff = pybel._forcefields["mmff94"] + success = ff.Setup(mol.OBMol) + if not success: + # should never happen, but just in case + sys.exit("MMFF94 force field setup failed") + + # 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]) + + # update the molecule geometry for the next energy + ff.SetCoordinates(mol.OBMol) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("MMFF94 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("pybel is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/uff.py b/avogadro/qtplugins/forcefield/scripts/uff.py new file mode 100644 index 0000000000..7d35226262 --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/uff.py @@ -0,0 +1,85 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + from openbabel import pybel + import numpy as np + + imported = True +except ImportError: + imported = False + + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "UFF", + "identifier": "UFF", + "description": "Calculate UFF energies and gradients", + "inputFormat": "cml", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": False, + "radical": False, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + mol = next(pybel.readfile("cml", filename)) + + ff = pybel._forcefields["uff"] + success = ff.Setup(mol.OBMol) + if not success: + # should never happen, but just in case + sys.exit("UFF force field setup failed") + + # 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]) + + # update the molecule geometry for the next energy + ff.SetCoordinates(mol.OBMol) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("UFF calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("pybel is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp b/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp index 2c789616d0..5fa58a1a6c 100644 --- a/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp +++ b/avogadro/qtplugins/scriptcharges/scriptchargemodel.cpp @@ -277,7 +277,7 @@ void ScriptChargeModel::readMetaData() if (!parseString(metaData, "identifier", identifierTmp)) { qWarning() << "Error parsing metadata for charge script:" << scriptFilePath() << "\n" - << "Error parsing required member 'operations'" + << "Error parsing required member 'identifier'" << "\n" << output; return;