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 #!/usr/bin/env python
2 #############################################################################
3 # Copyright (C) 2006, 2007, 2008 OpenEye Scientific Software, Inc.
4 #############################################################################
5 import os, sys
6 from openeye.oechem import *
7 from openeye.oezap import *
8
9 def PrintHeader(protTitle, ligTitle):
10 OEThrow.Info("\nBinding Energy and Wall Clock Time for %s and %s" %
11 (protTitle, ligTitle))
12 OEThrow.Info("%15s %15s %10s" % ("Focussed?", "Energy(kT)", "Time(s)"))
13
14 def PrintInfo(focussed, energy, time):
15 OEThrow.Info("%15s %15.3f %10.1f" % (focussed, energy, time))
16
17 def CalcBindingEnergy(zap, protein, ligand, cmplx):
18 stopwatch = OEStopwatch()
19 stopwatch.Start()
20
21 ppot = OEFloatArray(protein.GetMaxAtomIdx())
22 zap.SetMolecule(protein)
23 zap.CalcAtomPotentials(ppot)
24 proteinEnergy = 0.0
25 for atom in protein.GetAtoms():
26 proteinEnergy+=ppot[atom.GetIdx()]*atom.GetPartialCharge()
27 proteinEnergy *= 0.5
28
29 lpot = OEFloatArray(ligand.GetMaxAtomIdx())
30 zap.SetMolecule(ligand)
31 zap.CalcAtomPotentials(lpot)
32 ligandEnergy = 0.0
33 for atom in ligand.GetAtoms():
34 ligandEnergy+=lpot[atom.GetIdx()]*atom.GetPartialCharge()
35 ligandEnergy *= 0.5
36
37 cpot = OEFloatArray(cmplx.GetMaxAtomIdx())
38 zap.SetMolecule(cmplx)
39 zap.CalcAtomPotentials(cpot)
40 cmplxEnergy = 0.0
41 for atom in cmplx.GetAtoms():
42 cmplxEnergy+=cpot[atom.GetIdx()]*atom.GetPartialCharge()
43 cmplxEnergy *= 0.5
44
45 energy = cmplxEnergy - ligandEnergy - proteinEnergy
46 time = stopwatch.Elapsed()
47
48 if (zap.IsFocusTargetSet()):
49 focussed = "Yes"
50 else:
51 focussed = "No"
52
53 PrintInfo(focussed, energy, time);
54
55 def main(argv = [__name__]):
56 if len(argv) != 3:
57 OEThrow.Usage("%s <protein> <ligand>" % argv[0])
58
59 ifs = oemolistream()
60 if not ifs.open(argv[1]):
61 OEThrow.Fatal("Unable to open %s for reading" % argv[1])
62 protein = OEGraphMol()
63 OEReadMolecule(ifs, protein)
64
65 if not ifs.open(argv[2]):
66 OEThrow.Fatal("Unable to open %s for reading" % argv[2])
67 ligand = OEGraphMol()
68 OEReadMolecule(ifs, ligand)
69
70 OEAssignBondiVdWRadii(protein)
71 OEMMFFAtomTypes(protein)
72 OEMMFF94PartialCharges(protein)
73
74 OEAssignBondiVdWRadii(ligand)
75 OEMMFFAtomTypes(ligand)
76 OEMMFF94PartialCharges(ligand)
77
78 cmplx = OEGraphMol(protein)
79 OEAddMols(cmplx, ligand)
80
81 epsin = 1.0
82 spacing = 0.5
83 zap = OEZap()
84 zap.SetInnerDielectric(epsin)
85 zap.SetGridSpacing(spacing)
86
87 PrintHeader(protein.GetTitle(), ligand.GetTitle())
88
89 CalcBindingEnergy(zap, protein, ligand, cmplx)
90 zap.SetFocusTarget(ligand)
91 CalcBindingEnergy(zap, protein, ligand, cmplx)
92
93 return 0
94
95 if __name__ == "__main__":
96 sys.exit(main(sys.argv))