Examples: 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 compound properties, such as solvation free energy or free energy of selecting specific conformations out of the ensemble, the higher level API are 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. Optimization is done using the default MMFF94. Final energy results are available as SD tags in the returned molecule and optionally with a call of a method OESzybkiResults::Print.

Listing 1: Simple Ligand in a Vacuum

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

int main(int argc, char* argv[])
{
  if (argc!=3)
    OEThrow.Usage("%s <input> <output>", argv[0]);

  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

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

  OEGraphMol mol;
  OEReadMolecule(ifs, mol);

  OESzybkiOptions opts;
  OESzybki sz(opts);
  OESzybkiResults res;
  if (!sz(mol, res)) return 1;

  OEWriteMolecule(ofs, mol);
  res.Print(OEOUT);
  return 0;
}

Download code

simple_vacuum.cpp

Single ligand in vacuum using SMIRNOFF

The following example illustrates how to optimize a single ligand in vacuum using SMIRNOFF (see SMIRNOFF). 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(). Partial charges have to be preassigned to the input molecule. 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.

Listing 2: Simple Ligand in a Vacuum Using SMIRNOFF

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

int main(int argc, char* argv[])
{
  if (argc!=3)
    OEThrow.Usage("%s <input> <output>", argv[0]);

  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

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

  OEGraphMol mol;
  OEReadMolecule(ifs, mol);

  OESzybkiOptions opts;

  //selection of SMIRNOFF force field 
  //Note: Partial charges have to be preassigned to the input molecule for example in the mol2 file
  opts.GetGeneralOptions().SetForceFieldType(OEForceFieldType::SMIRNOFF99FROSST);

  OESzybki sz(opts);
  OESzybkiResults res;
  if (!sz(mol, res)) return 1;

  OEWriteMolecule(ofs, mol);
  res.Print(OEOUT);
  return 0;
}

Download code

simple_vacuum.cpp

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.

Listing 3: Optimization of a Set of Ligands

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oeszybki.h"
#include "ligand_multiple_options.itf"

using namespace OESystem;
using namespace OEPlatform;
using namespace OEChem;
using namespace OESz;
using namespace std;

