The Way of Zap

There follows a dozen or so examples of using ZAP. The list is always growing because the toolkit philosophy really does work. Most of the examples use OEInterface to provide consistent command-line argument handling.

Partial Charges and Atomic Radii

There are three essential quantities for any PB calculation: coordinates, dielectric model (including radii), and charge. Coordinates will come from the atomic coordinates of the input molecule, but partial charges and the dielectric model deserve some extra attention.

Charges and radii can be calculated at run-time using functions provided by OEChem or other OpenEye toolkits (AM1-BCC charges from Quacpac TK for example). The dielectric model consists of both the value of the inner and outer dielectric, as well as how the dielectric adjusts from one to the other. ZAP uses atomic radii as the basis for its dielectric algorithm. By default, we recommend using Bondi VdW radii which can be set on a molecule using the function OEAssignBondiVdWRadii Or you can set specific radii on an atom-by-atom basis using the OEAtomBase::SetRadius function.

The current recommended value for the inner and outer dielectrics for a molecule in an aqueous solution are 1.0 and 80.0. Therefore, these are the default values for the OEZap and OEBind implementations. The historic value for the inner dielectric is 2.0, and this value may be set using the SetInnerDielectric method. Additionally, the default value for the inner dielectric of the OEET implementation remains at 2.0 for the time being.

Accurate PB calculations also depend on quality atomic charges. In OEChem, we provide MMFF94 partial charges as the lowest level of charges that we consider usable in PB. The examples that follow will use these charges, so that they don’t depend on any toolkits beside ZAP and OEChem. But for production use, we recommend AM1-BCC charges

If your input file format can store partial charges and/or radii, you can use those values instead of calculating them at run time. OEB files can store radii and partial charges, so they serve as a useful intermediate file between many OpenEye programs. Note that charges and radii will only be written to an OEB file if they are present on the molecule when written.

Atomic charges can come in a PDB file, along with atomic radii, in the form of the extra fields added after the coordinates. This style originated with DelPhi. Note that to get OEChem to read charges and radii from these fields in a PDB file, you need to set a specific flavor on the input oemolistream before reading any molecules from the stream.

Example of using OEChem to read charges and radii from PDB

oemolistream ifs;
ifs.SetFlavor(OEFormat::PDB, OEIFlavor::PDB::Default | OEIFlavor::PDB::DELPHI);

Grids

ZAP uses regular cubic lattices, or grids, to solve the PB equation. The examples here illustrate how to retrieve such from the calculation and to manipulate and store the information held by each. Typically, grids are not used per se but information is extracted, such as the potential at particular points in space, as is illustrated in the next section. However, some programs, such as VIDA, can read grids and display their properties. Also, having direct access to grids allows for manipulations of the potential map that may not have been anticipated.

This first example shows the simplest method of generating a grid of potentials using ZAP. We save the grid in OpenEye GRD (*.grd) format which is a compact, binary format that can be visualized in VIDA. Alternatively, the grid can be written into GRASP format (*.phi), which means that the grid stored will be 65 cubed regardless of the size of the grid used in the calculation. A \(65^3\) grid is obtained by interpolation to a grid with that many points that fits over the largest dimension of the grid calculated.

Listing 1: Calculating a potential grid

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oegrid.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

int main(int argc, char *argv[])
{  
  if (argc!=2)
    OEThrow.Usage("%s <molfile>", argv[0]);
  
  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);

  float epsin = 1.0f;
  OEGraphMol mol;

  OEReadMolecule(ifs, mol);  
  OEAssignBondiVdWRadii(mol);
  OEMMFFAtomTypes(mol);
  OEMMFF94PartialCharges(mol);

  OEZap zap;
  zap.SetInnerDielectric(epsin);
  zap.SetMolecule(mol);

  OEScalarGrid grid;
  if (zap.CalcPotentialGrid(grid))
    OEWriteGrid("zap.grd", grid);

  return 0;
}

The next example is an elaboration of the previous simple version where we add in control of the parameters of the calculation. Options are provided to set the internal and external (solute and solvent) dielectric constants, the distance between the molecule and the edges of the grid (boundary spacing or buffer) and the grid spacing. A smaller grid spacing implies a more dense and accurate grid, but it does come with a larger memory footprint.

