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

ifs = oemolistream()
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 irregardless 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

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
# zap_grid1.py
#############################################################################

import sys
from openeye.oechem import *
from openeye.oegrid import *
from openeye.oezap import *


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

    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])
    mol = OEGraphMol()
    OEReadMolecule(ifs, mol)
    OEAssignBondiVdWRadii(mol)
    OEMMFFAtomTypes(mol)
    OEMMFF94PartialCharges(mol)

    epsin = 1.0
    zap = OEZap()
    zap.SetInnerDielectric(epsin)
    zap.SetMolecule(mol)

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


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

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

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
# zap_grid2.py
#############################################################################

import sys
from openeye.oechem import *
from openeye.oegrid import *
from openeye.oezap import *


def main(argv=[__name__]):
    itf = OEInterface()
    if not SetupInterface(argv, itf):
        return 1

    zap = OEZap()
    zap.SetInnerDielectric(itf.GetFloat("-epsin"))
    zap.SetOuterDielectric(itf.GetFloat("-epsout"))
    zap.SetGridSpacing(itf.GetFloat("-grid_spacing"))
    zap.SetBoundarySpacing(itf.GetFloat("-buffer"))

    mol = OEGraphMol()
    ifs = oemolistream()
    if not ifs.open(itf.GetString("-in")):
        OEThrow.Fatal("Unable to open %s for reading" % itf.GetString("-in"))
    OEReadMolecule(ifs, mol)
    OEAssignBondiVdWRadii(mol)
    OEMMFFAtomTypes(mol)
    OEMMFF94PartialCharges(mol)

    zap.SetMolecule(mol)

    grid = OEScalarGrid()
    if zap.CalcPotentialGrid(grid):
        if itf.GetBool("-mask"):
            OEMaskGridByMolecule(grid, mol)
        OEWriteGrid(itf.GetString("-out"), grid)
    return 0

InterfaceData = """
#zap_grid2 interface definition

!PARAMETER -in
  !TYPE string
  !BRIEF Input molecule file
  !REQUIRED true
!END

!PARAMETER -out
  !TYPE string
  !BRIEF Output grid file
  !REQUIRED true
!END

!PARAMETER -epsin
  !TYPE float
  !BRIEF Inner dielectric
  !DEFAULT 1.0
  !LEGAL_RANGE 0.0 100.0
!END

!PARAMETER -epsout
  !TYPE float
  !BRIEF Outer dielectric
  !DEFAULT 80.0
  !LEGAL_RANGE 0.0 100.0
!END

!PARAMETER -grid_spacing
  !TYPE float
  !DEFAULT 0.5
  !BRIEF Spacing between grid points (Angstroms)
  !LEGAL_RANGE 0.1 2.0
!END

!PARAMETER -buffer
  !TYPE float
  !DEFAULT 2.0
  !BRIEF Extra buffer outside extents of molecule.
  !LEGAL_RANGE 0.1 10.0
!END

!PARAMETER -mask
  !TYPE bool
  !DEFAULT false
  !BRIEF Mask potential grid by the molecule
!END
"""


def SetupInterface(argv, itf):
    OEConfigure(itf, InterfaceData)
    if OECheckHelp(itf, argv):
        return False
    if not OEParseCommandLine(itf, argv):
        return False
    infile = itf.GetString("-in")
    if not OEIsReadable(OEGetFileType(OEGetFileExtension(infile))):
        OEThrow.Warning("%s is not a readable input file" % infile)
        return False
    outfile = itf.GetString("-out")
    if not OEIsWriteableGrid(OEGetGridFileType(OEGetFileExtension(outfile))):
        OEThrow.Warning("%s is not a writable grid file" % outfile)
        return False
    return True

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

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

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
# zap_grid3.py
#############################################################################
import sys
from openeye.oechem import *
from openeye.oegrid import *
from openeye.oezap import *


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

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

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

    zap = OEZap()
    zap.SetInnerDielectric(1.0)
    zap.SetGridSpacing(0.5)
    zap.SetMolecule(mol)

    # calculate standard 2-dielectric grid
    grid1 = OEScalarGrid()
    zap.CalcPotentialGrid(grid1)

    # calculate grid with single dielectric
    grid2 = OEScalarGrid()
    zap.SetOuterDielectric(zap.GetInnerDielectric())
    zap.CalcPotentialGrid(grid2)

    # take the difference
    OESubtractScalarGrid(grid1, grid2)

    # mask out everything outside the molecule
    OEMaskGridByMolecule(grid1, mol, OEGridMaskType_GaussianMinus)

    OEWriteGrid("zap_diff.grd", grid1)


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

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

