Working with Szybki TK

The basic Szybki TK API is provided to the users in the OESzybki, OESzybkiOptions, OESzybkiResults and OESzybkiEnsembleResults classes. In addition, when the purpose is to calculate such compound properties like solvation free energy or free energy of selecting specific conformations out of the ensemble, the higher level API defined in the free functions OEEstimateSolvFreeEnergy and OEEstimateConfFreeEnergies.

Ligand Energetics and Optimization

Single ligand in vacuum

The following example illustrates how to optimize a single ligand in vacuum. As one can see assuming a molecule is successfully read, only two objects are needed to perform the optimization: OESzybki and OESzybkiResults. The latter is passed as the second parameter to the parenthesis operator OESzybki.operator(). Molecule with the optimized coordinates is returned as a first parameter. Final energy results are available as SD tags in the returned molecule and optionally with a call of a method OESzybkiResults.Print.

/******************************************************************************
 * Copyright 2004-2017 OpenEye Scientific Software, Inc.
 *****************************************************************************/
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class SimpleVacuum {

    public static void main(String[] args) {
        if (args.length!=2) {
            oechem.OEThrow.Usage("SimpleVacuum <molfile> <outfile>");
        }

        oemolistream ifs = new oemolistream();
        oemolostream ofs = new oemolostream();

        if (!ifs.open(args[0])) {
            oechem.OEThrow.Fatal("Unable to open " + args[0] + " for reading");
        }
        if (!ofs.open(args[1])) {
            oechem.OEThrow.Fatal("Unable to open " + args[1] + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);

        OESzybkiOptions opts = new OESzybkiOptions();
        OESzybki sz = new OESzybki(opts);
        OESzybkiResults res = new OESzybkiResults();

        if (!sz.call(mol, res)) {
            return;
        }

        oechem.OEWriteMolecule(ofs, mol);
        ofs.close();
        res.Print(oechem.getOeout());

        res.delete();
        sz.delete();
    }
}

Optimization of a set of ligands

The next example illustrates the usage of Szybki TK to optimize a set of compounds with the MMFF94 force field in vacuum or in solution using Sheffield solvation model. Optionally attractive VdW can be removed. The optimization is done by default in full Cartesian coordinates, however torsion space optimization or single point calculation could be done too. A group of atoms which belong to specified SMARTS pattern might be excluded from optimization so their positions will be fixed at their initial coordinates. Note that the OESzybki object is made with a constructor which takes the instance of the OESzybkiOptions class.

//*****************************************************************************
//*  Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class LigandMultipleOptions
{
    public static void main(String argv[])
    {
        OEInterface itf = new OEInterface();
        if (!SetupInterface("LigandMultipleOptions", itf, argv))
            return;

        oemolistream ifs = new oemolistream();
        if (!ifs.open(itf.GetString("-in")))
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-in") + " for reading");

        oemolostream ofs = new oemolostream();
        if (!ofs.open(itf.GetString("-out")))
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-out") + " for writing");

        oeostream logfile = oechem.oeout;
        if (itf.HasString("-log")) {
            if (!logfile.open(itf.GetString("-log")))
                oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-log") + " for writing");
        }

        // Szybki options
        OESzybkiOptions opts = new OESzybkiOptions();

        // select run type
        if (itf.GetBool("-t"))
            opts.SetRunType(OERunType.TorsionsOpt);

        if (itf.GetBool("-n"))
            opts.SetRunType(OERunType.SinglePoint);

        // apply solvent model
        if (itf.GetBool("-s"))
            opts.GetSolventOptions().SetSolventModel(OESolventModel.Sheffield);

        // remove attractive VdW forces
        if (itf.GetBool("-a"))
            opts.GetGeneralOptions().SetRemoveAttractiveVdWForces(true);

        // Szybki object;
        OESzybki sz = new OESzybki(opts);

        // fix atoms
        if (itf.HasString("-f"))
            if (!sz.FixAtoms(itf.GetString("-f")))
                oechem.OEThrow.Warning("Failed to fix atoms for " + itf.GetString("-f"));

        // process molecules
        OEMol mol = new OEMol();
        while (oechem.OEReadMolecule(ifs, mol)) {
            logfile.write("\nMolecule " + mol.GetTitle() + "\n");
            boolean no_res = true;
            for (OESzybkiResults results : sz.call(mol)) {
                results.Print(logfile);
                no_res = false;
            }

            if (no_res) {
                oechem.OEThrow.Warning("No results processing molecule: " + mol.GetTitle());
                continue;
            }
            oechem.OEWriteMolecule(ofs, mol);
        }
    }

    private static boolean SetupInterface(String progName, OEInterface itf, String argv[])
    {
        oechem.OEConfigure(itf, interfaceData);
        if (oechem.OECheckHelp(itf, argv, progName, false)) {
            return false;
        }
        if (!oechem.OEParseCommandLine(itf, argv, progName)) {
            return false;
        }
        if (!oechem.OEIsReadable(oechem.OEGetFileType(oechem.OEGetFileExtension(itf.GetString("-in"))))) {
            oechem.OEThrow.Warning(itf.GetString("-in") + " is not a readable input file");
            return false;
        }
        if (!oechem.OEIsWriteable(oechem.OEGetFileType(oechem.OEGetFileExtension(itf.GetString("-out"))))) {
            oechem.OEThrow.Warning(itf.GetString("-out") + " is not a writable output file");
            return false;
        }
        return true;
    }

    private static String interfaceData = 
"!PARAMETER -in\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Input molecule file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -out\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Output molecule file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -log\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF Log file name. Defaults to standard out.\n" +
"!END\n" +
"\n" +
"!PARAMETER -s\n" +
"!TYPE bool\n" +
"!DEFAULT false\n" +
"!REQUIRED false\n" +
"!BRIEF Optimization in solution.\n" +
"!END\n" +
"\n" +
"!PARAMETER -t\n" +
"!TYPE bool\n" +
"!DEFAULT false\n" +
"!REQUIRED false\n" +
"!BRIEF Optimization of torsions.\n" +
"!END\n" +
"\n" +
"!PARAMETER -a\n" +
"!TYPE bool\n" +
"!DEFAULT false\n" +
"!REQUIRED false\n" +
"!BRIEF No attractive VdW forces.\n" +
"!END\n" +
"\n" +
"!PARAMETER -n\n" +
"!TYPE bool\n" +
"!DEFAULT false\n" +
"!REQUIRED false\n" +
"!BRIEF Single point calculation.\n" +
"!END\n" +
"\n" +
"!PARAMETER -f\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF SMARTS pattern of fixed atoms.\n" +
"!END";
}

Optimization of a single ligand with the Newton-Raphson method

The next example in this section shows how to use Newton-Raphson optimization method, rather than the default BFGS:

//*****************************************************************************
// Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class SimpleNewton
{
    public static void main(String[] args) {
        if (args.length != 2)
            oechem.OEThrow.Usage("SimpleNewton input_molecule output_molecule");

        String inputFile  = args[0];
        String outputFile = args[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile))
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile))
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OESzybkiOptions opts = new OESzybkiOptions();
        opts.GetOptOptions().SetOptimizerType(OEOptType.NEWTON);
        opts.GetSolventOptions().SetSolventModel(OESolventModel.Sheffield);

        OESzybki sz = new OESzybki(opts);
        OESzybkiResults res = new OESzybkiResults();
        if (sz.call(mol, res)) {
            oechem.OEWriteMolecule(ofs, mol);
            res.Print(oechem.oeout);
        }
        ofs.close();
    }
}