Listing 2: Calculating potential grid with optional parameters

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oegrid.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

#include "zap_grid2.itf"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;
using namespace std;

bool SetupInterface(int argc, char *argv[], OEInterface &itf)
{
  OEConfigure(itf, InterfaceData);
  if(OECheckHelp(itf, argc, argv))
    return false;
  if (!OEParseCommandLine(itf, argc, argv))
    return false;

  string infile=itf.Get<string>("-in");
  if(!OEIsReadable(infile))
  { 
    OEThrow.Warning("%s is not a readable input file", infile.c_str());
    return false;
  }

  string outfile=itf.Get<string>("-out");
  if(!OEIsWriteableGrid(outfile))
  {
    OEThrow.Warning("%s is not a writable grid file", outfile.c_str());
    return false;
  }
  return true;
}

int main(int argc, char *argv[])
{  
  OEInterface itf;
  if (!SetupInterface(argc, argv, itf))
    return 1;

  OEZap zap;
  zap.SetInnerDielectric(itf.Get<float>("-epsin"));
  zap.SetOuterDielectric(itf.Get<float>("-epsout"));
  zap.SetGridSpacing(itf.Get<float>("-grid_spacing"));
  zap.SetBoundarySpacing(itf.Get<float>("-buffer"));

  oemolistream ifs;
  if (!ifs.open(itf.Get<std::string>("-in")))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);

  OEGraphMol mol;

  OEReadMolecule(ifs,mol);  
  OEAssignBondiVdWRadii(mol);
  OEMMFFAtomTypes(mol);
  OEMMFF94PartialCharges(mol);
  zap.SetMolecule(mol);

  OEScalarGrid grid;
  if (zap.CalcPotentialGrid(grid))
    OEWriteGrid(itf.Get<string>("-out"), grid);

  return 0;
}

Here we show how to calculate a difference-map, that is to say the potential difference between a standard, two dielectric, calculation and a single dielectric calculation (approximating pure Coulombic potentials). These difference potentials represent the electrostatic response of the solvent to the charges within the solute molecule. If we mask away everything outside the molecule, we can see the contributions from the charges inside the molecule.

Listing 3: Calculating potential grid with optional parameters

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oegrid.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

int main(int argc, char *argv[])
{  
  if (argc!=2)
    OEThrow.Usage("%s <molfile>", argv[0]);
  
  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);

  OEGraphMol mol;

  OEReadMolecule(ifs, mol);  
  OEAssignBondiVdWRadii(mol);
  OEMMFFAtomTypes(mol);
  OEMMFF94PartialCharges(mol);

  OEZap zap;
  zap.SetInnerDielectric(1.0f);
  zap.SetGridSpacing(0.5f);
  zap.SetMolecule(mol);

  OEScalarGrid grid1, grid2;
  // calculate standard 2-dielectric grid
  zap.CalcPotentialGrid(grid1);
  
  // calculate grid with single dielectric
  zap.SetOuterDielectric(zap.GetInnerDielectric());
  zap.CalcPotentialGrid(grid2);

  // take the difference
  OESubtractGrid(grid1, grid2);

  // mask out everything outside the molecule
  OEMaskGridByMolecule(grid1, mol, OEGridMaskType::GaussianMinus);

  OEWriteGrid("zap_diff.grd", grid1);

  return 0;
}

Atom Potentials

The potentials at any charge, as calculated by PB, can be decomposed into three forms

  1. Induced solvent potential: the potential produced by the polarization of the solvent.
  2. Inter-charge Coulombic energy: the potential that would be produced if the solvent had the same dielectric as the solute molecule from all other charges.
  3. Self potential: the potential produced by a charge at itself. Of course, as mentioned above, this is infinite analytically, but PB will produce an “approximation” to this infinity because the grid spacing is finite.

This example program shows how to calculate each of these quantities. In addition to command line flags to control dielectric, grid spacing, etc., there are three flags that affect the type of potentials calculated:

  • -calc_type solvent_only
  • -calc_type remove_self
  • -calc_type coulombic

