3.6 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 SetFocusTarget() function.

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 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 quadradic 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.

 1 /*******************************************************************
 2 Copyright (C) 2006, 2007, 2008 OpenEye Scientific Software, Inc.
 3 *******************************************************************/
 4 #include "openeye.h"
 5
 6 #include "oezap.h"
 7 #include "oechem.h"
 8 #include "oesystem.h"
 9 #include "oeplatform.h"
10
11 using namespace OEPB;
12 using namespace OEChem;
13 using namespace OESystem;
14 using namespace OEPlatform;
15
16 void PrintHeader(const char *protTitle,const char *ligTitle)
17 {
18   OEThrow.Info("\nBinding Energy and Wall Clock Time for %s and %s",
19                protTitle, ligTitle);
20   OEThrow.Info("%15s %15s %10s", "Focussed?", "Energy(kT)", "Time(s)");
21 }
22
23 void PrintInfo(const char *focussed, double energy, double time)
24 {
25   OEThrow.Info("%15s %15.3f %10.1f", focussed, energy, time);
26 }
27
28 void CalcBindingEnergy(OEZap &zap, const OEGraphMol &protein,
29                        const OEGraphMol &ligand, const OEGraphMol &cmplx)
30 {
31   double proteinEnergy, ligandEnergy, cmplxEnergy, energy, time;
32
33   OEStopwatch stopwatch;
34   stopwatch.Start();
35
36   float *ppot = new float[protein.GetMaxAtomIdx()];
37   zap.SetMolecule(protein);
38   zap.CalcAtomPotentials(ppot);
39   proteinEnergy = 0.0f;
40   for (OEIter<OEAtomBase> atom=protein.GetAtoms();atom;++atom)
41     proteinEnergy+=ppot[atom->GetIdx()]*(float)atom->GetPartialCharge();
42   proteinEnergy *= 0.5f;
43   delete ppot;
44
45   float *lpot = new float[ligand.GetMaxAtomIdx()];
46   zap.SetMolecule(ligand);
47   zap.CalcAtomPotentials(lpot);
48   ligandEnergy = 0.0f;
49   for (OEIter<OEAtomBase> atom=ligand.GetAtoms();atom;++atom)
50     ligandEnergy+=lpot[atom->GetIdx()]*(float)atom->GetPartialCharge();
51   ligandEnergy *= 0.5f;
52   delete lpot;
53
54   float *cpot = new float[cmplx.GetMaxAtomIdx()];
55   zap.SetMolecule(cmplx);
56   zap.CalcAtomPotentials(cpot);
57   cmplxEnergy = 0.0f;
58   for (OEIter<OEAtomBase> atom=cmplx.GetAtoms();atom;++atom)
59     cmplxEnergy+=cpot[atom->GetIdx()]*(float)atom->GetPartialCharge();
60   cmplxEnergy *= 0.5f;
61   delete cpot;
62
63   energy = cmplxEnergy - ligandEnergy - proteinEnergy;
64   time = stopwatch.Elapsed();
65
66   const char *focussed;
67   if (zap.IsFocusTargetSet())
68     focussed = "Yes";
69   else
70     focussed = "No";
71
72   PrintInfo(focussed, energy, time);
73 }
74
75 int main(int argc, char *argv[])
76 {
77   if (argc!=3)
78     OEThrow.Usage("%s <protein> <ligand>", argv[0]);
79
80   oemolistream ifs;
81   if (!ifs.open(argv[1]))
82     OEThrow.Fatal("Could not open %s for reading", argv[1]);
83   OEGraphMol protein;
84   OEReadMolecule(ifs,protein);
85
86   if (!ifs.open(argv[2]))
87     OEThrow.Fatal("Could not open %s for reading", argv[2]);
88   OEGraphMol ligand;
89   OEReadMolecule(ifs,ligand);
90
91   OEAssignBondiVdWRadii(protein);
92   OEMMFFAtomTypes(protein);
93   OEMMFF94PartialCharges(protein);
94
95   OEAssignBondiVdWRadii(ligand);
96   OEMMFFAtomTypes(ligand);
97   OEMMFF94PartialCharges(ligand);
98
99   OEGraphMol cmplx(protein);
100   OEAddMols(cmplx, ligand);
101
102   float epsin = 1.0f;
103   float spacing = 0.5f;
104   OEZap zap;
105   zap.SetInnerDielectric(epsin);
106   zap.SetGridSpacing(spacing);
107
108   PrintHeader(protein.GetTitle(), ligand.GetTitle());
109
110   CalcBindingEnergy(zap, protein, ligand, cmplx);
111   zap.SetFocusTarget(ligand);
112   CalcBindingEnergy(zap, protein, ligand, cmplx);
113
114   return 0;
115 }
116

Listing:3.9 ZAP Focussing Example