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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *


#////////////////////////////////////////////////////////////////////////////
#  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(OEFunc2):
    def __init__(self):
        OEFunc2.__init__(self)
        pass

    def begin(self):
        pass

    def NumVar(self):
        return 1

    def __call__(self, x, h=None, g=None):
        if isinstance(x, OEDoubleArray):
            if g is not None:
                g[0] = 2.0*x[0]-7.0 # gradient
                h[0] = 2.0  # hessien
                return True
            elif h is not None:
                h[0] = 2.0*x[0]-7.0 #gradient
                return x[0]*x[0]-7*x[0]+63
            else:
                return x[0]*x[0]-7*x[0]+63
        else:
            x1 = OEDoubleArray(x, self.NumVar(), False)
            if g is not None:
                g1 = OEDoubleArray(g, self.NumVar(), False)
                g1[0] = 2.0*x1[0]-7.0
                h1 = OEDoubleArray(h, self.NumVar(), False)
                h1[0] = 2.0
                return True
            elif h is not None:
                h1 = OEDoubleArray(h, self.NumVar(), False)
                h1[0] = 2.0*x1[0]-7.0
                return x1[0]*x1[0]-7*x1[0]+63
            else:
                return x1[0]*x1[0]-7*x1[0]+63

def main(args):
    if len(args) != 2:
        OEThrow.Usage("%s <initial guess>" % args[0])

    x = OEDoubleArray(1)
    try:
        x[0] = float(args[1])
    except ValueError:
        OEThrow.Usage("%s <initial guess (expecting a number)>" % args[0])

    func = Function()

    # Calculate function value at given x
    value = func(x);
    OEThrow.Info("Function value at x = %d is %d" % (x[0], value))

    # Optimize function using BFGS Optimizer and checkpoint
    xOpt = OEDoubleArray(1)
    optimizer = OENewtonOpt()
    value = optimizer(func, x, xOpt);
    OEThrow.Info("Function has a minimia at x = %d and the minimum value is %d" % (xOpt[0], value))

    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *


#////////////////////////////////////////////////////////////////////////////
#  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(OEFunc1):
    def __init__(self):
        OEFunc1.__init__(self)
        pass

    def begin(self):
        pass

    def NumVar(self):
        return 1

    def __call__(self, x, g=None):
        if isinstance(x, OEDoubleArray):
            if g is not None:
                g[0] = 2.0*x[0]-7.0
                return x[0]*x[0]-7*x[0]+63
            else:
                return x[0]*x[0]-7*x[0]+63
        else:
            x1 = OEDoubleArray(x, self.NumVar(), False)
            if g is not None:
                g1 = OEDoubleArray(g, self.NumVar(), False)
                g1[0] = 2.0*x1[0]-7.0
                return x1[0]*x1[0]-7*x1[0]+63
            else:
                return x1[0]*x1[0]-7*x1[0]+63

class ChkPt(OECheckpoint0):
    def __init__(self):
        OECheckpoint0.__init__(self)
        pass

    def __call__(self, iteration, nvar, fval, var, state):
        # Intermediate information during optimization
        var1 = OEDoubleArray(var, 1, False)
        OEThrow.Info("iteration: %d x: %d value: %d state: %d" % (iteration, var1[0], fval, state))

        # To demonstrate how to force quitting optimization from checkpoint
        # this returns false when iteration is 5
        if iteration >= 5:
          return False
        else:
          return True


def main(args):
    if len(args) != 2:
        OEThrow.Usage("%s <initial guess>" % args[0])

    x = OEDoubleArray(1)
    try:
        x[0] = float(args[1])
    except ValueError:
        OEThrow.Usage("%s <initial guess (expecting a number)>" % args[0])

    func = Function()

    # Calculate function value at given x
    value = func(x);
    OEThrow.Info("Function value at x = %d is %d" % (x[0], value))

    # Optimize function using BFGS Optimizer and checkpoint
    xOpt = OEDoubleArray(1)
    optimizer = OEBFGSOpt()
    value = optimizer(func, ChkPt(), x, xOpt);
    OEThrow.Info("Function has a minimia at x = %d and the minimum value is %d" % (xOpt[0], value))

    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