With none of these flags, the code executes a single PB calculation, interpolates the potentials from the grid produced, reports the sum of these potentials multiplied by the charge on the respective atoms and outputs these potentials and charges in a table. This potential corresponds to (a)+(b)+(b) above.

Using the -calc_type solvent_only option executes two PB calculations, the standard calculation plus a second calculation where the external dielectric has been set to the solute dielectric. The atom potentials are then formed from the difference of these two calculations such that the remaining potential is that produced by the solvent alone. The sum of this set of atomic potentials multiplied by atomic charges all multiplied by 0.5 is the electrostatic component of transferring from a media of the same dielectric as the solute to water. These potentials correspond to (a) above.

The -calc_type remove_self flag for the example program below executes an internal algorithm in the ZAP toolkit that extracts from each atom potential that potential produced by that atom, if it is charged. This it removes the artifactual energy from the grid, i.e. energy that is not actually physical but a manifestation of the finite resolution of the grids used. The remaining potential is then that from the solvent and that from all the other charges, that is (a)+(b) above.

The -calc_type coulombic flag actually prevents any PB calculation. It merely uses Coulomb’s Law to calculate the potential (b), the inner-atomic Coulombic potential.

Listing 4: Calculating atom potentials

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

#include "zap_atompot.itf"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

static void Output(const OEMolBase &mol, float *apot, bool table)
{ 
  OEThrow.Info("Title: %s", mol.GetTitle());
  if (table)
  {
    OEThrow.Info("Atom potentials");
    OEThrow.Info("Index  Elem    Charge     Potential");
  }
  double energy=0.0;
  for (OEIter<OEAtomBase> atom=mol.GetAtoms();atom;++atom)
  {
    energy += atom->GetPartialCharge()*apot[atom->GetIdx()]; 
    if (table)
      OEThrow.Info("%3d     %2s     %6.3f     %8.3f",
                   atom->GetIdx(),OEGetAtomicSymbol(atom->GetAtomicNum()),
                   atom->GetPartialCharge(),apot[atom->GetIdx()]);
  }
  OEThrow.Info("Sum of {Potential * Charge over all atoms * 0.5} in kT = %f\n",
               0.5*energy);
}

static bool SetupInterface(int argc, char *argv[], OEInterface &itf)
{
  OEConfigure(itf, InterfaceData);
  if(OECheckHelp(itf, argc, argv))
    return false;
  if (!OEParseCommandLine(itf, argc, argv))
    return false;
  std::string filename=itf.Get<std::string>("-in");
  if(!OEIsReadable(OEGetFileType(OEGetFileExtension(filename.c_str()))))
  { 
    OEThrow.Warning("%s is not a readable input file", filename.c_str());
    return false;
  }
  return true;
}
      
int main(int argc, char *argv[])
{  
  OEInterface itf;
  if (!SetupInterface(argc, argv, itf))
    return 1;

  OEGraphMol mol;

  oemolistream ifs;
  if (!ifs.open(itf.Get<std::string>("-in")))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);

  OEReadMolecule(ifs,mol);
  OEAssignBondiVdWRadii(mol);

  if (!itf.Get<bool>("-file_charges"))
  {
    OEMMFFAtomTypes(mol);
    OEMMFF94PartialCharges(mol);
  }

  OEZap zap;
  zap.SetMolecule(mol);
  zap.SetInnerDielectric(itf.Get<float>("-epsin"));
  zap.SetBoundarySpacing(itf.Get<float>("-boundary"));
  zap.SetGridSpacing(itf.Get<float>("-grid_spacing"));

  std::string calcType = itf.Get<std::string>("-calc_type");
  if (calcType=="default")
  {
    float *apot = new float[mol.GetMaxAtomIdx()];
    zap.CalcAtomPotentials(apot);
    Output(mol,apot,itf.Get<bool>("-atomtable"));
    delete [] apot;
  }
  else if (calcType=="solvent_only")
  {
    float *apot = new float[mol.GetMaxAtomIdx()];
    zap.CalcAtomPotentials(apot);

    float *apot2 = new float[mol.GetMaxAtomIdx()];
    zap.SetOuterDielectric(zap.GetInnerDielectric());
    zap.CalcAtomPotentials(apot2);

    // find the difference
    for (OEIter<OEAtomBase> atom=mol.GetAtoms();atom;++atom)
      apot[atom->GetIdx()] -= apot2[atom->GetIdx()];

    Output(mol,apot,itf.Get<bool>("-atomtable"));

    delete [] apot;
    delete [] apot2;
  }
  else if (calcType=="remove_self")
  {    
    float *apot = new float[mol.GetMaxAtomIdx()];
    zap.CalcAtomPotentials(apot,true);
    Output(mol,apot,itf.Get<bool>("-atomtable"));
    delete [] apot;
  }
  else if (calcType=="coulombic")
  {
    float epsin = itf.Get<float>("-epsin");
    float x = OECoulombicSelfEnergy(mol, epsin);
    OEThrow.Info("Coulombic Assembly Energy ");
    OEThrow.Info(" = Sum of {Potential * Charge over all atoms * 0.5} in kT "
                 "= %f",x);

    float *apot = new float[mol.GetMaxAtomIdx()];
    OECoulombicAtomPotentials(mol, epsin, apot);

    Output(mol,apot,itf.Get<bool>("-atomtable"));
    delete [] apot;
  }

  return 0;
}