#!/usr/bin/env python
#############################################################################
#  Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
from __future__ import print_function
import sys
from openeye.oechem import *
from openeye.oezap import *


def Output(mol, apot, showAtomTable):
    print("Title: %s" % mol.GetTitle())
    if showAtomTable:
        print("Atom potentials")
        print("Index  Elem    Charge     Potential")

    energy = 0.0
    for atom in mol.GetAtoms():
        energy += atom.GetPartialCharge()*apot[atom.GetIdx()]
        if showAtomTable:
            print("%3d     %2s     %6.3f     %8.3f" %
                  (atom.GetIdx(),
                   OEGetAtomicSymbol(atom.GetAtomicNum()),
                   atom.GetPartialCharge(),
                   apot[atom.GetIdx()]))

    print("Sum of {Potential * Charge over all atoms * 0.5} in kT = %f\n" %
          (0.5*energy))


def CalcAtomPotentials(itf):
    mol = OEGraphMol()

    ifs = oemolistream()
    if not ifs.open(itf.GetString("-in")):
        OEThrow.Fatal("Unable to open %s for reading" % itf.GetString("-in"))

    OEReadMolecule(ifs, mol)
    OEAssignBondiVdWRadii(mol)

    if not itf.GetBool("-file_charges"):
        OEMMFFAtomTypes(mol)
        OEMMFF94PartialCharges(mol)

    zap = OEZap()
    zap.SetMolecule(mol)
    zap.SetInnerDielectric(itf.GetFloat("-epsin"))
    zap.SetBoundarySpacing(itf.GetFloat("-boundary"))
    zap.SetGridSpacing(itf.GetFloat("-grid_spacing"))

    showAtomTable = itf.GetBool("-atomtable")
    calcType = itf.GetString("-calc_type")
    if calcType == "default":
        apot = OEFloatArray(mol.GetMaxAtomIdx())
        zap.CalcAtomPotentials(apot)
        Output(mol, apot, showAtomTable)

    elif calcType == "solvent_only":
        apot = OEFloatArray(mol.GetMaxAtomIdx())
        zap.CalcAtomPotentials(apot)

        apot2 = OEFloatArray(mol.GetMaxAtomIdx())
        zap.SetOuterDielectric(zap.GetInnerDielectric())
        zap.CalcAtomPotentials(apot2)

        # find the differences
        for atom in mol.GetAtoms():
            idx = atom.GetIdx()
            apot[idx] -= apot2[idx]

        Output(mol, apot, showAtomTable)

    elif calcType == "remove_self":
        apot = OEFloatArray(mol.GetMaxAtomIdx())
        zap.CalcAtomPotentials(apot, True)
        Output(mol, apot, showAtomTable)

    elif calcType == "coulombic":
        epsin = itf.GetFloat("-epsin")
        x = OECoulombicSelfEnergy(mol, epsin)
        print("Coulombic Assembly Energy")
        print("  = Sum of {Potential * Charge over all atoms * 0.5} in kT = %f" % x)
        apot = OEFloatArray(mol.GetMaxAtomIdx())
        OECoulombicAtomPotentials(mol, epsin, apot)
        Output(mol, apot, showAtomTable)

    return 0


def SetupInterface(itf, InterfaceData):
    OEConfigure(itf, InterfaceData)
    if OECheckHelp(itf, sys.argv):
        return False
    if not OEParseCommandLine(itf, sys.argv):
        return False
    return True


def main(InterfaceData):
    itf = OEInterface()
    if not SetupInterface(itf, InterfaceData):
        return 1

    return CalcAtomPotentials(itf)