Optimization of all conformers of a ligand

Finally, the last example in this section shows how to use Newton-Raphson optimization method on all conformers of a ligand. The current charges of the ligand are used and will not be changed during the optimization.

//*****************************************************************************
// Copyright (C) 2015-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oequacpac.*;
import openeye.oeszybki.*;

public class SimpleConformer
{
    public static void main(String[] args) {
        if (args.length != 2)
            oechem.OEThrow.Usage("SimpleConformer input_molecule output_molecule");

        String inputFile  = args[0];
        String outputFile = args[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile))
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile))
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");

        OESzybkiOptions opts = new OESzybkiOptions();
        opts.GetOptOptions().SetOptimizerType(OEOptType.NEWTON);
        opts.GetGeneralOptions().SetForceFieldType(OEForceFieldType.MMFF94S);
        opts.GetSolventOptions().SetSolventModel(OESolventModel.Sheffield);
        opts.GetSolventOptions().SetChargeEngine(new OEChargeEngineNoOp());

        OESzybki sz = new OESzybki(opts);
        OESzybkiResults res = new OESzybkiResults();
        OEMol mol = new OEMol();
        while(oechem.OEReadMolecule(ifs, mol))
        {
            for(OEConfBase conf : mol.GetConfs())
            {
                if (sz.call(conf, res))
                {
                    oechem.OESetSDData(conf, 
                        new OESDDataPair("Total_energy", 
                        String.valueOf(res.GetTotalEnergy())));
                }
            }
            oechem.OEWriteMolecule(ofs, mol);
        }
        
        ifs.close();
        ofs.close();
    }
}

Protein-Ligand Energetics and Optimization

Examples in this section show how to optimize a bound ligand.

Optimization of a single bound ligand

The simplest case is illustrated below. Notice the usage of the method OESzybki.SetProtein which tells Szybki that the ligand is placed inside the protein. Since no protein-ligand electrostatics have been specified, neither the coordinates types which should be used by the optimizer, the code below performs the optimization for a rigid ligand using 6 translational-rotational coordinates in the MMFF94 VdW potential field.