Solvation Energies: PBSA

A principle use of PB is to calculate the energy stored in the exterior dielectric, for example the partial alignment of water dipoles. This corresponds to the difference between a calculation with uniform dielectric and one with a different external dielectric. In the examples for atom potentials this would be the sum of the solvent potentials at each atom, multiplied by the charges on that atom all multiplied by 0.5.

However, we know water also has an energy component due to hydrophobicity, which is typically estimated as a constant (approximately 10-25 calories per square Angstrom) multiplied by the accessible area of the molecule. This value is an open scientific question. We recommend using 10 for vacuum-water transfer energies and 25 for protein calculations, but ultimately your own scientific judgment should be used. These examples include an area calculation. Taken together these numbers are an approximation of the transfer energy from a low dielectric medium (alkane solvent, protein interior) into water.

The following examples use the OEArea object to calculate the accessible surface area of the molecule. The area is calculated using the same grid based Gaussians that ZAP uses. Therefore, it will not give the same results as triangle summation methods for calculating surface area, which are used in OESpicoli.

Listing 5: Calculating solvation energies

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

static const float KCalsPerKT = 0.59f;
static const float KCalsPerSqAngstrom = 0.025f;

void PrintHeader()
{
  OEThrow.Info("%-12s %12s %12s %12s %14s", "Title", "Solv(kcal)", 
               "Area(Ang^2)", "Total(kcal)", "Int.Coul(kcal)");
}

void PrintLine(const char *title, float solv, float area, float coul)
{
  float total = KCalsPerKT*solv + KCalsPerSqAngstrom*area;
  OEThrow.Info("%-12s %12.2f %12.2f %12.2f %14.2f\n",
         title, KCalsPerKT*solv, area, total, KCalsPerKT*coul);
}

int main(int argc, char *argv[])
{  
  float epsin = 1.0f;

  if (argc!=2)
    OEThrow.Usage("%s <molfile>", argv[0]);
  
  OEGraphMol mol;
  OEZap zap;
  zap.SetInnerDielectric(epsin);
  zap.SetGridSpacing(0.5f);

  OEArea area;

  PrintHeader();
  float solv, aval, coul;

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

  while (OEReadMolecule(ifs,mol))
  {
    OEAssignBondiVdWRadii(mol);
    OEMMFFAtomTypes(mol);
    OEMMFF94PartialCharges(mol);
    zap.SetMolecule(mol);
    solv = zap.CalcSolvationEnergy();
    aval = area.GetArea(mol);
    coul = OECoulombicSelfEnergy(mol, epsin);
    PrintLine(mol.GetTitle(), solv, aval, coul);
  }

  return 0;
}

This next solvation example calculates the transfer energy from vacuum to water.

Listing 6: Calculating vacuum-water transfer energies

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

static const float KCalsPerKT = 0.59f;
static const float KCalsPerSqAngstrom = 0.010f;

