 |
 |
 | ZAP
Electrostatics at OpenEye C++ Programming Manual | ||||
 |
 |
 |
 |
||||
|
|
| ZAP
Electrostatics at OpenEye C++ Programming Manual | ||||
|
|
|
|
||||
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) 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 #include "zap_atompot.itf"
12
13 using namespace OEPB;
14 using namespace OEChem;
15 using namespace OESystem;
16 using namespace OEPlatform;
17
18 static void Output(const OEMolBase &mol, float *apot, bool table)
19 {
20 OEThrow.Info("Title: %s", mol.GetTitle());
21 if (table)
22 {
23 OEThrow.Info("Atom potentials");
24 OEThrow.Info("Index Elem Charge Potential");
25 }
26 double energy=0.0;
27 for (OEIter<OEAtomBase> atom=mol.GetAtoms();atom;++atom)
28 {
29 energy += atom->GetPartialCharge()*apot[atom->GetIdx()];
30 if (table)
31 OEThrow.Info("%3d %2s %6.3f %8.3f",
32 atom->GetIdx(),OEGetAtomicSymbol(atom->GetAtomicNum()),
33 atom->GetPartialCharge(),apot[atom->GetIdx()]);
34 }
35 OEThrow.Info("Sum of {Potential * Charge over all atoms * 0.5} in kT = %f\n",
36 0.5*energy);
37 }
38
39 static bool SetupInterface(int argc, char *argv[], OEInterface &itf)
40 {
41 OEConfigure(itf, InterfaceData);
42 if(OECheckHelp(itf, argc, argv))
43 return false;
44 if (!OEParseCommandLine(itf, argc, argv))
45 return false;
46 std::string filename=itf.Get<std::string>("-in");
47 if(!OEIsReadable(OEGetFileType(OEGetFileExtension(filename.c_str()))))
48 {
49 OEThrow.Warning("%s is not a readable input file", filename.c_str());
50 return false;
51 }
52 return true;
53 }
54
55 int main(int argc, char *argv[])
56 {
57 OEInterface itf;
58 if (!SetupInterface(argc, argv, itf))
59 return 1;
60
61 OEGraphMol mol;
62
63 oemolistream ifs;
64 if (!ifs.open(itf.Get<std::string>("-in")))
65 OEThrow.Fatal("Could not open %s for reading", argv[1]);
66
67 OEReadMolecule(ifs,mol);
68 OEAssignBondiVdWRadii(mol);
69
70 if (!itf.Get<bool>("-file_charges"))
71 {
72 OEMMFFAtomTypes(mol);
73 OEMMFF94PartialCharges(mol);
74 }
75
76 OEZap zap;
77 zap.SetMolecule(mol);
78 zap.SetInnerDielectric(itf.Get<float>("-epsin"));
79 zap.SetBoundarySpacing(itf.Get<float>("-boundary"));
80 zap.SetGridSpacing(itf.Get<float>("-grid_spacing"));
81
82 std::string calcType = itf.Get<std::string>("-calc_type");
83 if (calcType=="default")
84 {
85 float *apot = new float[mol.GetMaxAtomIdx()];
86 zap.CalcAtomPotentials(apot);
87 Output(mol,apot,itf.Get<bool>("-atomtable"));
88 delete [] apot;
89 }
90 else if (calcType=="solvent_only")
91 {
92 float *apot = new float[mol.GetMaxAtomIdx()];
93 zap.CalcAtomPotentials(apot);
94
95 float *apot2 = new float[mol.GetMaxAtomIdx()];
96 zap.SetOuterDielectric(zap.GetInnerDielectric());
97 zap.CalcAtomPotentials(apot2);
98
99 // find the difference
100 for (OEIter<OEAtomBase> atom=mol.GetAtoms();atom;++atom)
101 apot[atom->GetIdx()] -= apot2[atom->GetIdx()];
102
103 Output(mol,apot,itf.Get<bool>("-atomtable"));
104
105 delete [] apot;
106 delete [] apot2;
107 }
108 else if (calcType=="remove_self")
109 {
110 float *apot = new float[mol.GetMaxAtomIdx()];
111 zap.CalcAtomPotentials(apot,true);
112 Output(mol,apot,itf.Get<bool>("-atomtable"));
113 delete [] apot;
114 }
115 else if (calcType=="coulombic")
116 {
117 float epsin = itf.Get<float>("-epsin");
118 float x = OECoulombicSelfEnergy(mol, epsin);
119 OEThrow.Info("Coulombic Assembly Energy ");
120 OEThrow.Info(" = Sum of {Potential * Charge over all atoms * 0.5} in kT "
121 "= %f",x);
122
123 float *apot = new float[mol.GetMaxAtomIdx()];
124 OECoulombicAtomPotentials(mol, epsin, apot);
125
126 Output(mol,apot,itf.Get<bool>("-atomtable"));
127 delete [] apot;
128 }
129
130 return 0;
131 }
|
|
| ZAP
Electrostatics at OpenEye C++ Programming Manual | ||||
|
|
|
|
||||