3.3 Atom Potentials

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

(a): Induced solvent potential: the potential produced by the polarization of the solvent.
(b): 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.
(c): 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.

 1 #!/usr/bin/env python
 2 #############################################################################
 3 #  Copyright (C) 2006, 2007, 2008 OpenEye Scientific Software, Inc.
 4 #############################################################################
 5
 6 import os, sys
 7 from openeye.oechem import *
 8 from openeye.oezap import *
 9
10 def Output(mol, apot, showAtomTable):
11     print "Title: %s"%mol.GetTitle()
12     if showAtomTable:
13         OEThrow.Info("Atom potentials");
14         OEThrow.Info("Index  Elem    Charge     Potential");
15
16     energy=0.0
17     for atom in mol.GetAtoms():
18         energy += atom.GetPartialCharge()*apot[atom.GetIdx()]
19         if showAtomTable:
20             print "%3d     %2s     %6.3f     %8.3f"%(atom.GetIdx(),
21                    OEGetAtomicSymbol(atom.GetAtomicNum()),
22                    atom.GetPartialCharge(),
23                    apot[atom.GetIdx()])
24
25     print "Sum of {Potential * Charge over all atoms * 0.5} in kT = %f\n" % (0.5*energy)
26
27 def CalcAtomPotentials(itf):
28     mol = OEGraphMol()
29
30     ifs = oemolistream()
31     if not ifs.open(itf.GetString("-in")):
32         OEThrow.Fatal("Unable to open %s for reading" % itf.GetString("-in"))
33
34     OEReadMolecule(ifs,mol)
35     OEAssignBondiVdWRadii(mol)
36
37     if not itf.GetBool("-file_charges"):
38         OEMMFFAtomTypes(mol)
39         OEMMFF94PartialCharges(mol)
40
41     zap = OEZap()
42     zap.SetMolecule(mol)
43     zap.SetInnerDielectric(itf.GetFloat("-epsin"))
44     zap.SetBoundarySpacing(itf.GetFloat("-boundary"))
45     zap.SetGridSpacing(itf.GetFloat("-grid_spacing"))
46
47     showAtomTable = itf.GetBool("-atomtable")
48     calcType = itf.GetString("-calc_type")
49     if calcType=="default":
50         apot = OEFloatArray(mol.GetMaxAtomIdx())
51         zap.CalcAtomPotentials(apot)
52         Output(mol, apot, showAtomTable)
53
54     elif calcType == "solvent_only":
55         apot = OEFloatArray(mol.GetMaxAtomIdx())
56         zap.CalcAtomPotentials(apot)
57
58         apot2 = OEFloatArray(mol.GetMaxAtomIdx())
59         zap.SetOuterDielectric(zap.GetInnerDielectric())
60         zap.CalcAtomPotentials(apot2)
61
62         # find the differences
63         for atom in mol.GetAtoms():
64             idx=atom.GetIdx()
65             apot[idx] -= apot2[idx]
66
67         Output(mol, apot, showAtomTable)
68
69     elif calcType == "remove_self":
70         apot = OEFloatArray(mol.GetMaxAtomIdx())
71         zap.CalcAtomPotentials(apot, True)
72         Output(mol, apot, showAtomTable)
73
74     elif calcType == "coulombic":
75         epsin = itf.GetFloat("-epsin")
76         x = OECoulombicSelfEnergy(mol, epsin)
77         print "Coulombic Assembly Energy"
78         print "  = Sum of {Potential * Charge over all atoms * 0.5} in kT = %f"%x
79         apot = OEFloatArray(mol.GetMaxAtomIdx())
80         OECoulombicAtomPotentials(mol, epsin, apot)
81         Output(mol, apot, showAtomTable)
82
83     return 0
84
85 def SetupInterface(itf, InterfaceData):
86     OEConfigure(itf, InterfaceData)
87     if OECheckHelp(itf, sys.argv):
88         return False
89     if not OEParseCommandLine(itf, sys.argv):
90         return False
91     return True
92
93 def main(InterfaceData):
94     itf=OEInterface()
95     if not SetupInterface(itf, InterfaceData):
96         return 1
97
98     return CalcAtomPotentials(itf)
99
100 InterfaceData="""
101 #zap_atompot interface definition
102
103 !PARAMETER -in
104   !TYPE string
105   !BRIEF Input molecule file.
106   !REQUIRED true
107   !KEYLESS 1
108 !END
109
110 !PARAMETER -file_charges
111   !TYPE bool
112   !DEFAULT false
113   !BRIEF Use partial charges from input file rather than calculating with MMFF.
114 !END
115
116 !PARAMETER -calc_type
117   !TYPE string
118   !DEFAULT default
119   !LEGAL_VALUE default
120   !LEGAL_VALUE solvent_only
121   !LEGAL_VALUE remove_self
122   !LEGAL_VALUE coulombic
123   !LEGAL_VALUE breakdown
124   !BRIEF Choose type of atom potentials to calculate
125 !END
126
127 !PARAMETER -atomtable
128   !TYPE bool
129   !DEFAULT false
130   !BRIEF Output a table of atom potentials
131 !END
132
133 !PARAMETER -epsin
134   !TYPE float
135   !BRIEF Inner dielectric
136   !DEFAULT 1.0
137   !LEGAL_RANGE 0.0 100.0
138 !END
139
140 !PARAMETER -grid_spacing
141   !TYPE float
142   !DEFAULT 0.5
143   !BRIEF Spacing between grid points (Angstroms)
144   !LEGAL_RANGE 0.1 2.0
145 !END
146
147 !PARAMETER -boundary
148   !ALIAS -buffer
149   !TYPE float
150   !DEFAULT 2.0
151   !BRIEF Extra buffer outside extents of molecule.
152   !LEGAL_RANGE 0.1 10.0
153 !END
154 """
155
156 if __name__ == "__main__":
157     sys.exit(main(InterfaceData))

Listing:3.5 Calculating atom potentials

				ZAP Electrostatics at OpenEye Python Programming Manual

				ZAP Electrostatics at OpenEye Python Programming Manual