int main(int argc, char *argv[])
{  
  float epsin = 1.0f;

  if (argc!=2)
    OEThrow.Usage("%s <molfile>", argv[0]);
  
  OEGraphMol mol;
  OEZap zap;
  zap.SetInnerDielectric(epsin);
  zap.SetGridSpacing(0.5f);
  OEArea area;

  OEThrow.Info("Title              Vacuum->Water(kcal)\n");
  float solv, aval;

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

  while (OEReadMolecule(ifs,mol))
  {
    OEAssignBondiVdWRadii(mol);
    OEMMFFAtomTypes(mol);
    OEMMFF94PartialCharges(mol);
    zap.SetMolecule(mol);
    solv = zap.CalcSolvationEnergy();
    aval = area.GetArea(mol);
    OEThrow.Info("%-20s   %6.2f", mol.GetTitle(), 
                 KCalsPerKT*solv+KCalsPerSqAngstrom*aval);    
  }

  return 0;
}

Binding Energies

Binding energies and other binding related properties can be calculated using the OEBind class. The OEBind class serves two primary purposes: 1) to have a class in place for users to easily calculate binding related data, and 2) to show how binding related data may be calculated. OEBind is not particularly well-suited for extremely efficient calculations required by simulations. This is because some of the data available from OEBind may be considered uninteresting for the project and there is no need to spend time calculating it, or because OEBind calculates the unbound protein and ligand properties every time, which is unnecessary if the same protein and ligand appear in multiple calculations. The API section for OEBind is detailed with regard to implementation, so the user can see how to create their own class or set of function calls if needed.

There are two classes that store the results of a binding calculation, OEBindResults and OESimpleBindResults. OEBindResults is a superset of OESimpleBindResults. For an in-depth understanding of what is calculated, the electrostatics portion of the algorithm for the SimpleBind function is shown in the API section.

The following is a sample program for using the OEBind class that is provided in the distribution.

Listing 7: Bind example

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

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

  float epsin = 1.0f;  
  OEMol protein, ligand;

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

  OEReadMolecule(ifs, protein);
  OEAssignBondiVdWRadii(protein);
  OEMMFFAtomTypes(protein);
  OEMMFF94PartialCharges(protein);

  oeerr << "protein:   " << protein.GetTitle() << oeendl;

  OEBind bind;
  bind.GetZap().SetInnerDielectric(epsin);
  bind.SetProtein(protein);
  OEBindResults results;

  if (!ifs.open(argv[2]))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);

  ifs.SetConfTest(OEIsomericConfTest());
  while (OEReadMolecule(ifs,ligand))
  {
    OEAssignBondiVdWRadii(ligand);
    OEMMFFAtomTypes(ligand);
    OEMMFF94PartialCharges(ligand);

    oeerr << "ligand:  " << ligand.GetTitle() << " has " 
          << ligand.NumConfs() << " conformers" << oeendl;

    for (OEIter<OEConfBase> conf=ligand.GetConfs();conf;++conf)
    {
      bind.Bind(conf, results);
      oeerr << " conf# " << conf->GetIdx()
            << "   be = " << results.GetBindingEnergy() << oeendl;
    }
    oeerr << oeendl;
  }

  return 0;
}

Focussing

Focussing is a way to achieve the desired precision around a target (such as a ligand) while maintaining reasonable time and memory limits for the calculation. The target for focussing is set using the OEZap::SetFocusTarget method.

For a typical ZAP electrostatics calculation, a consistent grid spacing is used for the entire system, where the default value is 0.5 Angstroms but it may be set to a custom value using OEZap::SetGridSpacing. This method is entirely appropriate for certain calculations, such as solvation energy calculations for small molecules. For other types of calculations, such as a binding energy calculation for a protein and ligand, a consistent grid spacing may cause the calculation to be either prohibitively expensive or insufficiently precise around the binding area.

Focussing alleviates this problem by using a fine grid for the target volume, and coarser grids away from the target. The grid spacing setting for ZAP is applied to the target volume, and the grid spacing is doubled for each addition coarse grid surrounding the target. A quadratic interpolation is used for the grid intersections. The implementation of the bind uses the SetFocusTarget() method to set the ligand as the target for focussing.