//*****************************************************************************
//*  Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class SimpleProtein {
    public static void main(String argv[]) {
        if (argv.length != 3) {
            oechem.OEThrow.Usage("SimpleProtein <molfile> <protein> <outfile>");
        }

        String inputLFile = argv[0];
        String inputPFile = argv[1];
        String outputFile = argv[2];

        oemolistream lfs = new oemolistream();
        if (!lfs.open(inputLFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputLFile + " for reading");
        }

        oemolistream pfs = new oemolistream();
        if (!pfs.open(inputPFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputPFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEGraphMol mol = new OEGraphMol();
        oechem.OEReadMolecule(lfs, mol);
        lfs.close();

        OEGraphMol protein = new OEGraphMol();
        oechem.OEReadMolecule(pfs, protein);
        pfs.close();

        OESzybkiOptions opts = new OESzybkiOptions();
        OESzybki sz = new OESzybki(opts);
        sz.SetProtein(protein);
        OESzybkiResults res = new OESzybkiResults();
        if (sz.call(mol, res)) {
            oechem.OEWriteMolecule(ofs, mol);
            res.Print(oechem.oeout);
        }
        ofs.close();
    }
}

Optimization of a set of bound ligands in a rigid receptor

The next example illustrates the usage of Szybki TK to optimize a set of ligands with MMFF94 force field inside a protein receptor. By default only VdW protein-ligand interaction is used. Optionally exact or grid Coulomb potential as well as PB solvent screening potentials could be added. When grid potential is selected (either Coulomb or PB) optionally it could be saved or read in when the corresponding grid file is present in the specified directory. Notice that when the exact Coulomb electrostatics is chosen, also the exact VdW potential is chosen (method OESzybkiProteinOptions.SetExactVdWProteinLigand) which allows for tight gradients convergence (methods OESzybkiOptOptions.SetMaxIter and OESzybkiOptOptions.SetGradTolerance). By default ligand is treated as a solid body, that is only its translational and rotational degrees of freedom are optimized. Optionally also torsional degrees could be optimized. In this example protein receptor is rigid. Molecular input file should contain initial 3D coordinates of molecules in any format supported by OEChem TK. Output file is specified with the -out flag. In addition a log file containing energy data terms values is written to stdout or a file specified by -log.

/******************************************************************************
 * Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
 *****************************************************************************/
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class ProteinMultipleOptions
{
    public static void main(String argv[])
    {
        OEInterface itf = new OEInterface ();
        oechem.OEConfigure(itf, interfaceData);
        oechem.OEParseCommandLine (itf, argv, "ProteinMultipleOptions");

        oemolistream lfs = new oemolistream();
        if (!lfs.open(itf.GetString("-in"))) {
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-in") + " for reading");
        }

        oemolistream pfs = new oemolistream();
        if (!pfs.open(itf.GetString("-p"))) {
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-p") + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(itf.GetString("-out"))) {
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-out") + " for writing");
        }

        oeostream logfile = oechem.oeout;
        if (itf.HasString("-log")) {
            if (!logfile.open(itf.GetString("-log")))
                oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-log") + " for writing");
        }

        // Szybki options
        OESzybkiOptions opts = new OESzybkiOptions();

        // select optimization type
        if (itf.GetBool("-t")) {
            opts.SetRunType(OERunType.TorsionsOpt);
        }
        else {
            opts.SetRunType(OERunType.CartesiansOpt);
        }

        // select protein-electrostatic model
        String emodel = itf.GetString("-e");
        int elecModel = OEProteinElectrostatics.NoElectrostatics;
        if (emodel.equals("VdW"))
            elecModel = OEProteinElectrostatics.NoElectrostatics;
        else if (emodel.equals("PB"))
            elecModel = OEProteinElectrostatics.GridPB;
        else if (emodel.equals("Coulomb"))
            elecModel = OEProteinElectrostatics.GridCoulomb;
        else if (emodel.equals("ExactCoulomb"))
            elecModel = OEProteinElectrostatics.ExactCoulomb;
        opts.GetProteinOptions().SetProteinElectrostaticModel(elecModel);

        opts.GetProteinOptions().SetExactVdWProteinLigand(true);
        opts.GetOptOptions().SetMaxIter(1000);
        opts.GetOptOptions().SetGradTolerance(1e-6);

        // Szybki object
        OESzybki sz = new OESzybki(opts);

        // read and setup protein
        OEGraphMol protein = new OEGraphMol();
        oechem.OEReadMolecule(pfs, protein);
        sz.SetProtein(protein);

        // save or load grid potential
        if (emodel.equals("PB") || emodel.equals("Coulomb")) {
            if(itf.HasString("-s"))
                sz.SavePotentialGrid(itf.GetString("-s"));
            if(itf.HasString("-l"))
                sz.LoadPotentialGrid(itf.GetString("-l"));
        }

        // process molecules
        OEMol mol = new OEMol();
        while (oechem.OEReadMolecule(lfs, mol)) {
            logfile.write("\nMolecule " + mol.GetTitle() + "\n");
            boolean no_res = true;
            for (OESzybkiResults res : sz.call(mol)) {
                res.Print(logfile);
                no_res = false;
            }

            if (no_res) {
                oechem.OEThrow.Warning("No results processing molecule: " + mol.GetTitle());
                continue;
            }

            oechem.OEWriteMolecule(ofs, mol);
        }
    }
    private static String interfaceData =
"!PARAMETER -in\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Input molecule file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -p\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Input protein file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -out\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Output molecule file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -log\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF Log file name. Defaults to standard out.\n" +
"!END\n" +
"\n" +
"!PARAMETER -e\n" +
"!TYPE string\n" +
"!DEFAULT VdW\n" +
"!LEGAL_VALUE VdW\n" +
"!LEGAL_VALUE PB\n" +
"!LEGAL_VALUE Coulomb\n" +
"!LEGAL_VALUE ExactCoulomb\n" +
"!BRIEF Protein ligand electrostatic model.\n" +
"!END\n" +
"\n" +
"!PARAMETER -t\n" +
"!TYPE bool\n" +
"!DEFAULT false\n" +
"!REQUIRED false\n" +
"!BRIEF Torsions added to the optimized variables.\n" +
"!END\n" +
"\n" +
"!PARAMETER -l\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF File name of the potential grid to be read.\n" +
"!END\n" +
"\n" +
"!PARAMETER -s\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF File name of the potential grid to be saved.\n" +
"!END";
}

Optimization of a set of bound ligands in a partially flexible receptor

The next example is very similar with respect to the previous one, but the side chains of the protein which are within the specified range from the ligand, are made flexible during optimization (methods OESzybkiProteinOptions.SetProteinFlexibilityType and OESzybkiProteinOptions.SetProteinFlexibilityRange). Optionally, partially optimized structure can be saved to a file.

//*****************************************************************************
// Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FlexibleProtein {
    public static void main(String[] args) {
        OEInterface itf = new OEInterface ();
        oechem.OEConfigure(itf, interfaceData);
        oechem.OEParseCommandLine (itf, args, "FlexibleProtein");

        oemolistream lfs = new oemolistream();
        if (!lfs.open(itf.GetString("-in")))
          oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-in") + " for reading");

        oemolistream pfs = new oemolistream();
        if (!pfs.open(itf.GetString("-p")))
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-p") + " for reading");

        oemolostream olfs = new oemolostream();
        if (!olfs.open(itf.GetString("-out")))
            oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-out") + " for writing");

        oemolostream opfs = new oemolostream();
        if (itf.HasString("-s")) {
            if (!opfs.open(itf.GetString("-s")))
                oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-s") + " for writing");
        }

        oeostream logfile = oechem.oeout;
        if (itf.HasString("-log")) {
            if (!logfile.open(itf.GetString("-log")))
              oechem.OEThrow.Fatal("Unable to open " + itf.GetString("-log") + " for writing");
        }

        // Szybki options
        OESzybkiOptions opts = new OESzybkiOptions();

        // select optimization type
        String opt = itf.GetString("-opt");
        if (opt.equals("Cartesian"))
            opts.SetRunType(OERunType.CartesiansOpt);
        else if (opt.equals("Torsion"))
            opts.SetRunType(OERunType.TorsionsOpt);
        else if (opt.equals("SolidBody"))
            opts.SetRunType(OERunType.SolidBodyOpt);

        // select protein-electrostatic model
        String emodel = itf.GetString("-e");
        int elecModel = OEProteinElectrostatics.NoElectrostatics;
        if (emodel.equals("VdW"))
            elecModel = OEProteinElectrostatics.NoElectrostatics;
        else if (emodel.equals("PB"))
            elecModel = OEProteinElectrostatics.GridPB;
        else if (emodel.equals("Coulomb"))
            elecModel = OEProteinElectrostatics.GridCoulomb;
        else if (emodel.equals("ExactCoulomb"))
            elecModel = OEProteinElectrostatics.ExactCoulomb;
        opts.GetProteinOptions().SetProteinElectrostaticModel(elecModel);

        // use smooth potential and tight convergence
        if (emodel.equals("VdW") || emodel.equals("ExactCoulomb")) {
            opts.GetProteinOptions().SetExactVdWProteinLigand(true);
            opts.GetOptOptions().SetMaxIter(1000);
            opts.GetOptOptions().SetGradTolerance(1e-6);
        }

        // protein flexibility
        opts.GetProteinOptions().SetProteinFlexibilityType(OEProtFlex.SideChains);
        opts.GetProteinOptions().SetProteinFlexibilityRange(itf.GetDouble("-d"));

        // Szybki object
        OESzybki sz = new OESzybki(opts);

        // read and setup protein
        OEGraphMol protein = new OEGraphMol();
        OEGraphMol oprotein = new OEGraphMol();  // optimized protein
        oechem.OEReadMolecule(pfs, protein);
        sz.SetProtein(protein);
        pfs.close();

        // process molecules
        OEMol mol = new OEMol();
        while (oechem.OEReadMolecule(lfs, mol)) {
            logfile.write("\nMolecule " + mol.GetTitle() + "\n");
            for (OESzybkiResults res:  sz.call(mol))
                res.Print(logfile);

            oechem.OEWriteMolecule(olfs, mol);
            if (itf.HasString("-s")) {
                sz.GetProtein(oprotein);
                oechem.OEWriteMolecule(opfs, oprotein);
            }
        }
        opfs.close();
    }

    private static String interfaceData =
"!BRIEF -in input_molecule -p protein -out output_molecule\n" +
"!PARAMETER -in\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Input molecule file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -p\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Input protein file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -out\n" +
"!TYPE string\n" +
"!REQUIRED true\n" +
"!BRIEF Output molecule file name.\n" +
"!END\n" +
"\n" +
"!PARAMETER -log\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF Log file name. Defaults to standard out.\n" +
"!END\n" +
"\n" +
"!PARAMETER -e\n" +
"!TYPE string\n" +
"!DEFAULT VdW\n" +
"!LEGAL_VALUE VdW\n" +
"!LEGAL_VALUE PB\n" +
"!LEGAL_VALUE Coulomb\n" +
"!LEGAL_VALUE ExactCoulomb\n" +
"!BRIEF Protein ligand electrostatic model.\n" +
"!END\n" +
"\n" +
"!PARAMETER -opt\n" +
"!TYPE string\n" +
"!DEFAULT Cartesian\n" +
"!LEGAL_VALUE Cartesian\n" +
"!LEGAL_VALUE Torsion\n" +
"!LEGAL_VALUE SolidBody\n" +
"!BRIEF Optimization method\n" +
"!END\n" +
"\n" +
"!PARAMETER -d\n" +
"!TYPE double\n" +
"!DEFAULT 5.0\n" +
"!BRIEF Distance criteria from protein side-chains flexibility.\n" +
"!END\n" +
"\n" +
"!PARAMETER -s\n" +
"!TYPE string\n" +
"!REQUIRED false\n" +
"!BRIEF File name the partially optimized protein will be saved.\n"+
"!END";
}

