The potentials at any charge, as calculated by PB, can be decomposed into three forms:
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 }