#////////////////////////////////////////////////////////////////////////////
#  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(OEMolFunc1):
    def __init__(self):
        OEMolFunc1.__init__(self)
        self.natoms = 0
        pass

    def begin(self):
        pass

    def NumVar(self):
        return 3*self.natoms

    def Setup(self, mol):
        self.natoms = mol.GetMaxAtomIdx()
        return True

    def __call__(self, coord, grad=None):
        energy = 0.0
        if isinstance(coord, OEDoubleArray):
            if grad is not None:
                for i in range(0, self.natoms):
                    grad[3*i  ] += 2.0 * coord[3*i  ]; #x gradient
                    grad[3*i+1] += 2.0 * coord[3*i+1]; #y gradient
                    grad[3*i+2] += 2.0 * coord[3*i+2]; #z gradient
                    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
            else:
                for i in range(0, self.natoms):
                    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
        else:
            coord1 = OEDoubleArray(coord, self.NumVar(), False)
            if grad is not None:
                grad1 = OEDoubleArray(grad, self.NumVar(), False)
                for i in range(0, self.natoms):
                    grad1[3*i  ] += 2.0 * coord1[3*i  ]; #x gradient
                    grad1[3*i+1] += 2.0 * coord1[3*i+1]; #y gradient
                    grad1[3*i+2] += 2.0 * coord1[3*i+2]; #z gradient
                    energy += coord1[3*i  ] * coord1[3*i  ] #x distance from zero
                    energy += coord1[3*i+1] * coord1[3*i+1] #y distance from zero
                    energy += coord1[3*i+2] * coord1[3*i+2] #z distance from zero
            else:
                for i in range(0, self.natoms):
                    energy += coord1[3*i  ] * coord1[3*i  ] #x distance from zero
                    energy += coord1[3*i+1] * coord1[3*i+1] #y distance from zero
                    energy += coord1[3*i+2] * coord1[3*i+2] #z distance from zero
        return energy

def main(args):
    if len(args) != 2:
        OEThrow.Usage("%s <input>" % args[0])

    ifs = oemolistream()
    if not ifs.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading" % args[1])

    mol = OEGraphMol()
    OEReadMolecule(ifs, mol)
    vecCoords = OEDoubleArray(3*mol.NumAtoms())
    mol.GetCoords(vecCoords)

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

    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

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

def main(args):
    if len(args) != 3:
        OEThrow.Usage("%s <input> <output>" % args[0])

    ifs = oemolistream()
    if not ifs.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading" % args[1])

    ofs = oemolostream()
    if not ofs.open(args[2]):
        OEThrow.Fatal("Unable to open %s for writing" % args[2])

    mol = OEMol()
    while OEReadMolecule(ifs, mol):
        OEAddExplicitHydrogens(mol);

        mmff = OEMMFF()
        if not mmff.PrepMol(mol) or not mmff.Setup(mol):
            OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            OEWriteMolecule(ofs, mol)
            continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        for conf in mol.GetConfs():
            OEThrow.Info("Molecule: %s Conformer: %d" % (mol.GetTitle(), conf.GetIdx()+1)) 
            conf.GetCoords(vecCoords)

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

            optimizer = OEBFGSOpt()
            energy = optimizer(mmff, vecCoords, vecCoords)
            OEThrow.Info("Optimized energy: %d kcal/mol" % energy)
            conf.SetCoords(vecCoords);

        OEWriteMolecule(ofs, mol);

    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

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

def main(args):
    if len(args) != 2:
        OEThrow.Usage("%s <input>" % args[0])

    ifs = oemolistream()
    if not ifs.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading" % args[1])

    mol = OEGraphMol()
    while OEReadMolecule(ifs, mol):
        OEAddExplicitHydrogens(mol)

        mmff = OEMMFF()
        if not mmff.PrepMol(mol) or not mmff.Setup(mol):
          OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
          continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        mol.GetCoords(vecCoords)
        for fcomp in mmff.GetFComponents(vecCoords):
          OEThrow.Info("Molecule: %s Component: %s Energy: %d kcal/mol" % (mol.GetTitle(), fcomp.name, fcomp.value))  
  
    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

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

def main(args):
    if len(args) != 2:
        OEThrow.Usage("%s <input>" % args[0])

    ifs = oemolistream()
    if not ifs.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading" % args[1])

    mol = OEGraphMol()
    while OEReadMolecule(ifs, mol):
        OEAddExplicitHydrogens(mol);

        mmff = OEMMFF()
        if not mmff.PrepMol(mol) or not mmff.Setup(mol):
            OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            continue

        #print out interactions
        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        mol.GetCoords(vecCoords)
        OEThrow.Info("Molecule: %s" % mol.GetTitle())
        for intc in mmff.GetInteractions(mol, vecCoords):
            vecGrads = OEDoubleArray(intc.NumValues())
            OEThrow.Info("Interaction: %s Value: %d" % (intc.GetName(), intc.GetValues(vecGrads)))

            for atom in intc.GetAtoms():
                OEThrow.Info("Atom index: %d" % atom.GetIdx())            

    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

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

def GetBoundsFromMol(bounds, mol, pad):
    set = False
    c = OEDoubleArray(3)
    for atom in mol.GetAtoms():
        mol.GetCoords(atom,c);
        if not 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