The following example program computes the binding energy of a protein and ligand twice, once without focussing and once with the ligand as the focussing target. The values of the binding energy and the time it took to calculate them are displayed. For a focussed atom potential calculation, the full electrostatics for the protein and complex would each be computed with multiple electrostatics calculations. The first calculation would create the potential grid for the target volume with a grid spacing of 0.5 (assuming the default spacing is being used). An additional calculation with a grid spacing of 1.0 would be done on a larger volume volume, and then additional calculations with grid spacings of 2.0, 4.0, and so on would be done until the grid is large enough to contain the entire volume of the system.

Listing 8: Focussing example

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

void PrintHeader(const char *protTitle,const char *ligTitle)
{
  OEThrow.Info("\nBinding Energy and Wall Clock Time for %s and %s", 
               protTitle, ligTitle);
  OEThrow.Info("%15s %15s %10s", "Focussed?", "Energy(kT)", "Time(s)");
}

void PrintInfo(const char *focussed, double energy, double time)
{
  OEThrow.Info("%15s %15.3f %10.1f", focussed, energy, time);
}

void CalcBindingEnergy(OEZap &zap, const OEGraphMol &protein, 
                       const OEGraphMol &ligand, const OEGraphMol &cmplx)
{
  double proteinEnergy, ligandEnergy, cmplxEnergy, energy, time;

  OEStopwatch stopwatch;
  stopwatch.Start();

  float *ppot = new float[protein.GetMaxAtomIdx()];
  zap.SetMolecule(protein);
  zap.CalcAtomPotentials(ppot);
  proteinEnergy = 0.0f;
  for (OEIter<OEAtomBase> atom=protein.GetAtoms();atom;++atom)
    proteinEnergy+=ppot[atom->GetIdx()]*(float)atom->GetPartialCharge();
  proteinEnergy *= 0.5f;
  delete[] ppot;

  float *lpot = new float[ligand.GetMaxAtomIdx()];
  zap.SetMolecule(ligand);
  zap.CalcAtomPotentials(lpot);
  ligandEnergy = 0.0f;
  for (OEIter<OEAtomBase> atom=ligand.GetAtoms();atom;++atom)
    ligandEnergy+=lpot[atom->GetIdx()]*(float)atom->GetPartialCharge();
  ligandEnergy *= 0.5f;
  delete[] lpot;

  float *cpot = new float[cmplx.GetMaxAtomIdx()];
  zap.SetMolecule(cmplx);
  zap.CalcAtomPotentials(cpot);
  cmplxEnergy = 0.0f;
  for (OEIter<OEAtomBase> atom=cmplx.GetAtoms();atom;++atom)
    cmplxEnergy+=cpot[atom->GetIdx()]*(float)atom->GetPartialCharge();
  cmplxEnergy *= 0.5f;
  delete[] cpot;

  energy = cmplxEnergy - ligandEnergy - proteinEnergy;
  time = stopwatch.Elapsed();

  const char *focussed;
  if (zap.IsFocusTargetSet())
    focussed = "Yes";
  else
    focussed = "No";

  PrintInfo(focussed, energy, time);
}

int main(int argc, char *argv[])
{  
  if (argc!=3)
    OEThrow.Usage("%s <protein> <ligand>", argv[0]);
 
  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);
  OEGraphMol protein;
  OEReadMolecule(ifs,protein);

  if (!ifs.open(argv[2]))
    OEThrow.Fatal("Could not open %s for reading", argv[2]);
  OEGraphMol ligand;
  OEReadMolecule(ifs,ligand);

  OEAssignBondiVdWRadii(protein);
  OEMMFFAtomTypes(protein);
  OEMMFF94PartialCharges(protein);

  OEAssignBondiVdWRadii(ligand);
  OEMMFFAtomTypes(ligand);
  OEMMFF94PartialCharges(ligand);

  OEGraphMol cmplx(protein);
  OEAddMols(cmplx, ligand);

  float epsin = 1.0f;
  float spacing = 0.5f;
  OEZap zap;
  zap.SetInnerDielectric(epsin);
  zap.SetGridSpacing(spacing);

  PrintHeader(protein.GetTitle(), ligand.GetTitle());

  CalcBindingEnergy(zap, protein, ligand, cmplx);
  zap.SetFocusTarget(ligand);
  CalcBindingEnergy(zap, protein, ligand, cmplx);

  return 0;
}