Estimation of PB binding for a set of ligands

This example shows how to fast estimate binding for a set of ligands using PB electrostatics. Two OESzybki objects are instantiated: one for the optimization of bound ligands in VdW-Coulomb potential, and the second one which performs single-point PB calculations. Final results are attached as SD tags to the output molecules.

/**********************************************************************
 Copyright (C) 2009-2017 by OpenEye Scientific Software, Inc.
***********************************************************************/
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class ProteinLigandPB { 
    public static void main(String[] args)
    {
        if (args.length != 3) {
            oechem.OEThrow.Usage("ProteinLigandPB ligand_file protein_file output_file (SDF or OEB)");
        }

        String ligandFile = args[0];
        String proteinFile = args[1];
        String outFile = args[2];

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outFile + " for writing");
        }

        if (!oechem.OEIsSDDataFormat(ofs.GetFormat())) {
            oechem.OEThrow.Fatal("Output file does not support SD data used by this example");
        }

        // Szybki options for VdW-Coulomb calculations
        OESzybkiOptions optsC = new OESzybkiOptions();
        optsC.GetProteinOptions().SetProteinElectrostaticModel(OEProteinElectrostatics.ExactCoulomb);
        optsC.SetRunType(OERunType.CartesiansOpt);

        // Szybki options for PB calculations
        OESzybkiOptions optsPB = new OESzybkiOptions();
        optsPB.GetProteinOptions().SetProteinElectrostaticModel(OEProteinElectrostatics.SolventPBForces);
        optsPB.SetRunType(OERunType.SinglePoint);

        // Szybki objects
        OESzybki szC = new OESzybki(optsC);
        OESzybki szPB = new OESzybki(optsPB);

        // read and setup protein
        oemolistream pfs = new oemolistream();
        if (!pfs.open(proteinFile)) {
            oechem.OEThrow.Fatal("Unable to open " + proteinFile + " for reading");
        }
        OEGraphMol protein = new OEGraphMol();
        oechem.OEReadMolecule(pfs, protein);
        pfs.close();
        szC.SetProtein(protein);
        szPB.SetProtein(protein);

        // process molecules
        oemolistream lfs = new oemolistream();
        if (!lfs.open(ligandFile)) { 
            oechem.OEThrow.Fatal("Unable to open " + ligandFile + " for reading");
        }

        OEMol mol = new OEMol();
        while(oechem.OEReadMolecule(lfs, mol)) {
            OESzybkiResultsIter resultsiter = szC.call(mol); // optimize mol

            if (!resultsiter.IsValid()) {
                oechem.OEThrow.Warning("No results processing molecule: " + mol.GetTitle());
                continue;
            }

            // do single point with better electrostatics and output results
            resultsiter = szPB.call(mol);
            OEConfBaseIter confiter = mol.GetConfs();
            while(resultsiter.hasNext()) {
                OEConfBase conf = confiter.next();
                OESzybkiResults results =  resultsiter.next();
                for (int i = 0; i < OEPotentialTerms.Max; ++i) {
                    if (i == OEPotentialTerms.ProteinLigandInteraction || i == OEPotentialTerms.VdWProteinLigand || 
                        i == OEPotentialTerms.CoulombProteinLigand     || i == OEPotentialTerms.ProteinDesolvation ||
                        i == OEPotentialTerms.LigandDesolvation        || i == OEPotentialTerms.SolventScreening)
                    {
                        String energy = String.format("%9.4f", results.GetEnergyTerm(i));
                        oechem.OEAddSDData(conf, oeszybkilib.OEGetEnergyTermName(i), energy);
                    }
                }
            }
            oechem.OEWriteMolecule(ofs, mol);
        }
        ofs.close();
        lfs.close();
    }
}

