Skip to content

Commit

Permalink
More work
Browse files Browse the repository at this point in the history
Signed-off-by: Geoff Hutchison <[email protected]>
  • Loading branch information
ghutchis committed Oct 14, 2023
1 parent dfad831 commit 44a76c7
Show file tree
Hide file tree
Showing 8 changed files with 539 additions and 51 deletions.
6 changes: 5 additions & 1 deletion avogadro/calc/energymanager.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
85 changes: 42 additions & 43 deletions avogadro/qtplugins/forcefield/forcefield.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,10 @@
#include <QtWidgets/QAction>
#include <QtWidgets/QMessageBox>

#include <QProgressDialog>
#include <QWriteLocker>
#include <QMutex>
#include <QMutexLocker>
#include <QProgressDialog>
#include <QWriteLocker>

#include <avogadro/qtgui/avogadropython.h>
#include <avogadro/qtgui/filebrowsewidget.h>
Expand All @@ -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();

Expand All @@ -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);
Expand Down Expand Up @@ -106,7 +100,7 @@ void Forcefield::setMolecule(QtGui::Molecule* mol)

m_molecule = mol;

// @todo set molecule for a calculator
// @todo set the method
}

void Forcefield::refreshScripts()
Expand All @@ -123,9 +117,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<EnergyCalculator> solver;
cppoptlib::ConjugatedGradientDescentSolver<EnergyCalculator> solver;
cppoptlib::LbfgsSolver<EnergyCalculator> solver;

int n = m_molecule->atomCount();
Core::Array<Vector3> pos = m_molecule->atomPositions3d();
Expand All @@ -138,42 +130,49 @@ void Forcefield::optimize()

// Create a Criteria class to adjust stopping criteria
cppoptlib::Criteria<Real> crit = cppoptlib::Criteria<Real>::defaults();
// @todo allow criteria to be set
crit.iterations = 5;
crit.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);
Expand Down Expand Up @@ -212,5 +211,5 @@ void Forcefield::freezeSelected()
}
}

} // end QtPlugins
}
} // namespace QtPlugins
} // namespace Avogadro
21 changes: 14 additions & 7 deletions avogadro/qtplugins/forcefield/forcefield.h
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,11 @@ class QAction;
class QDialog;

namespace Avogadro {

namespace Calc {
class EnergyCalculator;
}

namespace QtPlugins {

/**
Expand Down Expand Up @@ -65,18 +70,20 @@ private slots:

private:
QList<QAction*> 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<QString, QString> m_forcefieldScripts;

const Io::FileFormat* m_outputFormat;
QString m_tempFileName;
};

}
}

Expand Down
102 changes: 102 additions & 0 deletions avogadro/qtplugins/forcefield/scripts/gfn1.py
Original file line number Diff line number Diff line change
@@ -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])
Loading

0 comments on commit 44a76c7

Please sign in to comment.