int main(int argc, char* argv[])
{
  //read and parse command line
  OEInterface itf(InterfaceData, argc, argv);

  //opening files
  oemolistream ifs;
  if (!ifs.open(itf.Get<string>("-in")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-in").c_str());

  oemolostream ofs;
  if (!ofs.open(itf.Get<string>("-out")))
    OEThrow.Fatal("Unable to open %s for writing", itf.Get<string>("-out").c_str());

  oeofstream logfile = oeout;
  if (itf.Has<string>("-log"))
    if (!logfile.open(itf.Get<string>("-log")))
      OEThrow.Fatal("Unable to open %s for logging",
                    itf.Get<string>("-log").c_str());

  //Szybki options
  OESzybkiOptions opts;

  //select run type
  if (itf.Get<bool>("-t"))
    opts.SetRunType(OERunType::TorsionsOpt);
  if (itf.Get<bool>("-n"))
    opts.SetRunType(OERunType::SinglePoint);

  //apply solvent model
  if (itf.Get<bool>("-s"))
    opts.GetSolventOptions().SetSolventModel(OESolventModel::Sheffield);

  //remove attractive VdW forces
  if (itf.Get<bool>("-a"))
	  opts.GetGeneralOptions().SetRemoveAttractiveVdWForces(true);

  //Szybki object
  OESzybki sz(opts);

  //fix atoms
  if (itf.Has<string>("-f"))
    if (!sz.FixAtoms(itf.Get<string>("-f")))
      OEThrow.Warning("Failed to fix atoms for pattern %s",
                      itf.Get<string>("-f").c_str());

  //process molecules
  OEMol mol;
  OEIter<OESzybkiResults> results;   //iterator to energy results
  while(OEReadMolecule(ifs, mol))
  {
    logfile <<"\nMolecule "<<mol.GetTitle() << "\n";
    bool no_res = true;
    for (results = sz(mol); results; ++results)
    {
      results->Print(logfile);
      no_res = false;
    }

    if (no_res)
    {
      OEThrow.Warning("No results processing molecule: %s", mol.GetTitle());
      continue;
    }

    OEWriteMolecule(ofs, mol);
  }
  return 0;
}

Download code

ligand_multiple_options.cpp

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:

Listing 4: Optimization of a Single Ligand with the Newton-Raphson Method

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include <oechem.h>
#include <oeszybki.h>

using namespace OEPlatform;
using namespace OEChem;
using namespace OESystem;
using namespace OESz;

int main(int argc,char** argv)
{
  if (argc != 3)
    OEThrow.Usage("%s input_molecule output_molecule", argv[0]);

  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

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

  OEGraphMol mol;
  OEReadMolecule(ifs, mol);

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

  OESzybki sz(opts);
  OESzybkiResults res;
  if (sz(mol, res))
  {
    OEWriteMolecule(ofs, mol);
    res.Print(oeout);
  }

  return 0;
}

Download code

simple_newton.cpp

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.

Listing 5: Optimization of All Conformers of a Ligand

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include <openeye.h>
#include <oechem.h>
#include <oeszybki.h>

using namespace OEPlatform;
using namespace OEChem;
using namespace OESystem;
using namespace OESz;

int main(int argc, char** argv)
{
  if (argc != 3)
    OEThrow.Usage("%s input_molecule output_molecule", argv[0]);

  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

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

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

  OESzybki sz(opts);
  OEMol mol;
  while (OEReadMolecule(ifs, mol))
  {
    for (OEIter<OEConfBase> conf = mol.GetConfs(); conf; ++conf)
    {
      OESzybkiResults res;
      if (sz(conf, res))
      {
        std::string energy = OENumberToString(res.GetTotalEnergy());
        OESetSDData(conf, "Total_energy", energy);
      }
    }
    OEWriteMolecule(ofs, mol);
  }

  return 0;
}

Download code

simple_conformer.cpp

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.

Listing 5: Optimization of a Single Bound Ligand

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

int main(int argc, char* argv[])
{
  if (argc!=4)
    OEThrow.Usage("%s <protein> <ligand> <output>", argv[0]);

  oemolistream lfs;
  if (!lfs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

  oemolistream pfs;
  if (!pfs.open(argv[2]))
    OEThrow.Fatal("Unable to open %s for reading", argv[2]);

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

  OEGraphMol mol;
  OEReadMolecule(lfs, mol);

  OEGraphMol protein;
  OEReadMolecule(pfs, protein);

  OESzybkiOptions opts;
  OESzybki sz(opts);
  sz.SetProtein(protein);
  OESzybkiResults res;
  if (!sz(mol, res)) return 1;

  OEWriteMolecule(ofs, mol);
  res.Print(OEOUT);

  return 0;
}

Download code

simple_protein.cpp

Protein-Ligand Energetics and Optimization

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

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.

Listing 6: Optimization of a Set of Bound Ligands in a Rigid Receptor

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oeszybki.h"
#include "protein_multiple_options.itf"

using namespace OESystem;
using namespace OEPlatform;
using namespace OESz;
using namespace OEChem;
using namespace std;

int main(int argc, char* argv[])
{
  //read and parse command line
  OEInterface itf(InterfaceData, argc, argv);

  //opening files
  oemolistream lfs;
  if (!lfs.open(itf.Get<string>("-in")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-in").c_str());

  oemolistream pfs;
  if (!pfs.open(itf.Get<string>("-p")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-p").c_str());

  oemolostream ofs;
  if (!ofs.open(itf.Get<string>("-out")))
    OEThrow.Fatal("Unable to open %s for writing", itf.Get<string>("-out").c_str());

  oeofstream logfile = oeout;
  if (itf.Has<string>("-log"))
    if (!logfile.open(itf.Get<string>("-log")))
      OEThrow.Fatal("Unable to open %s for logging",
                    itf.Get<string>("-log").c_str());

  //Szybki options object
  OESzybkiOptions opts;

  //select optimization type
  if (itf.Get<bool>("-t")) opts.SetRunType(OERunType::TorsionsOpt);
  else opts.SetRunType(OERunType::CartesiansOpt);

  //select protein-electrostatic model
  string emodel = itf.Get<string>("-e");
  unsigned int elecModel = OEProteinElectrostatics::NoElectrostatics;
  if (emodel == "VdW")
	  elecModel = OEProteinElectrostatics::NoElectrostatics;
  else if (emodel == "PB")
	  elecModel = OEProteinElectrostatics::GridPB;
  else if (emodel == "Coulomb")
	  elecModel = OEProteinElectrostatics::GridCoulomb;
  else if (emodel == "ExactCoulomb")
  {
	  elecModel = OEProteinElectrostatics::ExactCoulomb;
	  opts.GetProteinOptions().SetExactVdWProteinLigand(true);
	  opts.GetOptOptions().SetMaxIter(1000);
	  opts.GetOptOptions().SetGradTolerance(1e-6);
  }  
  opts.GetProteinOptions().SetProteinElectrostaticModel(elecModel);

  //Szybki object
  OESzybki sz(opts);

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

  //save or load grid potential
  if (emodel == "PB" || emodel == "Coulomb")
  {
    if (itf.Has<string>("-s")) sz.SavePotentialGrid(itf.Get<string>("-s"));
    if (itf.Has<string>("-l")) sz.LoadPotentialGrid(itf.Get<string>("-l"));
  }

  //process molecules
  OEMol mol;
  OEIter<OESzybkiResults> res;   //iterator to energy results
  while(OEReadMolecule(lfs, mol))
  {
    logfile << "\nMolecule " << mol.GetTitle() << "\n";
    bool no_res = true;
    for (res = sz(mol); res; ++res)
    {
      res->Print(logfile);
      no_res = false;
    }

    if (no_res)
    {
      OEThrow.Warning("No results processing molecule: %s", mol.GetTitle());
      continue;
    }

    OEWriteMolecule(ofs, mol);

  }
  return 0;
}

Download code

protein_multiple_options.cpp

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.

Listing 7: Optimization of a Set of Bound Ligands in a Partially Flexible Receptor

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oeszybki.h"
#include "flexible_protein.itf"

using namespace OESystem;
using namespace OEPlatform;
using namespace OESz;
using namespace OEChem;
using namespace std;

int main(int argc, char* argv[])
{
  //read and parse command line
  OEInterface itf(InterfaceData, argc, argv);

  //opening files
  oemolistream lfs;
  if (!lfs.open(itf.Get<string>("-in")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-in").c_str());

  oemolistream pfs;
  if (!pfs.open(itf.Get<string>("-p")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-p").c_str());

  oemolostream olfs;
  if (!olfs.open(itf.Get<string>("-out")))
    OEThrow.Fatal("Unable to open %s for writing", itf.Get<string>("-out").c_str());

  oemolostream opfs;
  if (itf.Has<string>("-s"))
    if (!opfs.open(itf.Get<string>("-s")))
      OEThrow.Error("Unable to open %s for protein output",
                    itf.Get<string>("-s").c_str());

  oeofstream logfile = oeout;
  if (itf.Has<string>("-log"))
    if (!logfile.open(itf.Get<string>("-log")))
      OEThrow.Error("Unable to open %s for logging",
                    itf.Get<string>("-log").c_str());

  //Szybki options
  OESzybkiOptions opts;

  //select optimization type
  string opt = itf.Get<string>("-opt");
  if (opt == "Cartesian")
    opts.SetRunType(OERunType::CartesiansOpt);
  if (opt == "Torsion")
    opts.SetRunType(OERunType::TorsionsOpt);
  if (opt == "SolidBody")
    opts.SetRunType(OERunType::SolidBodyOpt);

  //select protein-electrostatic model
  string emodel = itf.Get<string>("-e");
  unsigned int elecModel = OEProteinElectrostatics::NoElectrostatics;
  if (emodel == "VdW")
    elecModel = OEProteinElectrostatics::NoElectrostatics;
  else if (emodel == "PB")
	  elecModel = OEProteinElectrostatics::GridPB;
  else if (emodel == "Coulomb")
	  elecModel = OEProteinElectrostatics::GridCoulomb;
  else if (emodel == "ExactCoulomb")
	  elecModel = OEProteinElectrostatics::ExactCoulomb;
  opts.GetProteinOptions().SetProteinElectrostaticModel(elecModel);

  //use smooth potential and tight convergence
  if (emodel == "VdW" || emodel == "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.Get<double>("-d"));

  //Szybki object
  OESzybki sz(opts);

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

  //process molecules
  OEMol mol;
  OEIter<OESzybkiResults> results;   //iterator to energy results
  while (OEReadMolecule(lfs, mol))
  {
    logfile << "\nMolecule " << mol.GetTitle() << "\n";
    for (results = sz(mol); results; ++results)
      results->Print(logfile);

    OEWriteMolecule(olfs, mol);
    if (itf.Has<string>("-s"))
    {
      sz.GetProtein(oprotein);
      OEWriteMolecule(opfs, oprotein);
    }
  }
  return 0;
}

Download code

flexible_protein.cpp

Optimization of a bound ligand in a partially flexible receptor

The next example is very similar with respect to the previous one, but the list of residue numbers of the protein that are made flexible during optimization is specified by the -residues flag (see method OESzybkiProteinOptions::AddFlexibleResidue). Partially optimized ligand and protein structures are saved to files specified by the -outl and -outp flags, respectively.

Listing 8: Optimization of a Bound Ligand in a Partially Flexible Feceptor

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include <openeye.h>
#include <vector>
#include <oesystem.h>
#include <oechem.h>
#include <oeszybki.h>

#include "flex_residues.itf"

using namespace std;
using namespace OESystem;
using namespace OEChem;
using namespace OESz;


int main(int argc, char* argv[])
{
  OEInterface itf(InterfaceData, argc, argv);

  oemolistream ifs;
  if (!ifs.open(itf.Get<string>("-in")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-in").c_str());
  oemolistream pfs;
  if (!pfs.open(itf.Get<string>("-protein")))
    OEThrow.Fatal("Unable to open %s for reading", itf.Get<string>("-protein").c_str());
  oemolostream ofs;
  if (!ofs.open(itf.Get<string>("-outl")))
    OEThrow.Fatal("Unable to open %s for writing", itf.Get<string>("-outl").c_str());
  oemolostream opfs;
  if (!opfs.open(itf.Get<string>("-outp")))
    OEThrow.Fatal("Unable to open %s for writing", itf.Get<string>("-outp").c_str());

  OEGraphMol ligand;
  OEReadMolecule(ifs, ligand);
  OEGraphMol protein;
  OEReadMolecule(pfs, protein);

  //flexible residue numbers
  vector<int> res_num;
  for(OEIter<const int> i = itf.GetList<int>("-residues"); i; ++i)
    res_num.push_back(*i);

  OESzybkiOptions opts;
  opts.GetGeneralOptions().SetForceFieldType(OEForceFieldType::MMFF94S);
  opts.GetProteinOptions().SetProteinElectrostaticModel(OEProteinElectrostatics::ExactCoulomb);
  opts.GetProteinOptions().SetProteinFlexibilityType(OEProtFlex::SideChainsList);
  opts.SetRunType(OERunType::CartesiansOpt);
  opts.GetOptOptions().SetMaxIter(2000u);
  opts.GetOptOptions().SetGradTolerance(1.0e-6);

  for(size_t k = 0; k < res_num.size(); ++k)
  {
    for(OEIter<const OEAtomBase> atom = protein.GetAtoms(); atom; ++atom)
    {
      const OEResidue& res = OEAtomGetResidue(atom);
      if(res.GetResidueNumber() == res_num[k])
      {
        opts.AddFlexibleResidue(res);
        break;
      }
    }
  }

  //optimization
  OESzybki sz(opts);
  sz.SetProtein(protein);
  OESzybkiResults res;
  sz(ligand, res);

  OEWriteMolecule(ofs, ligand);
  sz.GetProtein(protein);
  OEWriteMolecule(opfs, protein);

  return 0;
}

Download code

flex_residues.cpp

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.

Listing 9: Estimation of PB Binding for a Set of Ligands

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEPlatform;
using namespace OESz;
using namespace OEChem;
using namespace std;

int main(int argc, char* argv[])
{
  if (argc != 4) OEThrow.Usage("%s ligand_file protein_file output_file (SDF or OEB)", argv[0]);

  oemolistream lfs;
  if (!lfs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

  oemolistream pfs;
  if (!pfs.open(argv[2]))
    OEThrow.Fatal("Unable to open %s for reading", argv[2]);

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

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

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

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

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

  //read and setup protein
  OEGraphMol protein;
  OEReadMolecule(pfs, protein);
  szC.SetProtein(protein);
  szPB.SetProtein(protein);

  //process molecules
  OEMol mol;
  OEIter<OESzybkiResults> results;   //iterators to energy results
  OEIter<OEConfBase>conf;            //iterator over conformations
  char buf[20];
  while(OEReadMolecule(lfs, mol))
  {
    //optimize mol
    bool no_res = true;
    for (results=szC(mol); results; ++results)
    {
      no_res = false;
      break;
    }

    if (no_res)
    {
      OEThrow.Warning("No results processing molecule: %s", mol.GetTitle());
      continue;
    }

    //do single point with better electrostatics and output results
    for (results = szPB(mol), conf = mol.GetConfs(); (bool)results && (bool)conf; ++results, ++conf)
    {
      for (unsigned 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)
        {
          snprintf(buf, 20, "%9.4f", results->GetEnergyTerm(i));
          OEAddSDData(conf, OEGetEnergyTermName(i), string(buf));
        }
      }
    }
    OEWriteMolecule(ofs, mol);
  }
  return 0;
}

Download code

protein_ligand_PB.cpp

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.

Listing 10: Optimization of a Bound Ligand Using Newton-Raphson Method

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OEPlatform;
using namespace OEChem;
using namespace OESystem;
using namespace OESz;

int main(int argc,char** argv)
{
  if (argc != 4)
    OEThrow.Usage("%s protein input_ligand output_ligand", argv[0]);

  oemolistream pfs;
  if (!pfs.open(argv[1]))
    OEThrow.Fatal("Unable to open %s for reading", argv[1]);

  oemolistream lfs;
  if (!lfs.open(argv[2]))
    OEThrow.Fatal("Unable to open %s for reading", argv[2]);

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

  OEGraphMol mol;
  OEGraphMol protein;
  OEReadMolecule(lfs, mol);
  OEReadMolecule(pfs, protein);

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

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

  OESzybkiResults res;
  if (sz(mol, res))
  {
    OEWriteMolecule(ofs, mol);
    res.Print(oeout);
  }

  return 0;
}

Download code

ligand_flexprotein_newton.cpp

DU Protein-Ligand Optimization with FF14SB-Parsley

Examples in this section show how to optimize a ligand in the protein active site using the OEFixedProteinLigandOptimizer and the OEFlexProteinLigandOptimizer classes. Both of these classes allow use of the FF14SB-Parsley forcefield, along with the MMFF and the MMFF-AMBER forcefields. These examples also demonstrate how to use design units as a source for the protein or the ligand. The first example optimizes an external ligand in a protein active site from a design unit, and the second example uses both the ligand and the protein from the same design unit.

Optimization of ligand in a rigid active site

This examples shows how to optimize a series of ligands in a rigid active site. The active site input is taken from a OEDesignUnit.

Optimization of ligand in a partially flexible active site

This examples shows how to optimize a ligand in a partially flexible active site. In this example, both the ligand and the active site is taken from a single OEDesignUnit.

Entropy estimation

Estimation of solution ligand entropy

The following code illustrates how to estimate compound entropy in solution with Szybki TK.

Listing 13: Ligand Entropy

/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OEPlatform;
using namespace OESystem;
using namespace OEChem;
using namespace OESz;

int main(int argc, char* argv[])
{
  OELigandEntropyOptions opts;
  OESimpleAppOptions options(opts,"LigandEntropy",OEFileStringType::Mol3D);
  if(OEConfigureOpts(options,argc,argv,false) == OEOptsConfigureStatus::Help)
    return 0;
  opts.UpdateValues(options);

  oemolistream ifs;
  if(!ifs.open(options.GetInFile()))
    OEThrow.Fatal("Unable to open %s for reading ", options.GetInFile().c_str());

  OEMol ensemble;
  while(OEReadMolecule(ifs,ensemble))
  {
    OEThrow.Info("Title: %s", ensemble.GetTitle());
    OEEntropyResults res;
    if(OELigandEntropy(res,ensemble,opts) == OESzybkiReturnCode::Success)
    {
      OEThrow.Info("Configurational entropy: %0.2f J/(mol K)",res.GetTotalEntropy());
      OEThrow.Info("Translational entropy  : %0.2f J/(mol K)",res.GetTranslationalEntropy());
      OEThrow.Info("VibRotEntropy entropy  : %0.2f J/(mol K)",res.GetVibrationalRotationalEntropy());
      OEThrow.Info("Solvation entropy      : %0.2f J/(mol K)",res.GetSolvationEntropy());
      OEThrow.Info("Total entropy          : %0.2f J/(mol K)",res.GetTotalEntropy());
    }
    else
    {
      OEThrow.Warning("Failed to evaluate entropy");
    }
  }
  return 0;
}

Download code

ligand_entropy.cpp

Estimation of bound ligand entropy

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

Listing 15: Estimation of Bound Ligand Entropy

/*
(C) 2022 Cadence Design Systems, Inc. (Cadence)
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oesystem.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OEPlatform;
using namespace OESystem;
using namespace OEChem;
using namespace OESz;

int main(int argc, char* argv[])
{
  OEBoundEntropyOptions opts;
  OERefInputAppOptions options(opts,"BoundEntropy",OEFileStringType::Mol3D,
                               OEFileStringType::DU,"-du");
  if(OEConfigureOpts(options,argc,argv,false) == OEOptsConfigureStatus::Help)
    return 0;
  opts.UpdateValues(options);

  oemolistream ifs;
  if(!ifs.open(options.GetInFile()))
    OEThrow.Fatal("Unable to open %s for reading ", options.GetInFile().c_str());

  oeifstream rfs;
  if(!rfs.open(options.GetRefFile()))
    OEThrow.Fatal("Unable to open %s for reading ", options.GetRefFile().c_str());

  OEBio::OEDesignUnit du;
  if(!OEReadDesignUnit(rfs,du))
    OEThrow.Fatal("Failed to read design unit");

  OEMol poses;
  while(OEReadMolecule(ifs,poses))
  {
    OEThrow.Info("Title: %s", poses.GetTitle());
    OEEntropyResults res;
    if(OEBoundLigandEntropy(res,du,OEBio::OEDesignUnitComponents::Protein,poses,opts) ==
       OESzybkiReturnCode::Success)
    {
      OEThrow.Info("Configurational entropy: %0.2f J/(mol K)",res.GetTotalEntropy());
      OEThrow.Info("Solvation entropy      : %0.2f J/(mol K)",res.GetSolvationEntropy());
      OEThrow.Info("Total entropy          : %0.2f J/(mol K)",res.GetTotalEntropy());
    }
    else
    {
      OEThrow.Warning("Failed to evaluate entropy");
    }
  }
  return 0;
}

Download code

bound_entropy.cpp

Solvation free energy estimation

Listing 16: Solvation Free Energy Estimation

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

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

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

  OEMol mol;
  OEReadMolecule(ifs, mol);

  OEFreeFormSolvOptions opts;
  opts.SetIonicState(OEFreeFormIonicState::Uncharged);
  OEFreeFormSolvResults res;

  OEGraphMol omol;
  if (!OEEstimateSolvFreeEnergy(res, omol, mol, opts))
    OEThrow.Error("Failed to calculate solvation free energy for molecule %s",
                  mol.GetTitle());

  double solvenergy = res.GetSolvationFreeEnergy();
  OEThrow.Info("Solvation Free energy for compound %s is %6.2f kcal/mol",
               mol.GetTitle(), solvenergy);

  OEWriteMolecule(ofs, omol);

  return 0;
}

Download code

freeformsolv.cpp

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.

Listing 17: Simple Free Energy Estimation

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

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

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

  OEMol mol;
  OEReadMolecule(ifs, mol);

  OEFreeFormConfOptions opts;
  OEFreeFormConf ffconf(opts);

  OEMol omol(mol);
  if(ffconf.EstimateFreeEnergies(omol)!=OEFreeFormReturnCode::Success)
    OEThrow.Error("Failed to estimate conformations free energies");

  OEFreeFormConfResults res(omol);
  OEThrow.Info("Number of unique conformations: %d", res.GetNumUniqueConfs());
  OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
  OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
  OEIter<OESingleConfResult> r;
  for (r = res.GetResultsForConformations(); r; ++r)
    OEThrow.Info("%2d %10.2f %14.2f", r->GetConfIdx(), r->GetDeltaG(), 
                 r->GetVibrationalEntropy());

  OEWriteMolecule(ofs, omol);

  return 0;
}

Download code

freeformconf.cpp

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.

Listing 18: Simple Restriction Energy Estimation

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

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

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

  OEMol mol;
  OEReadMolecule(ifs, mol);

  OEFreeFormConfOptions opts;
  OEFreeFormConf ffconf(opts);

  OEMol omol(mol);
  OEMol rmol(mol);
  if(ffconf.EstimateFreeEnergies(omol, rmol)!=OEFreeFormReturnCode::Success)
    OEThrow.Error("Failed to estimate conformations free energies");

  OEFreeFormConfResults res(omol);
  OEThrow.Info("Number of unique conformations: %d", res.GetNumUniqueConfs());
  OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
  OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
  OEIter<OESingleConfResult> r;
  for (r = res.GetResultsForConformations(); r; ++r)
    OEThrow.Info("%2d %10.2f %14.2f", r->GetConfIdx(), r->GetDeltaG(), 
                 r->GetVibrationalEntropy());

  OERestrictionEnergyResult rstrRes(rmol);
  OEThrow.Info("Global strain: %f", rstrRes.GetGlobalStrain());
  OEThrow.Info("Local strain: %f", rstrRes.GetLocalStrain());
  OEWriteMolecule(ofs, omol);

  return 0;
}

Download code

freeformconfrstr.cpp

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.

Listing 19: Advanced Free Energy Estimation

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

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

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

  OEMol mol;
  OEReadMolecule(ifs, mol);

  OEFreeFormConfOptions opts;
  OEFreeFormConfAdvanced ffconf(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(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)
    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)
    OEThrow.Error("Pre-optimization of the ensembles failed");

  // Remove duplicates from the pre-optimized ensemble
  if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode::Success)
    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)
    OEThrow.Error("Optimization of the ensembles failed");

  // Remove duplicates to obtain the set of minimum energy conformers
  if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode::Success)
    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)
    OEThrow.Error("Estimation of FreeFormConf energies failed");

  // Gather results of calculation into a results object for ease of viewing, etc.
  OEFreeFormConfResults res(omol);
  OEThrow.Info("Number of unique conformations: %d", res.GetNumUniqueConfs());
  OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
  OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
  OEIter<OESingleConfResult> r;
  for (r = res.GetResultsForConformations(); r; ++r)
    OEThrow.Info("%2d %10.2f %14.2f", r->GetConfIdx(), r->GetDeltaG(), 
                 r->GetVibrationalEntropy());

  OEWriteMolecule(ofs, omol);
  return 0;
}

Download code

freeformconfadvanced.cpp

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.

Listing 20: Advanced Restriction Energy Estimation

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

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

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

  OEMol mol;
  OEReadMolecule(ifs, mol);

  OEFreeFormConfOptions opts;
  OEFreeFormConfAdvanced ffconf(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(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(mol);
  OEMol fmol(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)
    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)
    OEThrow.Error("Pre-optimization of the ensembles failed");

  // Remove duplicates from the pre-optimized ensemble
  if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode::Success)
    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)
    OEThrow.Error("Optimization of the ensembles failed");

  // Remove duplicates to obtain the set of minimum energy conformers
  if(ffconf.RemoveDuplicates(omol)!=OEFreeFormReturnCode::Success)
    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)
    OEThrow.Error("Estimation of FreeFormConf energies failed");

  // Gather results of calculation into a results object for ease of viewing, etc.
  OEFreeFormConfResults res(omol);
  OEThrow.Info("Number of unique conformations: %d", res.GetNumUniqueConfs());
  OEThrow.Info("Conf.  Delta_G   Vibrational_Entropy");
  OEThrow.Info("      [kcal/mol]     [J/(mol K)]");
  OEIter<OESingleConfResult> r;
  for (r = res.GetResultsForConformations(); r; ++r)
    OEThrow.Info("%2d %10.2f %14.2f", 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)
	  OEThrow.Error("Identification 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)
	  OEThrow.Error("Restriction energy estimation failed");

  // Gather restriction energies into a results object for ease of viewing, etc.
  OERestrictionEnergyResult rstrRes(fmol);
  OEThrow.Info("Global strain: %f", rstrRes.GetGlobalStrain());
  OEThrow.Info("Local strain: %f", 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)
    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)
	  OEThrow.Error("Restriction energy estimation failed");

  // Gather restriction energies into a results object for ease of viewing, etc.
  OERestrictionEnergyResult rstrRes2(fmol);
  OEThrow.Info("Global strain: %f", rstrRes2.GetGlobalStrain());
  OEThrow.Info("Local strain: %f", rstrRes2.GetLocalStrain());
  
  OEWriteMolecule(ofs, omol);
  return 0;
}

Download code

freeformconfadvrstr.cpp

Finding similar conformers

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

Listing 21: Finding Similar Conformers

/* 
(C) 2022 Cadence Design Systems, Inc. (Cadence) 
All rights reserved.
TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of Cadence products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. Cadence claims no rights to Customer's
modifications. Modification of Sample Code is at Customer's sole and
exclusive risk. Sample Code may require Customer to have a then
current license or subscription to the applicable Cadence offering.
THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED.  OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"
#include "oechem.h"
#include "oeszybki.h"

using namespace OESystem;
using namespace OEChem;
using namespace OESz;

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

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

  OEMol mol;
  OEReadMolecule(ifs, mol);

  OEFreeFormConfOptions opts;
  OEFreeFormConf ffconf(opts);

  // Estimate free energies to ontain the minimum energy conformers
  OEMol omol(mol);
  if(ffconf.EstimateFreeEnergies(omol)!=OEFreeFormReturnCode::Success)
    OEThrow.Error("Failed to estimate conformations free energies");

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

  OEWriteMolecule(ofs, fmol);

  return 0;
}

Download code

findsimilarconfs.cpp

Torsion scan

The following code shows how to scan all torsions in the input molecule