def main(args):
    if len(args) != 5:
        OEThrow.Usage("%s <protein> <ligand> <bound-ligand> <output>" % args[0])

    ips = oemolistream()
    if not ips.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading protein" % args[1])

    ils = oemolistream()
    if not ils.open(args[2]):
        OEThrow.Fatal("Unable to open %s for reading ligand" % args[2])

    ibs = oemolistream()
    if not ibs.open(args[3]):
        OEThrow.Fatal("Unable to open %s for reading bound ligand" % args[3])

    ofs = oemolostream()
    if not ofs.open(args[4]):
        OEThrow.Fatal("Unable to open %s for writing" % args[4])


    # Load the protein molecule and Prepare with MMFF
    protein = OEGraphMol()
    OEReadMolecule(ips, protein)
    OEAddExplicitHydrogens(protein)
    mmff = OEMMFF()
    mmff.PrepMol(protein)

    # Get bounds for protein ligand interaction
    bounds = OEDoubleArray(6)
    boundsMol = OEGraphMol()
    OEReadMolecule(ibs, boundsMol)
    GetBoundsFromMol(bounds, boundsMol, 5.0)

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

    mol = OEMol() 
    while OEReadMolecule(ils, mol):
        OEAddExplicitHydrogens(mol)

        if (not mmff.PrepMol(mol)) or  (not mmff.Setup(mol)):
            OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            OEWriteMolecule(ofs, mol)
            continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        for conf in mol.GetConfs():
            OEThrow.Info("Molecule: %s Conformer: %d" % (mol.GetTitle(), conf.GetIdx()+1))
            conf.GetCoords(vecCoords)

            # Calculate energy
            energy = mmff(vecCoords)
            OEThrow.Info("Initial energy: %d kcal/mol" % energy)

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

    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

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

def main(args):
    if len(args) != 4:
        OEThrow.Usage("%s <protein> <ligand> <output>" % args[0])

    ips = oemolistream()
    if not ips.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading protein" % args[1])

    ils = oemolistream()
    if not ils.open(args[2]):
        OEThrow.Fatal("Unable to open %s for reading ligand" % args[2])

    ofs = oemolostream()
    if not ofs.open(args[3]):
        OEThrow.Fatal("Unable to open %s for writing" % args[3])

    # Load the protein molecule and Prepare with MMFF
    protein = OEGraphMol()
    OEReadMolecule(ips, protein)
    OEAddExplicitHydrogens(protein)
    mmff = OEMMFFAmber(protein)
    mmff.PrepMol(protein)

    mol = OEMol()
    while OEReadMolecule(ils, mol):
        OEAddExplicitHydrogens(mol)
        if (not mmff.PrepMol(mol)) or (not mmff.Setup(mol)):
            OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            OEWriteMolecule(ofs, mol);
            continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        for  conf in mol.GetConfs():
            OEThrow.Info("Molecule: %s Conformer: %d" % (mol.GetTitle(), conf.GetIdx()+1))
            conf.GetCoords(vecCoords)

            # Calculate energy
            energy = mmff(vecCoords);
            OEThrow.Info("Initial energy: %d kcal/mol" % energy)

            # Optimize the ligand
            optimizer = OEBFGSOpt()
            energy = optimizer(mmff, vecCoords, vecCoords)
            OEThrow.Info("Optimized energy: %d kcal/mol" % energy)
            conf.SetCoords(vecCoords)
    
        OEWriteMolecule(ofs, mol)
    return 0

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

#////////////////////////////////////////////////////////////////////////////
#  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.                                             //
#////////////////////////////////////////////////////////////////////////////

def main(args):
    if len(args) != 3:
        OEThrow.Usage("%s <input> <output>" % args[0])

    ifs = oemolistream()
    if not ifs.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading" % args[1])

    ofs = oemolostream()
    if not ofs.open(args[2]):
        OEThrow.Fatal("Unable to open %s for writing" % args[2])

    mol = OEMol() 
    while OEReadMolecule(ifs, mol):
        OEAddExplicitHydrogens(mol)

        mmff = OEMMFF()
        if not mmff.PrepMol(mol) or not mmff.Setup(mol):
          OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
          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.
        adaptor = OESubsetAdaptor(mmff, False, False)

        # Use a simple atoms predicate for the subset, followed by setup
        if (not adaptor.Set(OEIsCarbon())) or (not adaptor.Setup(mol)):
            OEThrow.Warning("Unable to process subset for molecule: title = '%s'" % mol.GetTitle())
            OEWriteMolecule(ofs, mol)
            continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        for conf in mol.GetConfs():
            OEThrow.Info("Molecule: %s Conformer: %d" % (mol.GetTitle(), conf.GetIdx()+1))
            conf.GetCoords(vecCoords)

            # Get adaptor variables set corresponding to the coordinates
            vecX = OEDoubleArray(adaptor.NumVar())
            adaptor.GetVar(vecX, vecCoords)

            # Calculate energy using adaptor
            energy = adaptor(vecX)
            OEThrow.Info("Initial energy: %d kcal/mol" % energy)

            # Optimize the adaptor
            optimizer = OEBFGSOpt()
            energy = optimizer(adaptor, vecX, vecX)
            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

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