Optimization of a bound ligand using Newton-Raphson method

The last example in this section illustrates how to use SzybkiTK to optimize a ligand in partially flexible protein with Newton-Raphson optimization method.

//*****************************************************************************
//* Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class LigandFlexProteinNewton
{
    public static void main(String argv[])
    {
        if (argv.length != 3)
            oechem.OEThrow.Usage("LigandFlexProteinNewton protein input_ligand output_ligand");

        String proteinFile = argv[0];
        String ligandFile  = argv[1];
        String outligandFile = argv[2];

        oemolistream pfs = new oemolistream();
        if (!pfs.open(proteinFile)) {
            oechem.OEThrow.Fatal("Unable to open " + proteinFile + " for reading");
        }

        oemolistream lfs = new oemolistream();
        if (!lfs.open(ligandFile)) {
            oechem.OEThrow.Fatal("Unable to open " + ligandFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outligandFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outligandFile + " for writing");
        }

        OEGraphMol mol = new OEGraphMol();
        OEGraphMol protein = new OEGraphMol();
        oechem.OEReadMolecule(lfs, mol);
        lfs.close();
        oechem.OEReadMolecule(pfs, protein);
        pfs.close();

        OESzybkiOptions opts = new OESzybkiOptions();
        opts.GetOptOptions().SetOptimizerType(OEOptType.NEWTON);
        opts.GetProteinOptions().SetProteinElectrostaticModel(OEProteinElectrostatics.ExactCoulomb);
        opts.GetProteinOptions().SetProteinFlexibilityType(OEProtFlex.Residues);
        opts.GetProteinOptions().SetProteinFlexibilityRange(2.0f);

        OESzybki sz = new OESzybki(opts);
        sz.SetProtein(protein);

        OESzybkiResults res = new OESzybkiResults();
        if (sz.call(mol, res)) {
            oechem.OEWriteMolecule(ofs, mol);
            res.Print(oechem.oeout);
        }
        ofs.close();
    }
}

Entropy estimation

Estimation of solution ligand entropy

The following 2 codes illustrate how to estimate compound entropy in solution with Szybki TK. In the second example partial charges used for Sheffield solvation are taken from the input molecule.

/**********************************************************************
 Copyright (C) 2014-2017 by OpenEye Scientific Software, Inc.
***********************************************************************/
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class LigandEntropy {
    public static void main(String[] args) {
        if (args.length != 1) { 
            oechem.OEThrow.Usage("LigandEntropy <lig>");
        }

        oemolistream ifs = new oemolistream();
        if (!ifs.open(args[0])) { 
            oechem.OEThrow.Fatal("Unable to open " + args[0] + " for reading");
        }

        OEMol lig = new OEMol();
        oechem.OEReadMolecule(ifs, lig);
        ifs.close();

        OESzybkiOptions opts = new OESzybkiOptions();
        OESzybki sz = new OESzybki(opts);
        OESzybkiEnsembleResults res = new OESzybkiEnsembleResults();
        double solutionEntropy = sz.GetEntropy(lig, res, OEEntropyMethod.Analytic);

        System.out.printf("Configurational entropy %10.2f\n", res.GetConfigurationalEntropy());
        System.out.printf("Solvation entropy       %10.2f\n", res.GetEnsembleLigSolvEntropy());
        System.out.printf("                            ======\n");
        System.out.printf("Total solution entropy  %10.2f J/(mol K)\n", solutionEntropy);
    }
}
//*****************************************************************************
//*  Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oequacpac.*;
import openeye.oeszybki.*;

