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 examples expects a number as input for the initial guess to solve the simple function.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#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 = %d is %d", 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 = %d and the minimum value is %d", xOpt[0], value);

  return 0;
}

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.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#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: %n nvar: %n x: %d value: %d state: %n", 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: %n nvar: %n x: %d value: %d gnorm: %d gradient: %d state: %n", 
        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 = %d is %d", 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 = %d and the minimum value is %d", xOpt[0], value);

  return 0;
}

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.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#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;
  }
};


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);
  OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
  mol.GetCoords(vecCoords);

  Harmonic hermonic;
  hermonic.Setup(mol);
  OEOpt::OEBFGSOpt optimizer;
  double energy = optimizer(hermonic, vecCoords, vecCoords);
  OESystem::OEThrow.Info("Optimized energy: %d", energy);

  return 0;
}

Single ligand and MMFF

Single ligand in vacuum

The following example illustrates how to calculate energy and optimize a single ligand in vacuum. The OEMMFF force field is used as is for this example. 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.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following function is an example to show how to evaluate energy     //
//  and optimize a small molecule using the MMFF force field.               //
//////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[])
{
  if (argc!=3)
    OESystem::OEThrow.Usage("%s <input> <output>", argv[0]);
    
  OEChem::oemolistream ifs;
  if (!ifs.open(argv[1]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", argv[1]);

  OEChem::oemolostream ofs;
  if (!ofs.open(argv[2]))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", argv[2]);

  OEChem::OEMol mol;
  while (OEChem::OEReadMolecule(ifs, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    OEMolMMFF::OEMMFF mmff;
    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'",mol.GetTitle());
      OEWriteMolecule(ofs, mol);
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESystem::OEThrow.Info("Molecule: %s Conformer: %n", mol.GetTitle(), conf->GetIdx()+1);
      conf->GetCoords(vecCoords);    

      double energy = mmff(vecCoords);
      OESystem::OEThrow.Info("Initial energy: %d kcal/mol", energy);

      OEOpt::OEBFGSOpt optimizer;
      energy = optimizer(mmff, vecCoords, vecCoords);
      OESystem::OEThrow.Info("Optimized energy: %d kcal/mol", energy);
      conf->SetCoords(vecCoords);
    }
    OEChem::OEWriteMolecule(ofs, mol);
  }

  return 0;
}

Energy components of single ligand in vacuum

The following example illustrates how to obtain various intermolecular energy components of a single ligand in vacuum, calculated using the OEMMFF force field.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following function is an example to show how to evaluate energies   //
//  of a small molecule and obtain various energy components.               //
//////////////////////////////////////////////////////////////////////////////

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;
  while (OEChem::OEReadMolecule(ifs, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    OEMolMMFF::OEMMFF mmff;
    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'",mol.GetTitle());
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    mol.GetCoords(vecCoords);
    OESystem::OEIter<OEOpt::OEFComponent> fcomp = mmff.GetFComponents(vecCoords);
    for (; fcomp; ++fcomp)
      OESystem::OEThrow.Info("Molecule: %s Component: %d Energy: %d kcal/mol", mol.GetTitle(), fcomp->name, fcomp->value);    
  }
  return 0;
}

MMFF interactions of single ligand

The following example illustrates how to obtain various interactions of a single ligand in the context of the OEMMFF force field.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following example demonstrates how to obtain interactions           //
//  of a small molecule in the context of MMFF force field.                 //
//////////////////////////////////////////////////////////////////////////////

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;
  while (OEChem::OEReadMolecule(ifs, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    OEMolMMFF::OEMMFF mmff;
    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'",mol.GetTitle());
      continue;
    }

    //print out interactions
    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    mol.GetCoords(vecCoords);
    OESystem::OEIter<OEMolPotential::OEInteraction> intc = mmff.GetInteractions(mol, vecCoords);
    OESystem::OEThrow.Info("Molecule: %s", mol.GetTitle());
    for(; intc; ++intc)
    {
      OEPlatform::OEMallocaPtr<double> vecGrads = OEMalloca(intc->NumValues());
      OESystem::OEThrow.Info("Interaction: %s Value: %d", intc->GetName(), intc->GetValues(vecGrads));
      for (OESystem::OEIter<OEChem::OEAtomBase>atom = intc->GetAtoms(); atom; ++atom)
        OESystem::OEThrow.Info("Atom index: %n", atom->GetIdx());
    }
  }

  return 0;
}

Protein-Ligand Complexes

Protein-ligand optimization with MMFF

The following example illustrates how to setup a ligand optimization within the context of a protein using the OEMMFF force field. Besides setting up the ligand MMFF interactions, protein preparation (OEForceField::PrepMol) followed by creation of the intermolecular interactions and addition of those to the force field are required to perform.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following example demonstrates how to perform ligand optimization   //
//  in the context of a protein using MMFF force field.                     //
//////////////////////////////////////////////////////////////////////////////

static void GetBoundsFromMol(double bounds[6],const OEChem::OEMolBase &mol,double pad)
{
  bool set=false;
  double c[3];
  OESystem::OEIter<OEChem::OEAtomBase> atom;
  for (atom = mol.GetAtoms();atom;++atom)
  {
    mol.GetCoords(atom,c);
    if (!set)
    {
      bounds[0] = bounds[3] = c[0];
      bounds[1] = bounds[4] = c[1];
      bounds[2] = bounds[5] = c[2];
      set = true;
    }
    else
    {
      if (c[0] < bounds[0]) bounds[0] = c[0];
      if (c[1] < bounds[1]) bounds[1] = c[1];
      if (c[2] < bounds[2]) bounds[2] = c[2];
      if (c[0] > bounds[3]) bounds[3] = c[0];
      if (c[1] > bounds[4]) bounds[4] = c[1];
      if (c[2] > bounds[5]) bounds[5] = c[2];
    }
  }

  bounds[0] -= pad;
  bounds[1] -= pad;
  bounds[2] -= pad;
  bounds[3] += pad;
  bounds[4] += pad;
  bounds[5] += pad;
}

int main(int argc, char *argv[])
{
  if (argc!=5)
    OESystem::OEThrow.Usage("%s <protein> <ligand> <bound-ligand> <output>", argv[0]);
    
  OEChem::oemolistream ips;
  if (!ips.open(argv[1]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading protein", argv[1]);

  OEChem::oemolistream ils;
  if (!ils.open(argv[2]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading ligands", argv[2]);

  OEChem::oemolistream ibs;
  if (!ibs.open(argv[3]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading bound ligand", argv[3]);

  OEChem::oemolostream ofs;
  if (!ofs.open(argv[4]))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", argv[4]);

  // Load the protein molecule and Prepare with MMFF
  OEChem::OEGraphMol protein;
  OEChem::OEReadMolecule(ips, protein);
  OEChem::OEAddExplicitHydrogens(protein);
  OEMolMMFF::OEMMFF mmff;
  mmff.PrepMol(protein);

  // Get bounds for protein ligand interaction
  double bounds[6];
  OEChem::OEGraphMol boundsMol;
  OEChem::OEReadMolecule(ibs, boundsMol);
  GetBoundsFromMol(bounds, boundsMol, 5.0);

  // Create intermolecular interaction functions and add to MMFF
  OEMolMMFF::OEMMFFParams params;
  OEMolMMFF::OEMMFFInterVdwNN interVdw(protein, params, bounds);
  OEMolMMFF::OEMMFFInterCoulomb interCoul(protein, bounds, 1.0, 1.0, 0.5);
  mmff.Add(interVdw);
  mmff.Add(interCoul);

  OEChem::OEMol mol; 
  while (OEChem::OEReadMolecule(ils, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'", mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESystem::OEThrow.Info("Molecule: %s Conformer: %n", mol.GetTitle(), conf->GetIdx()+1); 
      conf->GetCoords(vecCoords);

      // Calculate energy using adaptor
      double energy = mmff(vecCoords);
      OESystem::OEThrow.Info("Initial energy: %d kcal/mol", energy);

      // Optimize the ligand
      OEOpt::OEBFGSOpt optimizer;
      energy = optimizer(mmff, vecCoords, vecCoords);
      OESystem::OEThrow.Info("Optimized energy: %d kcal/mol", energy);
    }
    OEWriteMolecule(ofs, mol);
  }

  return 0;
}

Protein-ligand optimization with MMFFAmber

The following example illustrates how to setup a ligand optimization within the context of a protein using the OEMMFFAmber force field. Using an intermolecular forcefield like OEMMFFAmber for protein-ligand complexes is simpler compared to using OEMMFF, as the forcefield is already equipped to combine a ligand with a protein.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oeff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following example demonstrates how to perform ligand optimization   //
//  in the context of a protein using MMFFAmber force field.               //
//////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[])
{
  if (argc!=4)
    OESystem::OEThrow.Usage("%s <protein> <ligand> <output>", argv[0]);
    
  OEChem::oemolistream ips;
  if (!ips.open(argv[1]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading protein", argv[1]);

  OEChem::oemolistream ils;
  if (!ils.open(argv[2]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading ligands", argv[2]);

  OEChem::oemolostream ofs;
  if (!ofs.open(argv[3]))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", argv[3]);

  // Load the protein molecule and Prepare with MMFF
  OEChem::OEGraphMol protein;
  OEChem::OEReadMolecule(ips, protein);
  OEChem::OEAddExplicitHydrogens(protein);
  OEFF::OEMMFFAmber mmff(protein);
  mmff.PrepMol(protein);

  OEChem::OEMol mol; 
  while (OEChem::OEReadMolecule(ils, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'", mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESystem::OEThrow.Info("Molecule: %s Conformer: %n", mol.GetTitle(), conf->GetIdx()+1);     
      conf->GetCoords(vecCoords);

      // Calculate energy using adaptor
      double energy = mmff(vecCoords);
      OESystem::OEThrow.Info("Initial energy: %d kcal/mol", energy);

      // Optimize the ligand
      OEOpt::OEBFGSOpt optimizer;
      energy = optimizer(mmff, vecCoords, vecCoords);
      OESystem::OEThrow.Info("Optimized energy: %d kcal/mol", energy);
    }
    OEChem::OEWriteMolecule(ofs, mol);
  }

  return 0;
}

Using adaptors

Optimizing 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 OEMMFF 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.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following function is an example to show how to evaluate energy     //
//  and optimize a small molecule, with keeping a subset of atoms fixed,    //
//  using the MMFF force field.                                             //
//////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[])
{
  if (argc!=3)
    OESystem::OEThrow.Usage("%s <input> <output>", argv[0]);
    
  OEChem::oemolistream ifs;
  if (!ifs.open(argv[1]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", argv[1]);

  OEChem::oemolostream ofs;
  if (!ofs.open(argv[2]))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", argv[2]);

  OEChem::OEMol mol; 
  while (OEChem::OEReadMolecule(ifs, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    OEMolMMFF::OEMMFF mmff;
    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'",mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    // Setup adaptor. The first (false) means not to pass ownership of mmff,
    // and the second (false) means not to exclude interactions related
    // to the subset which would be fixed for calculations.
    OEMolPotential::OESubsetAdaptor adaptor(mmff, false, false);

    // Use a simple atoms predicate for the subset, followed by setup
    if (!adaptor.Set(OEChem::OEIsCarbon()) || !adaptor.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process subset for molecule: title = '%s'",mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESystem::OEThrow.Info("Molecule: %s Conformer: %n", mol.GetTitle(), conf->GetIdx()+1); 
      conf->GetCoords(vecCoords);

      // Get adaptor variables set corresponding to the coordinates
      OEPlatform::OEMallocaPtr<double> vecX = OEMalloca(adaptor.NumVar()*sizeof(double));
      adaptor.GetVar(vecX, vecCoords);

      // Calculate energy using adaptor
      double energy = adaptor(vecX);
      OESystem::OEThrow.Info("Initial energy: %d kcal/mol", energy);

      // Optimize the adaptor
      OEOpt::OEBFGSOpt optimizer;
      energy = optimizer(adaptor, vecX, vecX);
      OESystem::OEThrow.Info("Optimized energy: %d kcal/mol", energy);

      // Get optimized coordinates corresponding to the adaptor optimized variables
      adaptor.AdaptVar(vecCoords, vecX);
      conf->SetCoords(vecCoords);
    }
    OEWriteMolecule(ofs, mol);
  }

  return 0;
}

Optimizing single ligand with fixed torsions

The following example illustrates how to fix a set of torsions using the OETorAdaptor during a single ligand optimization in vacuum. The adaptor is initialized with the OEMMFF 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.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following function is an example to show how to evaluate energy     //
//  and optimize a small molecule, with keeping a set of torsions fixed,    //
//  using the MMFF force field.                                             //
//////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[])
{
  if (argc!=3)
    OESystem::OEThrow.Usage("%s <input> <output>", argv[0]);
    
  OEChem::oemolistream ifs;
  if (!ifs.open(argv[1]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading", argv[1]);

  OEChem::oemolostream ofs;
  if (!ofs.open(argv[2]))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", argv[2]);

  OEChem::OEMol mol; 
  while (OEReadMolecule(ifs, mol))
  {
    OEAddExplicitHydrogens(mol);

    OEMolMMFF::OEMMFF mmff;
    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'",mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    // Setup adaptor. The first (false) means not to pass ownership of mmff,
    // and the second (false) means not to exclude interactions related
    // to the subset which would be fixed for calculations.
    OEMolPotential::OETorAdaptor adaptor(mmff, false, false);

    // Use a simple atoms predicate for the torsions, followed by setup
    if (!adaptor.Set(OEChem::OEIsRotor()) || !adaptor.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process torsions for molecule: title = '%s'",mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESystem::OEThrow.Info("Molecule: %s Conformer: %n", mol.GetTitle(), conf->GetIdx()+1); 
      conf->GetCoords(vecCoords);

      // Get adaptor variables set corresponding to the coordinates
      OEPlatform::OEMallocaPtr<double> vecX = OEMalloca(adaptor.NumVar()*sizeof(double));
      adaptor.GetVar(vecX, vecCoords);

      // Calculate energy using adaptor
      double energy = adaptor(vecX);
      OESystem::OEThrow.Info("Initial energy: %d kcal/mol", energy);

      // Optimize the adaptor
      OEOpt::OEBFGSOpt optimizer;
      energy = optimizer(adaptor, vecX, vecX);
      OESystem::OEThrow.Info("Optimized energy: %d kcal/mol", energy);

      // Get optimized coordinates corresponding to the adaptor optimized variables
      adaptor.AdaptVar(vecCoords, vecX);
      conf->SetCoords(vecCoords);
    }
    OEChem::OEWriteMolecule(ofs, mol);
  }

  return 0;
}

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.

/**********************************************************************
Copyright (C) 2017
OpenEye Scientific Software, Inc.
***********************************************************************/

#include "openeye.h"
#include "oechem.h"
#include "oeopt.h"
#include "oemolpotential.h"
#include "oemmff.h"
#include "oeff.h"

//////////////////////////////////////////////////////////////////////////////
//  The following example demonstrates how to perform rigid optimization of //
//  a ligand in the context of a protein using MMFFAmber force field. The  //
//  OEQuatAdaptor it used for rigid optimization.                           //
//////////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[])
{
  if (argc!=4)
    OESystem::OEThrow.Usage("%s <protein> <ligand> <output>", argv[0]);
    
  OEChem::oemolistream ips;
  if (!ips.open(argv[1]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading protein", argv[1]);

  OEChem::oemolistream ils;
  if (!ils.open(argv[2]))
    OESystem::OEThrow.Fatal("Unable to open %s for reading ligands", argv[2]);

  OEChem::oemolostream ofs;
  if (!ofs.open(argv[3]))
    OESystem::OEThrow.Fatal("Unable to open %s for writing", argv[3]);

  // Load the protein molecule and Prepare with MMFF
  OEChem::OEGraphMol protein;
  OEChem::OEReadMolecule(ips, protein);
  OEChem::OEAddExplicitHydrogens(protein);
  OEFF::OEMMFFAmber mmff(protein);
  mmff.PrepMol(protein);

  OEChem::OEMol mol; 
  while (OEChem::OEReadMolecule(ils, mol))
  {
    OEChem::OEAddExplicitHydrogens(mol);

    if (!mmff.PrepMol(mol) || !mmff.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process molecule: title = '%s'", mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    // Setup adaptor. The first (false) means not to pass ownership of mmff,
    // and the second (false) means not to exclude interactions related
    // to the subset which would be fixed for calculations.
    OEMolPotential::OEQuatAdaptor adaptor(mmff, false, false);
    if (!adaptor.Setup(mol))
    {
      OESystem::OEThrow.Warning("Unable to process subset for molecule: title = '%s'",mol.GetTitle());
      OEChem::OEWriteMolecule(ofs, mol);
      continue;
    }

    OEPlatform::OEMallocaPtr<double> vecCoords = OEMalloca(3*mol.NumAtoms()*sizeof(double));
    for (OESystem::OEIter<OEChem::OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESystem::OEThrow.Info("Molecule: %s Conformer: %n", mol.GetTitle(), conf->GetIdx()+1); 
      conf->GetCoords(vecCoords);

      // Get adaptor variables set corresponding to the coordinates
      OEPlatform::OEMallocaPtr<double> vecX = OEMalloca(adaptor.NumVar()*sizeof(double));
      adaptor.GetVar(vecX, vecCoords);

      // Calculate energy using adaptor
      double energy = adaptor(vecX);
      OESystem::OEThrow.Info("Initial energy: %d kcal/mol", energy);

      // Optimize the ligand
      OEOpt::OEBFGSOpt optimizer;
      energy = optimizer(adaptor, vecX, vecX);
      OESystem::OEThrow.Info("Optimized energy: %d kcal/mol", energy);

      // Get optimized coordinates corresponding to the adaptor optimized variables
      adaptor.AdaptVar(vecCoords, vecX);
      conf->SetCoords(vecCoords);
    }
    OEChem::OEWriteMolecule(ofs, mol);
  }

  return 0;
}