#////////////////////////////////////////////////////////////////////////////
#  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.                                             //
#////////////////////////////////////////////////////////////////////////////

def main(args):
    if len(args) != 3:
        OEThrow.Usage("%s <input> <output>" % args[0])

    ifs = oemolistream()
    if not ifs.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading" % args[1])

    ofs = oemolostream()
    if not ofs.open(args[2]):
        OEThrow.Fatal("Unable to open %s for writing" % args[2])

    mol = OEMol() 
    while OEReadMolecule(ifs, mol):
        OEAddExplicitHydrogens(mol)

        mmff = OEMMFF()
        if not mmff.PrepMol(mol) or not mmff.Setup(mol):
          OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
          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.
        adaptor = OETorAdaptor(mmff, False, False)

        # Use a simple atoms predicate for the subset, followed by setup
        if (not adaptor.Set(OEIsRotor())) or (not adaptor.Setup(mol)):
            OEThrow.Warning("Unable to process torsions for molecule: title = '%s'" % mol.GetTitle())
            OEWriteMolecule(ofs, mol)
            continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        for conf in mol.GetConfs():
            OEThrow.Info("Molecule: %s Conformer: %d" % (mol.GetTitle(), conf.GetIdx()+1))
            conf.GetCoords(vecCoords)

            # Get adaptor variables set corresponding to the coordinates
            vecX = OEDoubleArray(adaptor.NumVar())
            adaptor.GetVar(vecX, vecCoords)

            # Calculate energy using adaptor
            energy = adaptor(vecX)
            OEThrow.Info("Initial energy: %d kcal/mol" % energy)

            # Optimize the adaptor
            optimizer = OEBFGSOpt()
            energy = optimizer(adaptor, vecX, vecX)
            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

if __name__ == "__main__":
    sys.exit(main(sys.argv))

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.

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2017 OpenEye Scientific Software, Inc.
#############################################################################
import sys

from openeye.oechem import *
from openeye.oeff import *

#////////////////////////////////////////////////////////////////////////////
#  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.                           //
#////////////////////////////////////////////////////////////////////////////

def main(args):
    if len(args) != 4:
        OEThrow.Usage("%s <protein> <ligand> <output>" % args[0])

    ips = oemolistream()
    if not ips.open(args[1]):
        OEThrow.Fatal("Unable to open %s for reading protein" % args[1])

    ils = oemolistream()
    if not ils.open(args[2]):
        OEThrow.Fatal("Unable to open %s for reading ligand" % args[2])

    ofs = oemolostream()
    if not ofs.open(args[3]):
        OEThrow.Fatal("Unable to open %s for writing" % args[3])

    # Load the protein molecule and Prepare with MMFF
    protein = OEGraphMol()
    OEReadMolecule(ips, protein)
    OEAddExplicitHydrogens(protein)
    mmff = OEMMFFAmber(protein)
    mmff.PrepMol(protein)

    mol = OEMol()
    while OEReadMolecule(ils, mol):
        OEAddExplicitHydrogens(mol);
        if (not mmff.PrepMol(mol)) or (not mmff.Setup(mol)):
            OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            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.
        adaptor = OEQuatAdaptor(mmff, False, False)
        if not adaptor.Setup(mol):
            OEThrow.Warning("Unable to process subset for molecule: title = '%s'" % mol.GetTitle())
            OEWriteMolecule(ofs, mol)
            continue

        vecCoords = OEDoubleArray(3*mol.NumAtoms())
        for  conf in mol.GetConfs():
            OEThrow.Info("Molecule: %s Conformer: %d" % (mol.GetTitle(), conf.GetIdx()+1))
            conf.GetCoords(vecCoords)

            # Get adaptor variables set corresponding to the coordinates
            vecX = OEDoubleArray(adaptor.NumVar())
            adaptor.GetVar(vecX, vecCoords);           

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

            # Optimize the ligand
            optimizer = OEBFGSOpt()
            energy = optimizer(adaptor, vecX, vecX)
            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

if __name__ == "__main__":
    sys.exit(main(sys.argv))