OEFF Examples¶
Simple Functions and Optimization¶
Solving a simple equation¶
The following example illustrates how to define a simple objective function. By deriving the objective function from OEFunc2, we can find the roots of the simple quadratic equation using OENewtonOpt optimizer. A class derived from OEFunc2 must contain analytical gradients and Hessians. This example expects a number as input for the initial guess to solve the simple function.
Listing 1: Define a simple objective function
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include <string>
//////////////////////////////////////////////////////////////////////////////
// The following function is a contrived example to show how to //
// write a user-defined objective function and optimize it. The simple //
// function defines a quadradic equation, and contains expressions for //
// gradient and second derivative of the function. //
//////////////////////////////////////////////////////////////////////////////
class Function : public OEOpt::OEFunc2
{
public:
unsigned int NumVar() const {return 1;}
double operator()(const double* x) {return (x[0]*x[0]-7*x[0]+63);}
double operator()(const double* x, double* g)
{
g[0] = 2.0*x[0]-7.0;
return operator()(x);
}
bool operator()(const double* x, double* h, double* g)
{
g[0] = 2.0*x[0]-7.0;
h[0] = 2.0;
return true;
}
OESystem::OEIterBase<OEOpt::OEFComponent> *GetFComponents(const double*)
{
//-Cannot provide a meaningful implementation
return 0;
}
};
int main(int argc, char *argv[])
{
if (argc!=2)
OESystem::OEThrow.Usage("%s <initial guess>", argv[0]);
double x[1];
try
{
x[0] = atof(argv[1]);
}
catch(...)
{
OESystem::OEThrow.Usage("%s <initial guess (expecting a number)>", argv[0]);
}
Function func;
// Calculate function value at given x
double value = func(x);
OESystem::OEThrow.Info("Function value at x = %f is %f", x[0], value);
// Optimize function using Newton Optimizer
double xOpt[1];
OEOpt::OENewtonOpt optimizer;
value = optimizer(func, x, xOpt);
OESystem::OEThrow.Info("Function has a minimia at x = %f and the minimum value is %f", xOpt[0], value);
return 0;
}
Download code
Using checkpoints for optimization monitoring¶
The following example illustrates how to define checkpoints and use then along with optimizers to monitor progress during optimization. For simplicity, a simple quadratic equation is defined as objective function and derived from OEFunc1. The quadratic equation is solved using the OEBFGSOpt optimizer. A class derived from OEFunc1 must contain analytical gradients. This examples expects a number as input for the initial guess to solve the simple function.
Listing 2: Define and use checkpoints
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include <string>
//////////////////////////////////////////////////////////////////////////////
// The following function displays how to use a checkpoint during function //
// optimization. A simple quadradic equation is minimized for which the //
// analytical gradient is provided but the second derivative is not. //
//////////////////////////////////////////////////////////////////////////////
class Function : public OEOpt::OEFunc1
{
public:
unsigned int NumVar() const {return 1;}
double operator()(const double* x) {return (x[0]*x[0]-7*x[0]+63);}
double operator()(const double* x, double* g)
{
g[0] = 2.0*x[0]-7.0;
return operator()(x);
}
OESystem::OEIterBase<OEOpt::OEFComponent> *GetFComponents(const double*)
{
//-Cannot provide a meaningful implementation
return 0;
}
};
class ChkPt: public OEOpt::OECheckpoint1
{
public:
bool operator()(unsigned int iteration, unsigned int nvar,
double fval, const double* var,
unsigned int state)
{
// Intermediate information during optimization
OESystem::OEThrow.Info("iteration: %d nvar: %d x: %f value: %f state: %d", iteration, nvar, var[0], fval, state);
// To demonstrate how to force quitting optimization from checkpoint
// this returns false when iteration is 5
return (iteration >= 5? false:true);
}
bool operator()(unsigned int iteration, unsigned int nvar,
double fval, double gnorm, const double* var,const double* grad, unsigned int state)
{
// Intermediate information during optimization
OESystem::OEThrow.Info("iteration: %d nvar: %d x: %f value: %f gnorm: %f gradient: %f state: %d",
iteration, nvar, var[0], fval, gnorm, grad[0], state);
// To demonstrate how to force quitting optimization from checkpoint
// this returns false when iteration is 5
return (iteration >= 5? false:true);
}
};
int main(int argc, char *argv[])
{
if (argc!=2)
OESystem::OEThrow.Usage("%s <initial guess>", argv[0]);
double x[1];
try
{
x[0] = atof(argv[1]);
}
catch(...)
{
OESystem::OEThrow.Usage("%s <initial guess(expecting a number)>", argv[0]);
}
Function func;
// Calculate function value at given x
double value = func(x);
OESystem::OEThrow.Info("Function value at x = %f is %f", x[0], value);
// Optimize function using BFGS Optimizer and checkpoint
double xOpt[1];
OEOpt::OEBFGSOpt optimizer;
ChkPt checkpoint;
value = optimizer(func, &checkpoint, x, xOpt);
OESystem::OEThrow.Info("Function has a minimia at x = %f and the minimum value is %f", xOpt[0], value);
return 0;
}
Download code
Molecule Functions¶
User defined molecule function¶
The following example illustrates how to define an objective function within the context of a molecule. Generally speaking, a molecule function (OEMolFunc) defines some sort of interaction involving a part or all of a molecule. For simplicity, a simple energy function is defined that consists sum of square of distance of all the atoms in the molecule. The molecule function is optimized using the OEBFGSOpt optimizer. A class derived from OEMolFunc1 must contain analytical gradients.
Listing 3: Define an objective function within the context of a molecule
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
//////////////////////////////////////////////////////////////////////////////
// The following function is a contrived example to show how to //
// write a user-defined function. The energy function is the square //
// of the distance from the origin. The derivative then is two times the //
// coordinate component. The function will attempt to place all atoms of //
// a molecule at the origin. //
//////////////////////////////////////////////////////////////////////////////
class Harmonic : public OEMolPotential::OEMolFunc1
{
private:
unsigned int natoms;
public:
unsigned int NumVar() const {return 3*natoms;}
bool Setup(const OEChem::OEMolBase& mol)
{
natoms = mol.GetMaxAtomIdx();
return true;
}
double operator()(const double *coord)
{
double energy = 0.0;
for (unsigned int i = 0;i < natoms;++i)
{
energy += coord[3*i ] * coord[3*i ]; //x distance from zero
energy += coord[3*i+1] * coord[3*i+1]; //y distance from zero
energy += coord[3*i+2] * coord[3*i+2]; //z distance from zero
}
return energy;
}
double operator()(const double *coord, double* g)
{
double energy = 0.0;
for (unsigned int i = 0;i < natoms;++i)
{
g[3*i ] += 2.0 * coord[3*i ]; //x gradient
g[3*i+1] += 2.0 * coord[3*i+1]; //y gradient
g[3*i+2] += 2.0 * coord[3*i+2]; //z gradient
energy += (coord[3*i ]*coord[3*i ]) +
(coord[3*i+1]*coord[3*i+1]) +
(coord[3*i+2]*coord[3*i+2]);
}
return energy;
}
OESystem::OEIterBase<OEOpt::OEFComponent> *GetFComponents(const double*)
{
//-Cannot provide a meaningful implementation
return 0;
}
Harmonic* CreateCopy() const
{
return new Harmonic(*this);
}
};
int main(int argc, char *argv[])
{
if (argc!=2)
OESystem::OEThrow.Usage("%s <input>", argv[0]);
OEChem::oemolistream ifs;
if (!ifs.open(argv[1]))
OESystem::OEThrow.Fatal("Unable to open %s for reading", argv[1]);
OEChem::OEGraphMol mol;
OEChem::OEReadMolecule(ifs, mol);
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
mol.GetCoords(&vecCoords[0]);
Harmonic hermonic;
hermonic.Setup(mol);
OEOpt::OEBFGSOpt optimizer;
double energy = optimizer(hermonic, &vecCoords[0], &vecCoords[0]);
OESystem::OEThrow.Info("Optimized energy: %f", energy);
return 0;
}
Download code
Small Molecule Force fields¶
Single ligand in vacuum¶
The following example illustrates how to calculate energy and optimize a single
ligand in vacuum. The molecule needs to be prepared (OEForceField::PrepMol),
followed by a call to OEMolFunc::Setup to create the interactions. Optimization is carried out using the OEBFGSOpt optimizer.
Listing 4: Calculate energy and optimize a single ligand in vacuum
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
/////////////////////////////////////////////////////////////////////////////
// The following is an example to show how to evaluate energies //
// of a small molecule and obtain various energy components. //
/////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "parsley");
ffType.AddLegalValue("parsley");
ffType.AddLegalValue("mmff94");
ffType.AddLegalValue("mmff94s");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEMolPotential::OEForceField* GetFF() const
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff94")
return new OEMolMMFF::OEMMFF();
else if (ff == "mmff94s")
return new OEMolMMFF::OEMMFF(OEMolMMFF::OEMMFF94sParams());
else
return new OEMolSmirnoff::OEParsley();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OESimpleAppOptions opts(ffOpts, "FFOptimize", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D);
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
std::unique_ptr<OEMolPotential::OEForceField> ff(ffOpts.GetFF());
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
if (!ff->PrepMol(mol) || !ff->Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
// Calculate energy
double energy = (*ff)(&vecCoords[0]);
OESystem::OEThrow.Info("Initial energy: %f kcal/mol", energy);
// Optimize the ligand
OEOpt::OEBFGSOpt optimizer;
energy = optimizer(*ff, &vecCoords[0], &vecCoords[0]);
OESystem::OEThrow.Info("Optimized energy: %f kcal/mol", energy);
conf->SetCoords(&vecCoords[0]);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
Energy components of a single ligand in vacuum¶
The following example illustrates how to obtain various intermolecular energy components of a single ligand in vacuum.
Listing 5: Obtain intermolecular energy components of a single ligand in vacuum
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
/////////////////////////////////////////////////////////////////////////////
// The following is an example to show how to evaluate energies //
// of a small molecule and obtain various energy components. //
/////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "parsley");
ffType.AddLegalValue("parsley");
ffType.AddLegalValue("mmff94");
ffType.AddLegalValue("mmff94s");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEMolPotential::OEForceField* GetFF() const
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff94")
return new OEMolMMFF::OEMMFF();
else if (ff == "mmff94s")
return new OEMolMMFF::OEMMFF(OEMolMMFF::OEMMFF94sParams());
else
return new OEMolSmirnoff::OEParsley();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OESimpleAppOptions opts(ffOpts, "FFEnergy", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D);
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
std::unique_ptr<OEMolPotential::OEForceField> ff(ffOpts.GetFF());
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
if (!ff->PrepMol(mol) || !ff->Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
OESystem::OEIter<OEOpt::OEFComponent> fcomp = ff->GetFComponents(&vecCoords[0]);
for (; fcomp; ++fcomp)
OESystem::OEThrow.Info("Component: %s Energy: %f kcal/mol", fcomp->name.c_str(), fcomp->value);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
MMFF interactions of a single ligand¶
The following example illustrates how to obtain various interactions of a single ligand in the context of a force field.
Listing 6: Obtain interactions of a single ligand in the context of a force field
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
/////////////////////////////////////////////////////////////////////////////
// The following example demonstrates how to obtain interactions //
// of a small molecule in the context of a force field. //
/////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "parsley");
ffType.AddLegalValue("parsley");
ffType.AddLegalValue("mmff94");
ffType.AddLegalValue("mmff94s");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEMolPotential::OEForceField* GetFF() const
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff94")
return new OEMolMMFF::OEMMFF();
else if (ff == "mmff94s")
return new OEMolMMFF::OEMMFF(OEMolMMFF::OEMMFF94sParams());
else
return new OEMolSmirnoff::OEParsley();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OESimpleAppOptions opts(ffOpts, "FFInteractions", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D);
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
std::unique_ptr<OEMolPotential::OEForceField> ff(ffOpts.GetFF());
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
if (!ff->PrepMol(mol) || !ff->Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
OESystem::OEIter<OEMolPotential::OEInteraction> intc = ff->GetInteractions(mol, &vecCoords[0]);
for(; intc; ++intc)
{
std::vector<double> vecGrads(intc->NumValues());
OESystem::OEThrow.Info("Interaction: %s Value: %f", intc->GetName().c_str(), intc->GetValues(&vecGrads[0]));
for (OESystem::OEIter<OEChem::OEAtomBase>atom = intc->GetAtoms(); atom; ++atom)
OESystem::OEThrow.Info("Atom index: %d", atom->GetIdx());
}
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
Protein-Ligand Complexes¶
Protein-ligand optimization¶
The following example illustrates how to set up a ligand optimization within the context of a protein contained in a design unit.
Listing 7: Set up a ligand optimization
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
//////////////////////////////////////////////////////////////////////////////
// The following example demonstrates how to perform ligand optimization //
// in the context of a protein using a complex force field. //
//////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "ff14sb_parsley");
ffType.AddLegalValue("ff14sb_parsley");
ffType.AddLegalValue("mmff_amber");
ffType.AddLegalValue("mmff");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEFF::OEComplexFF* GetFF() const
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff")
return new OEFF::OEMMFFComplex();
else if (ff == "mmff_amber")
return new OEFF::OEMMFFAmberComplex();
else
return new OEFF::OEFF14SBParsleyComplex();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OERefInputAppOptions opts(ffOpts, "FFComplexOpt", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::DU, "-protein");
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEPlatform::oeifstream rfs;
if (!rfs.open(opts.GetRefFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetRefFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
OEBio::OEDesignUnit du;
if(!OEBio::OEReadDesignUnit(rfs, du))
OESystem::OEThrow.Fatal("Failed to read design unit");
OEChem::OEGraphMol protein;
du.GetProtein(protein);
std::unique_ptr<OEFF::OEComplexFF> ff(ffOpts.GetFF());
if(!ff->SetupHost(protein))
OESystem::OEThrow.Fatal("Failed to setup protein: %s", protein.GetTitle());
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
OESystem::OEThrow.Info("Optimizing %s", mol.GetTitle());
if (!ff->Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
// Calculate energy
double energy = (*ff)(&vecCoords[0]);
OESystem::OEThrow.Info("Initial energy: %f kcal/mol", energy);
// Optimize the ligand
OEOpt::OEBFGSOpt optimizer;
energy = optimizer(*ff, &vecCoords[0], &vecCoords[0]);
OESystem::OEThrow.Info("Optimized energy: %f kcal/mol", energy);
conf->SetCoords(&vecCoords[0]);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
Using adaptors¶
Optimizing a single ligand with fixed atoms¶
The following example illustrates how to fix a subset of atoms using the OESubsetAdaptor during a single ligand optimization in vacuum. The adaptor is initialized with the force field interactions of the ligand, and passed as the objective function to be optimized. Methods of the adaptor are used to convert between the Cartesian coordinates of the ligand and the adaptor variables.
Listing 8: Fix a subset of atoms
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
/////////////////////////////////////////////////////////////////////////////
// The following is an example to show how to evaluate energy //
// and optimize a small molecule, with keeping a subset of atoms fixed. //
/////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "parsley");
ffType.AddLegalValue("parsley");
ffType.AddLegalValue("mmff94");
ffType.AddLegalValue("mmff94s");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEMolPotential::OEForceField* GetFF() const
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff94")
return new OEMolMMFF::OEMMFF();
else if (ff == "mmff94s")
return new OEMolMMFF::OEMMFF(OEMolMMFF::OEMMFF94sParams());
else
return new OEMolSmirnoff::OEParsley();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OESimpleAppOptions opts(ffOpts, "FFSubset", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D);
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
std::unique_ptr<OEMolPotential::OEForceField> ff(ffOpts.GetFF());
// Setup adaptor. The first (false) means not to pass ownership of ff,
// and the second (false) means not to exclude interactions related
// to the subset which would be fixed for calculations.
OEMolPotential::OESubsetAdaptor adaptor(*ff, false, false);
// Use a simple atoms predicate for the subset
adaptor.Set(OEChem::OEIsCarbon());
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
if (!ff->PrepMol(mol) || !adaptor.Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
// Get adaptor variables set corresponding to the coordinates
std::vector<double> vecX(adaptor.NumVar());
adaptor.GetVar(&vecX[0], &vecCoords[0]);
// Calculate energy
double energy = adaptor(&vecX[0]);
OESystem::OEThrow.Info("Initial energy: %f kcal/mol", energy);
// Optimize the ligand
OEOpt::OEBFGSOpt optimizer;
energy = optimizer(adaptor, &vecX[0], &vecX[0]);
OESystem::OEThrow.Info("Optimized energy: %f kcal/mol", energy);
// Get optimized coordinates corresponding to the adaptor optimized variables
adaptor.AdaptVar(&vecCoords[0], &vecX[0]);
conf->SetCoords(&vecCoords[0]);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
Optimizing a single ligand torsions¶
The following example illustrates how to optimize a set of torsion angles using the OETorAdaptor for a single ligand in vacuum. The adaptor is initialized with the force field interactions of the ligand, and passed as the objective function to be optimized. Methods of the adaptor are used to convert between the Cartesian coordinates of the ligand and the adaptor variables.
Listing 9: Optimize a set of torsion angles
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
/////////////////////////////////////////////////////////////////////////////
// The following is an example to show how to evaluate energy //
// and optimize a set of torsions in a small molecule //
/////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "parsley");
ffType.AddLegalValue("parsley");
ffType.AddLegalValue("mmff94");
ffType.AddLegalValue("mmff94s");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEMolPotential::OEForceField* GetFF() const
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff94")
return new OEMolMMFF::OEMMFF();
else if (ff == "mmff94s")
return new OEMolMMFF::OEMMFF(OEMolMMFF::OEMMFF94sParams());
else
return new OEMolSmirnoff::OEParsley();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OESimpleAppOptions opts(ffOpts, "FFTorAdaptor", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D);
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
std::unique_ptr<OEMolPotential::OEForceField> ff(ffOpts.GetFF());
// Torsion adaptor. The first (false) means not to pass ownership of ff,
// and the second (false) means not to exclude any interactions.
OEMolPotential::OETorAdaptor adaptor(*ff, false, false);
// Use a simple rotor predicate for the torsions
adaptor.Set(OEChem::OEIsRotor());
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
if (!ff->PrepMol(mol) || !adaptor.Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
// Get adaptor variables set corresponding to the coordinates
std::vector<double> vecX(adaptor.NumVar());
adaptor.GetVar(&vecX[0], &vecCoords[0]);
// Calculate energy
double energy = adaptor(&vecX[0]);
OESystem::OEThrow.Info("Initial energy: %f kcal/mol", energy);
// Optimize the ligand
OEOpt::OEBFGSOpt optimizer;
energy = optimizer(adaptor, &vecX[0], &vecX[0]);
OESystem::OEThrow.Info("Optimized energy: %f kcal/mol", energy);
// Get optimized coordinates corresponding to the adaptor optimized variables
adaptor.AdaptVar(&vecCoords[0], &vecX[0]);
conf->SetCoords(&vecCoords[0]);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
Optimizing rigid ligand in protein¶
The following example illustrates how to perform rigid optimization of a ligand in the context of a protein using the OEQuatAdaptor. The adaptor is initialized with the protein-ligand interactions, and passed as the objective function to be optimized. Methods of the adaptor are used to convert between the Cartesian coordinates of the ligand and the adaptor variables.
Listing 10: Perform rigid optimization of a ligand
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
//////////////////////////////////////////////////////////////////////////////
// The following example demonstrates how to perform rigid optimization of //
// a ligand in the context of a protein. The OEQuatAdaptor it used for //
// rigid optimization. //
//////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "ff14sb_parsley");
ffType.AddLegalValue("ff14sb_parsley");
ffType.AddLegalValue("mmff_amber");
ffType.AddLegalValue("mmff");
ffType.SetBrief("Force field type");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
OEFF::OEComplexFF* GetFF()
{
std::string ff = m_fftype->GetStringValue();
if (ff.empty())
ff = m_fftype->GetStringDefault();
if (ff == "mmff")
return new OEFF::OEMMFFComplex();
else if (ff == "mmff_amber")
return new OEFF::OEMMFFAmberComplex();
else
return new OEFF::OEFF14SBParsleyComplex();
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OERefInputAppOptions opts(ffOpts, "FFComplexQuat", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D, OEChem::OEFileStringType::DU, "-protein");
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEPlatform::oeifstream rfs;
if (!rfs.open(opts.GetRefFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetRefFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
OEBio::OEDesignUnit du;
if(!OEBio::OEReadDesignUnit(rfs, du))
OESystem::OEThrow.Fatal("Failed to read design unit");
OEChem::OEGraphMol protein;
du.GetProtein(protein);
OEFF::OEComplexFF* ff = ffOpts.GetFF();
if(!ff->SetupHost(protein))
OESystem::OEThrow.Fatal("Failed to setup protein: %s", protein.GetTitle());
OEMolPotential::OEQuatAdaptor adap(*ff);
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
OESystem::OEThrow.Info("Optimizing %s", mol.GetTitle());
if (!adap.Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
std::vector<double> vecX(adap.NumVar());
adap.GetVar(&vecX[0], &vecCoords[0]);
// Calculate energy
double energy = adap(&vecX[0]);
OESystem::OEThrow.Info("Initial energy: %f kcal/mol", energy);
// Optimize the ligand
OEOpt::OEBFGSOpt optimizer;
energy = optimizer(adap, &vecX[0], &vecX[0]);
OESystem::OEThrow.Info("Optimized energy: %f kcal/mol", energy);
adap.AdaptVar(&vecCoords[0], &vecX[0]);
conf->SetCoords(&vecCoords[0]);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code
Custom Force Fields¶
Custom Smirnoff¶
The following example illustrates how to load custom SMIRNOFF force field parameters from an OFFXML file and use that to calculate energy and optimize a single ligand in vacuum.
Listing 11: Load custom SMIRNOFF force field parameters
/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/
#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oebio.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oesmirnoff.h"
#include "oeff.h"
/////////////////////////////////////////////////////////////////////////////
// The following example demonstrates how to load a custom Smirnoff //
// forcefield, and use that to optiiza a ligand. //
/////////////////////////////////////////////////////////////////////////////
class FFOptions : public OESystem::OEOptions
{
public:
FFOptions(std::string name = "FFOptions")
: OESystem::OEOptions(name)
{
OESystem::OEStringParameter ffType("-ff", "");
ffType.SetBrief("Force field XML file");
m_fftype = AddParameter(ffType);
}
FFOptions* CreateCopy() const { return new FFOptions(*this); }
std::string GetXML() const
{
std::string ffxml = m_fftype->GetStringValue();
if(ffxml.empty())
ffxml = m_fftype->GetStringDefault();
return ffxml;
}
private:
OESystem::OEParameter* m_fftype;
};
int main(int argc, char* argv[])
{
FFOptions ffOpts;
OEChem::OESimpleAppOptions opts(ffOpts, "FFCustomSmirnoff", OEChem::OEFileStringType::Mol3D,
OEChem::OEFileStringType::Mol3D);
if (OESystem::OEConfigureOpts(opts, argc, argv, false) == OESystem::OEOptsConfigureStatus::Help)
return 0;
ffOpts.UpdateValues(opts);
OEChem::oemolistream ifs;
if (!ifs.open(opts.GetInFile()))
OESystem::OEThrow.Fatal("Unable to open %s for reading", opts.GetInFile().c_str());
OEChem::oemolostream ofs;
if (!ofs.open(opts.GetOutFile()))
OESystem::OEThrow.Fatal("Unable to open %s for writing", opts.GetOutFile().c_str());
OEMolSmirnoff::OESmirnoffParams params;
if (!ffOpts.GetXML().empty())
if (!params.Load(ffOpts.GetXML()))
OESystem::OEThrow.Fatal("Unable to load force field parameters: %s", ffOpts.GetXML().c_str());
OEMolSmirnoff::OESmirnoff ff(params);
OEChem::OEMol mol;
while (OEChem::OEReadMolecule(ifs, mol))
{
if (!ff.PrepMol(mol) || !ff.Setup(mol))
{
OESystem::OEThrow.Warning("Unable to process molecule: title = %s", mol.GetTitle());
OEChem::OEWriteMolecule(ofs, mol);
continue;
}
std::vector<double> vecCoords(3*mol.GetMaxAtomIdx());
for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
{
OESystem::OEThrow.Info("Molecule: %s Conformer: %d", mol.GetTitle(), conf->GetIdx()+1);
conf->GetCoords(&vecCoords[0]);
// Calculate energy
double energy = ff(&vecCoords[0]);
OESystem::OEThrow.Info("Initial energy: %f kcal/mol", energy);
// Optimize the ligand
OEOpt::OEBFGSOpt optimizer;
energy = optimizer(ff, &vecCoords[0], &vecCoords[0]);
OESystem::OEThrow.Info("Optimized energy: %f kcal/mol", energy);
conf->SetCoords(&vecCoords[0]);
}
OEChem::OEWriteMolecule(ofs, mol);
}
return 0;
}
Download code