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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff


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

    def begin(self):
        pass

    def NumVar(self):
        return 1

    def __call__(self, x, h=None, g=None):
        if isinstance(x, oechem.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 = oechem.OEDoubleArray(x, self.NumVar(), False)
            if g is not None:
                g1 = oechem.OEDoubleArray(g, self.NumVar(), False)
                g1[0] = 2.0*x1[0]-7.0
                h1 = oechem.OEDoubleArray(h, self.NumVar(), False)
                h1[0] = 2.0
                return True
            elif h is not None:
                h1 = oechem.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:
        oechem.OEThrow.Usage("%s <initial guess>" % args[0])

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

    func = Function()

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

    # Optimize function using BFGS Optimizer and checkpoint
    xOpt = oechem.OEDoubleArray(1)
    optimizer = oeff.OENewtonOpt()
    value = optimizer(func, x, xOpt)
    oechem.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))

Download code

Func2Opt.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff


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

    def begin(self):
        pass

    def NumVar(self):
        return 1

    def __call__(self, x, g=None):
        if isinstance(x, oechem.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 = oechem.OEDoubleArray(x, self.NumVar(), False)
            if g is not None:
                g1 = oechem.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(oeff.OECheckpoint0):
    def __init__(self):
        oeff.OECheckpoint0.__init__(self)
        pass

    def __call__(self, iteration, nvar, fval, var, state):
        # Intermediate information during optimization
        var1 = oechem.OEDoubleArray(var, 1, False)
        oechem.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:
        oechem.OEThrow.Usage("%s <initial guess>" % args[0])

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

    func = Function()

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

    # Optimize function using BFGS Optimizer and checkpoint
    xOpt = oechem.OEDoubleArray(1)
    optimizer = oeff.OEBFGSOpt()
    value = optimizer(func, ChkPt(), x, xOpt)
    oechem.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))

Download code

Func1Chk.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

# ///////////////////////////////////////////////////////////////////////////
# 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(oeff.OEMolFunc1):
    def __init__(self):
        oeff.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, oechem.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 = oechem.OEDoubleArray(coord, self.NumVar(), False)
            if grad is not None:
                grad1 = oechem.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:
        oechem.OEThrow.Usage("%s <input>" % args[0])

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

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

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

    return 0


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

Download code

Molfunc1Opt.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

# //////////////////////////////////////////////////////////////////////////
# The following example demonstrates how to perform ligand optimization   //
# //////////////////////////////////////////////////////////////////////////


class FFOptions(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "parsley")
        ffType.AddLegalValue("parsley")
        ffType.AddLegalValue("mmff94")
        ffType.AddLegalValue("mmff94s")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff94":
            return oeff.OEMMFF()
        elif ff == "mmff94s":
            return oeff.OEMMFF(oeff.OEMMFF94sParams())
        else:
            return oeff.OEParsley()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OESimpleAppOptions(ffOpts, "FFOptimize", oechem.OEFileStringType_Mol3D,
                                     oechem.OEFileStringType_Mol3D)
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    ff = ffOpts.GetFF()
    for mol in ifs.GetOEMols():
        print("Optimizing", mol.GetTitle())
        if not ff.PrepMol(mol) or not ff.Setup(mol):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

            # Calculate energy
            energy = ff(vecCoords)
            oechem.OEThrow.Info("Initial energy: %0.2f kcal/mol" % energy)

            # Optimize the ligand
            optimizer = oeff.OEBFGSOpt()
            energy = optimizer(ff, vecCoords, vecCoords)
            oechem.OEThrow.Info("Optimized energy: %0.2f kcal/mol" % energy)
            conf.SetCoords(vecCoords)

        oechem.OEWriteMolecule(ofs, mol)

    return 0


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

Download code

FFOptimize.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

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


class FFOptions(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "parsley")
        ffType.AddLegalValue("parsley")
        ffType.AddLegalValue("mmff94")
        ffType.AddLegalValue("mmff94s")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff94":
            return oeff.OEMMFF()
        elif ff == "mmff94s":
            return oeff.OEMMFF(oeff.OEMMFF94sParams())
        else:
            return oeff.OEParsley()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OESimpleAppOptions(ffOpts, "FFEnergy", oechem.OEFileStringType_Mol3D,
                                     oechem.OEFileStringType_Mol3D)
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    ff = ffOpts.GetFF()
    for mol in ifs.GetOEMols():
        if not ff.PrepMol(mol) or not ff.Setup(mol):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

            for fcomp in ff.GetFComponents(vecCoords):
                oechem.OEThrow.Info("Component: %s Energy: %0.2f kcal/mol"
                                    % (fcomp.name, fcomp.value))

        oechem.OEWriteMolecule(ofs, mol)

    return 0


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

Download code

FFEnergy.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

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