public class SolutionEntropy
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("SolutionEntropy input_molecule output_molecule");
        }

        String inputFile = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        // input set of conformations : solution;
        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OESzybkiOptions opts = new OESzybkiOptions();
        opts.GetSolventOptions().SetChargeEngine(new OEChargeEngineNoOp());

        OESzybki sz = new OESzybki(opts);
        OESzybkiEnsembleResults eres = new OESzybkiEnsembleResults();
        double entropy = sz.GetEntropy(mol, eres, OEEntropyMethod.Analytic);

        
        DecimalFormat df = new DecimalFormat("####.#");
        oechem.OEThrow.Info("Estimated molar solution entropy of the input compound is: " + df.format(entropy) + " J/(mol*K)");
        oechem.OEThrow.Info("Vibrational entropies (in J/(mol*K)) for all conformations:");

        OEConfBaseIter conf = mol.GetConfs();
         for (OESzybkiResults res : eres.GetResultsForConformations()) {
             if (!conf.IsValid()) {
                continue;
             }

            System.out.printf("%2d %5.1f\n", conf.Target().GetIdx(), res.GetVibEntropy());
            conf.Increment();
        }

        // ensemble of unique conformations;
        oechem.OEWriteMolecule(ofs, mol);
        ofs.close();
    }
}

Estimation of bound ligand entropy

The next example below shows how to calculate entropy of a bound ligand.

//*****************************************************************************
//*  Copyright (C) 2014-2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class BoundEntropy
{
    public static void main(String argv[])
    {
        if (argv.length != 3) {
            oechem.OEThrow.Usage("BoundEntropy bound_ligand protein opt_ligand");
        }

        String ligandFile    = argv[0];
        String proteinFile   = argv[1];
        String optligandFile = argv[2];

        oemolistream lfs = new oemolistream();
        if (!lfs.open(ligandFile)) {
            oechem.OEThrow.Fatal("Unable to open " + ligandFile + " for reading");
        }

        oemolistream pfs = new oemolistream();
        if (!pfs.open(proteinFile)) {
            oechem.OEThrow.Fatal("Unable to open " + proteinFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(optligandFile)) {
            oechem.OEThrow.Fatal("Unable to open " + optligandFile + " for writing");
        }

        OEMol ligand = new OEMol();
        oechem.OEReadMolecule(lfs, ligand);
        lfs.close();

        OEGraphMol protein = new OEGraphMol();
        oechem.OEReadMolecule(pfs, protein);
        pfs.close();

        OESzybkiOptions opts = new OESzybkiOptions();
        OESzybki sz = new OESzybki(opts);
        sz.SetProtein(protein);

        OESzybkiEnsembleResults eres = new OESzybkiEnsembleResults();
        double entropy = sz.GetEntropy(ligand, eres, OEEntropyMethod.Analytic, OEEnvType.Protein);

        System.err.printf("Estimated entropy of the bound ligand is: %5.1f J/(mol*K)\n", entropy);

        oechem.OEWriteMolecule(ofs, ligand);
        ofs.close();
    }
}

Solvation free energy estimation

//*****************************************************************************
//* Copyright (C) 2014 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FreeFormSolv
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("FreeFormSolv <input> <output>");

        }
        String inputFile  = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OEFreeFormSolvOptions opts = new OEFreeFormSolvOptions();
        opts.SetIonicState(OEFreeFormIonicState.Uncharged);
        OEFreeFormSolvResults res = new OEFreeFormSolvResults();

        OEGraphMol omol = new OEGraphMol();
        if (!oeszybkilib.OEEstimateSolvFreeEnergy(res, omol, mol, opts)) {
            oechem.OEThrow.Error("Failed to calculate solvation free energy for molecule " + mol.GetTitle());
        }

        DecimalFormat df = new DecimalFormat("####.##");
        double solvenergy = res.GetSolvationFreeEnergy();
        oechem.OEThrow.Info("Solvation free energy for compound " + mol.GetTitle() + " is " + df.format(solvenergy) + " kcal/mol");

        oechem.OEWriteMolecule(ofs, omol);
        ofs.close();
    }
}

Conformations free energy estimation

Warning

This capability should not be used on 32-bit platforms because the memory requirements are too high.

Simple free energy estimation

The following code illustrates how to use the high level commands from OEFreeFormConf to estimate conformer free energies in solution with Szybki TK.

//*****************************************************************************
//* Copyright (C) 2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FreeFormConf
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("FreeFormConf <input> <output>");
        }

        String inputFile  = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OEFreeFormConfOptions opts = new OEFreeFormConfOptions();
        OEFreeFormConf ffconf = new OEFreeFormConf(opts);

        OEMol omol = new OEMol(mol);       
        if (!(ffconf.EstimateFreeEnergies(omol) == OEFreeFormReturnCode.Success)) {
            oechem.OEThrow.Error("Failed to estimate conformational free energies for molecule " + mol.GetTitle());
        }

        OEFreeFormConfResults res = new OEFreeFormConfResults(omol);
        oechem.OEThrow.Info("Number of unique conformations: " + res.GetNumUniqueConfs());
        oechem.OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
        oechem.OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
        for (OESingleConfResult r : res.GetResultsForConformations()) {
            System.err.printf("%2d %10.2f %14.2f\n", r.GetConfIdx(), r.GetDeltaG(), r.GetVibrationalEntropy());
        }
        oechem.OEWriteMolecule(ofs, omol);
        ofs.close();
    }
}

