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:

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 /*******************************************************************************
 2  * Copyright (C) 2006, 2008 OpenEye Scientific Software, Inc.
 3  ******************************************************************************/
 4 package openeye.examples.oezap;
 5
 6 import openeye.oechem.*;
 7 import openeye.oezap.*;
 8
 9 public class ZapAtoms1 {
10   private static String InterfaceString = "!PARAMETER -in \n"
11       + "  !TYPE string\n"
12       + "  !BRIEF Input molecule file\n"
13       + "  !REQUIRED true\n"
14       + "!END\n"
15       + "!PARAMETER -epsin\n"
16       + "  !TYPE float\n"
17       + "  !BRIEF Inner dielectric\n"
18       + "  !DEFAULT 1.0\n"
19       + "  !LEGAL_RANGE 0.0 100.0\n"
20       + "!END\n"
21       + "!PARAMETER -epsout\n"
22       + "  !TYPE float\n"
23       + "  !BRIEF Outer dielectric\n"
24       + "  !DEFAULT 80.0\n"
25       + "  !LEGAL_RANGE 0.0 100.0\n"
26       + "!END\n"
27       + "!PARAMETER -grid_spacing\n"
28       + "  !TYPE float\n"
29       + "  !DEFAULT 0.5\n"
30       + "  !BRIEF Spacing between grid points (Angstroms)\n"
31       + "  !LEGAL_RANGE 0.1 2.0\n"
32       + "!END\n"
33       + "!PARAMETER -boundary\n"
34       + "  !ALIAS bubber\n"
35       + "  !TYPE float\n"
36       + "  !DEFAULT 2.0\n"
37       + "  !BRIEF Extra buffer outside extents of molecule.\n"
38       + "  !LEGAL_RANGE 0.1 10.0\n"
39       + "!END\n"
40       + "!PARAMETER -file_charges\n"
41       + "!TYPE bool\n"
42       + "!DEFAULT false\n"
43       + "!BRIEF Use partial charges from input file rather than calculating with MMFF.\n"
44       + "!END\n" + "!PARAMETER -calc_type\n" + "  !TYPE string\n"
45       + "  !DEFAULT default\n" + "  !LEGAL_VALUE default\n"
46       + "  !LEGAL_VALUE solvent_only\n" + "  !LEGAL_VALUE remove_self\n"
47       + "  !LEGAL_VALUE coulombic\n" + "  !LEGAL_VALUE breakdown\n"
48       + "  !BRIEF Choose type of atom potentials to calculate\n" + "!END\n"
49       + "!PARAMETER -atomtable\n" + "  !TYPE bool\n" + "  !DEFAULT false\n"
50       + "  !BRIEF Output a table of atom potentials\n" + "!END\n";
51
52   private static void Output(OEMolBase mol, OEFloatArray apot, boolean atomtable) {
53     System.out.println("Title: " + mol.GetTitle());
54     if (atomtable) {
55       System.out.println("Atom potentials");
56       System.out.println("Index\t Elem\t Charge\t Potential");
57     }
58     double energy = 0.0f;
59     OEAtomBaseIter aiter = mol.GetAtoms();
60     while (aiter.hasNext()) {
61       OEAtomBase atom = aiter.next();
62       energy += atom.GetPartialCharge() * apot.getItem(atom.GetIdx());
63       if (atomtable) {
64         System.out.println(atom.GetIdx() + "\t"
65             + oechem.OEGetAtomicSymbol(atom.GetAtomicNum()) + "\t"
66             + atom.GetPartialCharge() + "\t" + apot.getItem(atom.GetIdx()));
67       }
68     }
69     System.out.println("Sum of (Potential * Charge * 0.5) in kT = " + energy
70         * 0.5f);
71   }
72
73   public static void main(String argv[]) {
74     OEInterface itf = new OEInterface();
75     if (!oechem.OEConfigure(itf, InterfaceString))
76       System.exit(1);
77
78     if (oechem.OECheckHelp(itf, argv))
79       System.exit(0);
80
81     oechem.OEParseCommandLine(itf, argv);
82     OEGraphMol mol = new OEGraphMol();
83     oemolistream ifs = new oemolistream();
84     if (!ifs.open(itf.GetString("-in"))) {
85       System.err.println("Unable to open for reading: " + argv[0]);
86       System.exit(1);
87     }
88     oechem.OEReadMolecule(ifs, mol);
89     oechem.OEAssignBondiVdWRadii(mol);
90     if (!itf.GetBool("-file_charges")) {
91       oechem.OEMMFFAtomTypes(mol);
92       oechem.OEMMFF94PartialCharges(mol);
93     }
94     ifs.close();
95
96     OEZap zap = new OEZap();
97     zap.SetMolecule(mol);
98     zap.SetInnerDielectric(itf.GetFloat("-epsin"));
99     zap.SetBoundarySpacing(itf.GetFloat("-boundary"));
100     zap.SetGridSpacing(itf.GetFloat("-grid_spacing"));
101
102     String calcType = itf.GetString("-calc_type");
103     if (calcType.equals("default")) {
104       OEFloatArray apot = new OEFloatArray(mol.GetMaxAtomIdx());
105       zap.CalcAtomPotentials(apot.ptrCast());
106       System.out.println("Before Output");
107       Output(mol, apot, itf.GetBool("-atomtable"));
108     } else if (calcType.equals("solvent_only")) {
109       OEFloatArray apot = new OEFloatArray(mol.GetMaxAtomIdx());
110       zap.CalcAtomPotentials(apot.ptrCast());
111       OEFloatArray apot2 = new OEFloatArray(mol.GetMaxAtomIdx());
112       zap.SetOuterDielectric(zap.GetInnerDielectric());
113       zap.CalcAtomPotentials(apot2.ptrCast());
114       OEAtomBaseIter aiter = mol.GetAtoms();
115       float potdiff;
116       while (aiter.hasNext()) {
117         OEAtomBase atom = aiter.next();
118         potdiff = apot.getItem(atom.GetIdx()) - apot2.getItem(atom.GetIdx());
119         apot.setItem(atom.GetIdx(), potdiff);
120       }
121       Output(mol, apot, itf.GetBool("-atomtable"));
122     } else if (calcType.equals("remove_self")) {
123       OEFloatArray apot = new OEFloatArray(mol.GetMaxAtomIdx());
124       zap.CalcAtomPotentials(apot.ptrCast());
125       Output(mol, apot, itf.GetBool("-atomtable"));
126     } else if (calcType.equals("coulombic")) {
127       float x = oezaplib.OECoulombicSelfEnergy(mol, itf.GetFloat("-epsin"));
128       System.out.println("Coulombic Assembly Energy ");
129       System.out
130           .println("Sum of {Potential * Charge over all atoms * 0.5} in kT = "
131               + x);
132     }
133   }
134 }

Listing:3.4 Calculating atom potentials