InterfaceData = """
#zap_atompot interface definition

!PARAMETER -in
  !TYPE string
  !BRIEF Input molecule file.
  !REQUIRED true
  !KEYLESS 1
!END

!PARAMETER -file_charges
  !TYPE bool
  !DEFAULT false
  !BRIEF Use partial charges from input file rather than calculating with MMFF.
!END

!PARAMETER -calc_type
  !TYPE string
  !DEFAULT default
  !LEGAL_VALUE default
  !LEGAL_VALUE solvent_only
  !LEGAL_VALUE remove_self
  !LEGAL_VALUE coulombic
  !LEGAL_VALUE breakdown
  !BRIEF Choose type of atom potentials to calculate
!END

!PARAMETER -atomtable
  !TYPE bool
  !DEFAULT false
  !BRIEF Output a table of atom potentials
!END

!PARAMETER -epsin
  !TYPE float
  !BRIEF Inner dielectric
  !DEFAULT 1.0
  !LEGAL_RANGE 0.0 100.0
!END

!PARAMETER -grid_spacing
  !TYPE float
  !DEFAULT 0.5
  !BRIEF Spacing between grid points (Angstroms)
  !LEGAL_RANGE 0.1 2.0
!END

!PARAMETER -boundary
  !ALIAS -buffer
  !TYPE float
  !DEFAULT 2.0
  !BRIEF Extra buffer outside extents of molecule.
  !LEGAL_RANGE 0.1 10.0
!END
"""

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

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

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
from __future__ import print_function
import sys
from openeye.oechem import *
from openeye.oezap import *

KCalsPerKT = 0.59
KCalsPerSqAngstrom = 0.025


def PrintHeader():
    print("Title        Solv(kcal)   Area(Ang^2)   Total(kcal)  Int.Coul(kcal)\n")


def PrintLine(title, solv, area, coul):
    total = KCalsPerKT*solv + KCalsPerSqAngstrom*area
    print("%-12s   %6.2f       %6.2f       %6.2f          %6.2f" %
          (title, KCalsPerKT*solv, area, total, KCalsPerKT*coul))


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

    epsin = 1.0

    zap = OEZap()
    zap.SetInnerDielectric(epsin)
    zap.SetGridSpacing(0.5)

    area = OEArea()

    PrintHeader()
    mol = OEGraphMol()
    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to 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

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

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

Listing 6: Calculating vacuum-water transfer energies

#!/usr/bin/env python
#############################################################################
#  Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
import sys
from openeye.oechem import *
from openeye.oezap import *

KCalsPerKT = 0.59
KCalsPerSqAngstrom = 0.010


def main(argv=[__name__]):

    if len(argv) != 2:
        OEThrow.Usage("%s <molfile>" % argv[0])

    epsin = 1.0

    zap = OEZap()
    zap.SetInnerDielectric(epsin)
    zap.SetGridSpacing(0.5)

    area = OEArea()

    mol = OEGraphMol()
    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])
    OEThrow.Info("%-20s %6s\n" % ("Title", "Vacuum->Water(kcal)"))
    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

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

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

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
# bind.py
#############################################################################
from __future__ import print_function
import sys

from openeye.oechem import *
from openeye.oezap import *


def main(argv=[__name__]):
    if len(argv) != 3:
        OEThrow.Usage("%s <protein> <ligands>" % argv[0])

    protein = OEMol()

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

    OEAssignBondiVdWRadii(protein)
    OEMMFFAtomTypes(protein)
    OEMMFF94PartialCharges(protein)
    print("protein:   " + protein.GetTitle())

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

    if not ifs.open(argv[2]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[2])
    ifs.SetConfTest(OEIsomericConfTest())

    ligand = OEMol()
    while OEReadMolecule(ifs, ligand):
        OEAssignBondiVdWRadii(ligand)
        OEMMFFAtomTypes(ligand)
        OEMMFF94PartialCharges(ligand)
        print("ligand:  %s has %d conformers" %
              (ligand.GetTitle(), ligand.NumConfs()))

        for conf in ligand.GetConfs():
            bind.Bind(conf, results)
            print(" conf# %d   be = %f" %
                  (conf.GetIdx(), results.GetBindingEnergy()))
    return 0

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

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

#!/usr/bin/env python
#############################################################################
#  Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
import sys
from openeye.oechem import *
from openeye.oezap import *


def PrintHeader(protTitle, ligTitle):
    OEThrow.Info("\nBinding Energy and Wall Clock Time for %s and %s" %
                 (protTitle, ligTitle))
    OEThrow.Info("%15s %15s %10s" % ("Focused?", "Energy(kT)", "Time(s)"))


def PrintInfo(focussed, energy, time):
    OEThrow.Info("%15s %15.3f %10.1f" % (focussed, energy, time))