Simple restriction energy estimation

The following code illustrates how to use the high level commands from OEFreeFormConf to estimate restriction energies on conformers, along with conformer free energies in solution, with Szybki TK.

//*****************************************************************************
//* Copyright (C) 2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FreeFormConfRstr
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("FreeFormConfRstr <input> <output>");
        }

        String inputFile  = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OEFreeFormConfOptions opts = new OEFreeFormConfOptions();
        OEFreeFormConf ffconf = new OEFreeFormConf(opts);

        OEMol omol = new OEMol(mol);
        OEMol rmol = new OEMol(mol);
        if (ffconf.EstimateFreeEnergies(omol, rmol) != OEFreeFormReturnCode.Success) {
            oechem.OEThrow.Error("Failed to estimate conformational free energies for molecule " + mol.GetTitle());
        }

        OEFreeFormConfResults res = new OEFreeFormConfResults(omol);
        oechem.OEThrow.Info("Number of unique conformations: " + res.GetNumUniqueConfs());
        oechem.OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
        oechem.OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
        for (OESingleConfResult r : res.GetResultsForConformations()) {
            System.err.printf("%2d %10.2f %14.2f\n", r.GetConfIdx(), r.GetDeltaG(), r.GetVibrationalEntropy());
        }

        OERestrictionEnergyResult rstrRes = new OERestrictionEnergyResult(rmol);
        oechem.OEThrow.Info("Global strain: " + rstrRes.GetGlobalStrain());
        oechem.OEThrow.Info("Local strain: " + rstrRes.GetLocalStrain());

        oechem.OEWriteMolecule(ofs, omol);
        ofs.close();
    }
}

Advanced free energy estimation

The following code illustrates how to use the low level commands from OEFreeFormConfAdvanced to estimate the conformer free energies in solution with Szybki TK. These low level methods of estimation gives an advantage over the high level methods of OEFreeFormConf in that these gives the user control over better managing certain expensive parts of the calculation.

//*****************************************************************************
//* Copyright (C) 2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FreeFormConfAdvanced
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("FreeFormConfAdvanced <input> <output>");
        }

        String inputFile  = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OEFreeFormConfOptions opts = new OEFreeFormConfOptions();
        OEFreeFormConfAdvanced ffconf = new OEFreeFormConfAdvanced(opts);

        // Make a copy of our MCMol.  We will execute the FreeFormConf commands on
        // the copied molecule so that our original molecule stays intact. 
        OEMol omol = new OEMol(mol);

        // Prepare a comprehensive ensemble of molecule conformers. This will
		// generate a comprehensive set of conformers, assign solvent charges on the molecule
		// and check that the ensemble is otherwise ready for FreeFormConf calculations.
        if(ffconf.PrepareEnsemble(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Failed to prepare ensemble for FreeFormConf calculations");
        }

        // Perform loose optimization of the ensemble conformers.  We will remove
        // duplicates based on the loose optimization, to reduce the time needed for
        // tighter, more stricter optimization
        if(ffconf.PreOptimizeEnsemble(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Pre-optimization of the ensembles failed");
        }

        // Remove duplicates from the pre-optimized ensemble
        if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Duplicate removal from the ensembles failed");
        }

        // Perform the desired optimization.  This uses a stricter convergence
        // criteria in the default settings.
        if(ffconf.Optimize(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Optimization of the ensembles failed");
        }

        // Remove duplicates to obtain the set of minimum energy conformers
        if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Duplicate removal from the ensembles failed");
        }

        // Perform FreeFormConf free energy calculations.  When all the above steps
        // have already been performed on the ensemble, this energy calculation
        // step is fast.
        if(ffconf.EstimateEnergies(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Estimation of FreeFormConf energies failed");
        }

        // Gather results of calculation into a results object for ease of viewing, etc.
        OEFreeFormConfResults res = new OEFreeFormConfResults(omol);
        oechem.OEThrow.Info("Number of unique conformations: " + res.GetNumUniqueConfs());
        oechem.OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
        oechem.OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
        for (OESingleConfResult r : res.GetResultsForConformations()) {
            System.err.printf("%2d %10.2f %14.2f\n", r.GetConfIdx(), r.GetDeltaG(), r.GetVibrationalEntropy());
        }
        
        oechem.OEWriteMolecule(ofs, omol);
        ofs.close();
    }
}

Advanced restriction energy estimation

The following code illustrates how to use the low level commands from OEFreeFormConfAdvanced to estimate the restriction energies of conformers, along with conformer free energies in solution, with Szybki TK. These low level methods of estimation gives an advantage over the high level methods of OEFreeFormConf in that these gives the user control over better managing certain expensive parts of the calculation.