class FFOptions(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "parsley")
        ffType.AddLegalValue("parsley")
        ffType.AddLegalValue("mmff94")
        ffType.AddLegalValue("mmff94s")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff94":
            return oeff.OEMMFF()
        elif ff == "mmff94s":
            return oeff.OEMMFF(oeff.OEMMFF94sParams())
        else:
            return oeff.OEParsley()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OESimpleAppOptions(ffOpts, "FFInteractions", oechem.OEFileStringType_Mol3D,
                                     oechem.OEFileStringType_Mol3D)
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    ff = ffOpts.GetFF()
    for mol in ifs.GetOEMols():
        if not ff.PrepMol(mol) or not ff.Setup(mol):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

            for intc in ff.GetInteractions(mol, vecCoords):
                vecGrads = oechem.OEDoubleArray(intc.NumValues())
                oechem.OEThrow.Info("Interaction: %s Value: %0.2f"
                                    % (intc.GetName(), intc.GetValues(vecGrads)))

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

        oechem.OEWriteMolecule(ofs, mol)

    return 0


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

Download code

FFInteractions.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

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


class FFOptions(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "ff14sb_parsley")
        ffType.AddLegalValue("ff14sb_parsley")
        ffType.AddLegalValue("mmff_amber")
        ffType.AddLegalValue("mmff")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff":
            return oeff.OEMMFFComplex()
        elif ff == "mmff_amber":
            return oeff.OEMMFFAmberComplex()
        else:
            return oeff.OEFF14SBParsleyComplex()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OERefInputAppOptions(ffOpts, "FFComplexOpt", oechem.OEFileStringType_Mol3D,
                                       oechem.OEFileStringType_Mol3D, oechem.OEFileStringType_DU, "-protein")
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    rfs = oechem.oeifstream()
    if not rfs.open(opts.GetRefFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetRefFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    du = oechem.OEDesignUnit()
    if not oechem.OEReadDesignUnit(rfs, du):
        oechem.OEThrow.Fatal("Failed to read design unit")
    protein = oechem.OEGraphMol()
    du.GetProtein(protein)

    ff = ffOpts.GetFF()
    if not ff.SetupHost(protein):
        oechem.OEThrow.Fatal("Failed to setup protein: %s" % protein.GetTitle())

    for mol in ifs.GetOEMols():
        print("Optimizing", mol.GetTitle())
        if (not ff.Setup(mol)):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

            # Calculate energy
            energy = ff(vecCoords)
            oechem.OEThrow.Info("Initial energy: %0.2f kcal/mol" % energy)

            # Optimize the ligand
            optimizer = oeff.OEBFGSOpt()
            energy = optimizer(ff, vecCoords, vecCoords)
            oechem.OEThrow.Info("Optimized energy: %0.2f kcal/mol" % energy)
            conf.SetCoords(vecCoords)

        oechem.OEWriteMolecule(ofs, mol)

    return 0


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

Download code

FFComplexOpt.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

# //////////////////////////////////////////////////////////////////////////
# 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(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "parsley")
        ffType.AddLegalValue("parsley")
        ffType.AddLegalValue("mmff94")
        ffType.AddLegalValue("mmff94s")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff94":
            return oeff.OEMMFF()
        elif ff == "mmff94s":
            return oeff.OEMMFF(oeff.OEMMFF94sParams())
        else:
            return oeff.OEParsley()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OESimpleAppOptions(ffOpts, "FFSubset", oechem.OEFileStringType_Mol3D,
                                     oechem.OEFileStringType_Mol3D)
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    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.
    adaptor = oeff.OESubsetAdaptor(ff, False, False)

    # Use a simple atoms predicate for the subset
    adaptor.Set(oechem.OEIsCarbon())

    for mol in ifs.GetOEMols():
        print("Optimizing", mol.GetTitle())
        if not ff.PrepMol(mol) or not adaptor.Setup(mol):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

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

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

            # Optimize the adaptor
            optimizer = oeff.OEBFGSOpt()
            energy = optimizer(adaptor, vecX, vecX)
            oechem.OEThrow.Info("Optimized energy: %0.2f 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


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

Download code

FFSubset.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

# //////////////////////////////////////////////////////////////////////////
# The following is an example to show how to evaluate energy              //
# and optimize a set of torsions in a small molecule                      //
# //////////////////////////////////////////////////////////////////////////


class FFOptions(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "parsley")
        ffType.AddLegalValue("parsley")
        ffType.AddLegalValue("mmff94")
        ffType.AddLegalValue("mmff94s")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff94":
            return oeff.OEMMFF()
        elif ff == "mmff94s":
            return oeff.OEMMFF(oeff.OEMMFF94sParams())
        else:
            return oeff.OEParsley()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OESimpleAppOptions(ffOpts, "FFTorAdaptor", oechem.OEFileStringType_Mol3D,
                                     oechem.OEFileStringType_Mol3D)
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    ff = ffOpts.GetFF()

    # 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 = oeff.OETorAdaptor(ff, False, False)

    # Use a simple predicate for the subset of torsions to optimize
    adaptor.Set(oechem.OEIsRotor())

    for mol in ifs.GetOEMols():
        print("Optimizing", mol.GetTitle())
        if not ff.PrepMol(mol) or not adaptor.Setup(mol):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

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

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

            # Optimize the adaptor
            optimizer = oeff.OEBFGSOpt()
            energy = optimizer(adaptor, vecX, vecX)
            oechem.OEThrow.Info("Optimized energy: %0.2f 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


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

Download code

FFTorAdaptor.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

# //////////////////////////////////////////////////////////////////////////
# 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(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "ff14sb_parsley")
        ffType.AddLegalValue("ff14sb_parsley")
        ffType.AddLegalValue("mmff_amber")
        ffType.AddLegalValue("mmff")
        ffType.SetBrief("Force field type")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetFF(self):
        ff = self._fftype.GetStringValue()
        if ff == "":
            ff = self._fftype.GetStringDefault()

        if ff == "mmff":
            return oeff.OEMMFFComplex()
        elif ff == "mmff_amber":
            return oeff.OEMMFFAmberComplex()
        else:
            return oeff.OEFF14SBParsleyComplex()


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OERefInputAppOptions(ffOpts, "FFComplexQuat", oechem.OEFileStringType_Mol3D,
                                       oechem.OEFileStringType_Mol3D, oechem.OEFileStringType_DU, "-protein")
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    rfs = oechem.oeifstream()
    if not rfs.open(opts.GetRefFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetRefFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    du = oechem.OEDesignUnit()
    if not oechem.OEReadDesignUnit(rfs, du):
        oechem.OEThrow.Fatal("Failed to read design unit")
    protein = oechem.OEGraphMol()
    du.GetProtein(protein)

    ff = ffOpts.GetFF()
    if not ff.SetupHost(protein):
        oechem.OEThrow.Fatal("Failed to setup protein: %s" % protein.GetTitle())
    adap = oeff.OEQuatAdaptor(ff)

    for mol in ifs.GetOEMols():
        print("Optimizing", mol.GetTitle())

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

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

            # Calculate energy
            energy = adap(vecX)
            oechem.OEThrow.Info("Initial energy: %0.2f kcal/mol" % energy)

            # Optimize the ligand
            optimizer = oeff.OEBFGSOpt()
            energy = optimizer(adap, vecX, vecX)
            oechem.OEThrow.Info("Optimized energy: %0.2f kcal/mol" % energy)
            adap.AdaptVar(vecCoords, vecX)
            conf.SetCoords(vecCoords)

        oechem.OEWriteMolecule(ofs, mol)
    return 0


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

Download code

FFComplexQuat.py

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

#!/usr/bin/env python
# (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.

import sys
from openeye import oechem
from openeye import oeff

# //////////////////////////////////////////////////////////////////////////
# The following example demonstrates how to load a custom Smirnoff        //
# forcefield, and use that to optiiza a ligand.                           //
# //////////////////////////////////////////////////////////////////////////


class FFOptions(oechem.OEOptions):
    def __init__(self):
        oechem.OEOptions.__init__(self, "FFOptions")

        ffType = oechem.OEStringParameter("-ff", "")
        ffType.SetBrief("Force field XML file")
        self._fftype = self.AddParameter(ffType)
        pass

    def CreateCopy(self):
        return self

    def GetXML(self):
        ffxml = self._fftype.GetStringValue()
        if ffxml == "":
            ffxml = self._fftype.GetStringDefault()
        return ffxml


def main(argv=[__name__]):
    ffOpts = FFOptions()
    opts = oechem.OESimpleAppOptions(ffOpts, "FFCustomSmirnoff", oechem.OEFileStringType_Mol3D,
                                     oechem.OEFileStringType_Mol3D)
    if oechem.OEConfigureOpts(opts, argv, False) == oechem.OEOptsConfigureStatus_Help:
        return 0
    ffOpts.UpdateValues(opts)

    ifs = oechem.oemolistream()
    if not ifs.open(opts.GetInFile()):
        oechem.OEThrow.Fatal("Unable to open %s for reading" % opts.GetInFile())

    ofs = oechem.oemolostream()
    if not ofs.open(opts.GetOutFile()):
        oechem.OEThrow.Fatal("Unable to open %s for writing" % opts.GetOutFile())

    params = oeff.OESmirnoffParams()
    if ffOpts.GetXML() != "":
        if not params.Load(ffOpts.GetXML()):
            oechem.OEThrow.Fatal("Unable to load force field parameters: %s" % ffOpts.GetXML())
    ff = oeff.OESmirnoff(params)

    for mol in ifs.GetOEMols():
        print("Optimizing", mol.GetTitle())
        if not ff.PrepMol(mol) or not ff.Setup(mol):
            oechem.OEThrow.Warning("Unable to process molecule: title = '%s'" % mol.GetTitle())
            oechem.OEWriteMolecule(ofs, mol)
            continue

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

            # Calculate energy
            energy = ff(vecCoords)
            oechem.OEThrow.Info("Initial energy: %0.2f kcal/mol" % energy)

            # Optimize the ligand
            optimizer = oeff.OEBFGSOpt()
            energy = optimizer(ff, vecCoords, vecCoords)
            oechem.OEThrow.Info("Optimized energy: %0.2f kcal/mol" % energy)
            conf.SetCoords(vecCoords)

        oechem.OEWriteMolecule(ofs, mol)

    return 0


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

Download code

FFCustomSmirnoff.py