def CalcBindingEnergy(zap, protein, ligand, cmplx):
    stopwatch = OEStopwatch()
    stopwatch.Start()

    ppot = OEFloatArray(protein.GetMaxAtomIdx())
    zap.SetMolecule(protein)
    zap.CalcAtomPotentials(ppot)
    proteinEnergy = 0.0
    for atom in protein.GetAtoms():
        proteinEnergy += ppot[atom.GetIdx()]*atom.GetPartialCharge()
    proteinEnergy *= 0.5

    lpot = OEFloatArray(ligand.GetMaxAtomIdx())
    zap.SetMolecule(ligand)
    zap.CalcAtomPotentials(lpot)
    ligandEnergy = 0.0
    for atom in ligand.GetAtoms():
        ligandEnergy += lpot[atom.GetIdx()]*atom.GetPartialCharge()
    ligandEnergy *= 0.5

    cpot = OEFloatArray(cmplx.GetMaxAtomIdx())
    zap.SetMolecule(cmplx)
    zap.CalcAtomPotentials(cpot)
    cmplxEnergy = 0.0
    for atom in cmplx.GetAtoms():
        cmplxEnergy += cpot[atom.GetIdx()]*atom.GetPartialCharge()
    cmplxEnergy *= 0.5

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

    if zap.IsFocusTargetSet():
        focused = "Yes"
    else:
        focused = "No"

    PrintInfo(focused, energy, time)


def main(argv=[__name__]):
    if len(argv) != 3:
        OEThrow.Usage("%s <protein> <ligand>" % argv[0])

    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])
    protein = OEGraphMol()
    OEReadMolecule(ifs, protein)

    if not ifs.open(argv[2]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[2])
    ligand = OEGraphMol()
    OEReadMolecule(ifs, ligand)

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

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

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

    epsin = 1.0
    spacing = 0.5
    zap = OEZap()
    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

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

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

#!/usr/bin/env python
#############################################################################
#  Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
import sys
from openeye.oechem import *
from openeye.oezap import *


def PrintHeader(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"))


def PrintForces(mol, forces):
    for atom in mol.GetAtoms():
        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]))


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

    epsin = 1.0
    zap = OEZap()
    zap.SetInnerDielectric(epsin)

    mol = OEGraphMol()
    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])
    while OEReadMolecule(ifs, mol):
        PrintHeader(mol.GetTitle())
        forces = OEFloatArray(mol.GetMaxAtomIdx()*3)
        OEAssignBondiVdWRadii(mol)
        OEMMFFAtomTypes(mol)
        OEMMFF94PartialCharges(mol)
        zap.SetMolecule(mol)
        zap.CalcForces(forces)
        PrintForces(mol, forces)

    return 0

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

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

#!/usr/bin/env python
#############################################################################
# Copyright (C) 2006-2014 OpenEye Scientific Software, Inc.
#############################################################################
# calc_et.py
#############################################################################
import sys
from openeye.oechem import *
from openeye.oezap import *


def main(argv=[__name__]):
    if len(argv) != 3:
        OEThrow.Usage("calc_et.py <reffile> <fitfile>")

    refmol = OEGraphMol()

    ifs = oemolistream()
    if not ifs.open(argv[1]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[1])
    OEReadMolecule(ifs, refmol)
    OEAssignBondiVdWRadii(refmol)
    OEMMFFAtomTypes(refmol)
    OEMMFF94PartialCharges(refmol)

    et = OEET()
    et.SetRefMol(refmol)

    OEThrow.Info("dielectric: %.4f" % et.GetDielectric())
    OEThrow.Info("inner mask: %.4f" % et.GetInnerMask())
    OEThrow.Info("outer mask: %.4f" % et.GetOuterMask())
    OEThrow.Info("salt conc : %.4f" % et.GetSaltConcentration())
    OEThrow.Info("join      : %d" % et.GetJoin())

    if not ifs.open(argv[2]):
        OEThrow.Fatal("Unable to open %s for reading" % argv[2])

    fitmol = OEGraphMol()
    while OEReadMolecule(ifs, fitmol):
        OEAssignBondiVdWRadii(fitmol)
        OEMMFFAtomTypes(fitmol)
        OEMMFF94PartialCharges(fitmol)
        OEThrow.Info("Title: %s, ET %4.2f" %
                     (fitmol.GetTitle(), et.Tanimoto(fitmol)))
    return 0

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