//*****************************************************************************
//* Copyright (C) 2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FreeFormConfAdvRstr
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("FreeFormConfAdvRstr <input> <output>");
        }

        String inputFile  = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OEFreeFormConfOptions opts = new OEFreeFormConfOptions();
        OEFreeFormConfAdvanced ffconf = new OEFreeFormConfAdvanced(opts);

        // Make a copy of our MCMol.  We will execute the FreeFormConf commands on
        // the copied molecule so that our original molecule stays intact. 
        OEMol omol = new OEMol(mol);

		// Make further copies of our original molecule.  The copied molecule(s) would be used
		// as source on which retriction energies would be calculated
		OEMol rmol = new OEMol(mol);
		OEMol fmol = new OEMol(mol);

        // Prepare a comprehensive ensemble of molecule conformers.  For calculation
		// of restriction energies we want to make sure that all the corresponding free
		// conformers are also part of the comprehensive ensemble.  This will also
		// assign solvent charges on the molecule and check that the ensemble is
		// otherwise ready for FreeFormConf calculations. The resulting `fmol`
		// contains the correspondig free conformers.
        if(ffconf.PrepareEnsemble(omol, rmol, fmol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Failed to prepare ensemble for FreeFormConf calculations");
        }

        // Perform loose optimization of the ensemble conformers.  We will remove
        // duplicates based on the loose optimization, to reduce the time needed for
        // tighter, more stricter optimization
        if(ffconf.PreOptimizeEnsemble(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Pre-optimization of the ensembles failed");
        }

        // Remove duplicates from the pre-optimized ensemble
        if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Duplicate removal from the ensembles failed");
        }

        // Perform the desired optimization.  This uses a stricter convergence
        // criteria in the default settings.
        if(ffconf.Optimize(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Optimization of the ensembles failed");
        }

        // Remove duplicates to obtain the set of minimum energy conformers
        if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Duplicate removal from the ensembles failed");
        }

        // Perform FreeFormConf free energy calculations.  When all the above steps
        // have already been performed on the ensemble, this energy calculation
        // step is fast.
        if(ffconf.EstimateEnergies(omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Estimation of FreeFormConf energies failed");
        }

        // Gather results of calculation into a results object for ease of viewing, etc.
        OEFreeFormConfResults res = new OEFreeFormConfResults(omol);
        oechem.OEThrow.Info("Number of unique conformations: " + res.GetNumUniqueConfs());
        oechem.OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
        oechem.OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
        for (OESingleConfResult r : res.GetResultsForConformations()) {
            System.err.printf("%2d %10.2f %14.2f\n", r.GetConfIdx(), r.GetDeltaG(), r.GetVibrationalEntropy());
        }

		// Identify the corresponding conformer(s) to the free minimized conformer(s).
        // If identified, the corresponding (Conf)Free energy information is also
        // copied to the free conformers
        if(ffconf.IdentifyConformer(fmol, omol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Identifican of free conformer(s) failed");
        }

        // Estimate restriction energies. Since both restricted and free conformer
        // energy components are already available, this operation is fast.
        if(ffconf.EstimateRestrictionEnergy(fmol, rmol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Restriction energy estimation failed");
        }

        // Gather restriction energies into a results object for ease of viewing, etc.
        OERestrictionEnergyResult rstrRes = new OERestrictionEnergyResult(fmol);
        oechem.OEThrow.Info("Global strain: " + rstrRes.GetGlobalStrain());
        oechem.OEThrow.Info("Local strain: " + rstrRes.GetLocalStrain());

        // Optionally it is desired to perform a restrained optimization of the
        // restricted conformer(s) to brush out any energy differences due to
        // force field constaints or the sources of coonformer coordinates.  Note: The
        // high level EstimateFreeEnergy method does not perform this opertion.
        if(ffconf.OptimizeRestraint(rmol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Restraint optimization of the conformer(s) failed");
        }

        // Estimate restriction energies. Since both restricted and free conformer
        // energy components are already available, this operation is fast.
        if(ffconf.EstimateRestrictionEnergy(fmol, rmol)!=OEFreeFormReturnCode.Success){
            oechem.OEThrow.Error("Restriction energy estimation failed");
        }

        // Gather restriction energies into a results object for ease of viewing, etc.
        OERestrictionEnergyResult rstrRes2 = new OERestrictionEnergyResult(fmol);
        oechem.OEThrow.Info("Global strain: " + rstrRes2.GetGlobalStrain());
        oechem.OEThrow.Info("Local strain: " + rstrRes2.GetLocalStrain());
        
        oechem.OEWriteMolecule(ofs, omol);
        ofs.close();
    }
}

Finding similar conformers

The following code illustrates how to find similar conformers to the ones at hand, from a pool of minimum energy conformers.

//*****************************************************************************
//* Copyright (C) 2017 OpenEye Scientific Software, Inc.
//*****************************************************************************
package openeye.examples.oeszybki;

import java.text.DecimalFormat;

import openeye.oechem.*;
import openeye.oeszybki.*;

public class FindSimilarConfs
{
    public static void main(String argv[])
    {
        if (argv.length != 2) {
            oechem.OEThrow.Usage("FindSimilarConfs <input> <output>");
        }

        String inputFile  = argv[0];
        String outputFile = argv[1];

        oemolistream ifs = new oemolistream();
        if (!ifs.open(inputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + inputFile + " for reading");
        }

        oemolostream ofs = new oemolostream();
        if (!ofs.open(outputFile)) {
            oechem.OEThrow.Fatal("Unable to open " + outputFile + " for writing");
        }

        OEMol mol = new OEMol();
        oechem.OEReadMolecule(ifs, mol);
        ifs.close();

        OEFreeFormConfOptions opts = new OEFreeFormConfOptions();
        OEFreeFormConf ffconf = new OEFreeFormConf(opts);

        // Estimate free energies to ontain the minimum energy conformers
        OEMol omol = new OEMol(mol);
        if (!(ffconf.EstimateFreeEnergies(omol) == OEFreeFormReturnCode.Success)) {
            oechem.OEThrow.Error("Failed to estimate conformational free energies for molecule " + mol.GetTitle());
        }

        // Find similar conformers to the ones we started with, from the
        // pool of minimum energy conformers
        OEMol fmol = new OEMol(mol);
        for (OEConfBase conf : mol.GetConfs())
            ffconf.FindSimilarConfs(fmol, omol, conf, new OESimilarByRMSD(0.05));

        oechem.OEWriteMolecule(ofs, fmol);
        ofs.close();
    }
}