Gradients/Forces

The solvent forces acting on the atoms in a molecule may be calculated using the OEZap::CalcForces method. The components of the forces are set to a float array of length 3N, where N is the number of atoms. The components of the gradient are equal in magnitude to the force, but have the opposite sign.

Listing 9: Calculating solvation forces

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

void PrintHeader(const char *title)
{
  OEThrow.Info("\nForce Components for %s, in kT/Angstrom",title);
  OEThrow.Info("%6s %8s %8s %8s %8s", "Index", "Element", "-dE/dx", 
               "-dE/dy", "-dE/dz");
}

void PrintForces(const OEMolBase &mol, const float *forces)
{
  OEIter<OEAtomBase> atom;
  for(atom=mol.GetAtoms(); atom; ++atom)
  {
    OEThrow.Info("%6d %8s %8.2f %8.2f %8.2f",
           atom->GetIdx(), OEGetAtomicSymbol(atom->GetAtomicNum()), 
           forces[atom->GetIdx()], forces[atom->GetIdx()+1], 
           forces[atom->GetIdx()+2]);
  }
}

int main(int argc, char *argv[])
{  
  if (argc!=2)
    OEThrow.Usage("%s <molfile>", argv[0]);
 
  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);

  float epsin = 1.0f;
  OEGraphMol mol;
  OEZap zap;
  zap.SetInnerDielectric(epsin);

  while (OEReadMolecule(ifs,mol))
  {
    PrintHeader(mol.GetTitle());
    float *forces = new float[mol.GetMaxAtomIdx()*3];
    OEAssignBondiVdWRadii(mol);
    OEMMFFAtomTypes(mol);
    OEMMFF94PartialCharges(mol);
    zap.SetMolecule(mol);
    zap.CalcForces(forces);
    PrintForces(mol, forces);
    delete [] forces;
  }
  return 0;
}

Electrostatic Similarity

Electrostatic similarity may be calculated using the OEET class. The following example program shows how to obtain the electrostatic Tanimoto between a reference molecule and a trial molecule. The electrostatic Tanimoto is affected by the partial charges of the molecules as well as the three-dimensional structure, including the tautomer state, spacial orientation, and conformation. In the example program shown below, MMFF charges are assigned to the molecules, but it is assumed that they have already been spatially aligned.

Listing 10: Calculating electrostatic tanimoto

/*
(C) 2017 OpenEye Scientific Software Inc. All rights reserved.

TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
provided to current licensees or subscribers of OpenEye products or
SaaS offerings (each a "Customer").
Customer is hereby permitted to use, copy, and modify the Sample Code,
subject to these terms. OpenEye 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 OpenEye 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 OpenEye be
liable for any damages or liability in connection with the Sample Code
or its use.
*/

#include "openeye.h"

#include "oezap.h"
#include "oechem.h"
#include "oesystem.h"
#include "oeplatform.h"

using namespace OEPB;
using namespace OEChem;
using namespace OESystem;
using namespace OEPlatform;

int main(int argc, char *argv[])
{  
  if (argc!=3)
    OEThrow.Usage("%s <reffile> <fitfile>", argv[0]);
  
  OEGraphMol refmol, fitmol;

  oemolistream ifs;
  if (!ifs.open(argv[1]))
    OEThrow.Fatal("Could not open %s for reading", argv[1]);
  OEReadMolecule(ifs,refmol);

  OEMMFFAtomTypes(refmol);
  OEMMFF94PartialCharges(refmol);
  OEAssignBondiVdWRadii(refmol);

  OEET et;
  et.SetRefMol(refmol);

  float tanimoto;

  if (!ifs.open(argv[2]))
    OEThrow.Fatal("Could not open %s for reading", argv[2]);

  while (OEReadMolecule(ifs, fitmol))
  {
    OEMMFFAtomTypes(fitmol);
    OEMMFF94PartialCharges(fitmol);
    OEAssignBondiVdWRadii(fitmol);

    tanimoto=et.Tanimoto(fitmol);
    OEThrow.Info("Title: %s, ET %4.2f", fitmol.GetTitle(), tanimoto);
  }

  return 0;
}