3.7 Gradients/Forces

The solvent forces acting on the atoms in a molecule may be calculated using the CalcForces method. The components of the forces are set to a float array of length 3N, where N is the number of atoms. The components of the gradient are equal in magnitude to the force, but have the opposite sign.

 1 /*******************************************************************
 2     Copyright (C) 2006, 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 using namespace OEPB;
12 using namespace OEChem;
13 using namespace OESystem;
14 using namespace OEPlatform;
15
16 void PrintHeader(const char *title)
17 {
18   OEThrow.Info("\nForce Components for %s, in kT/Angstrom",title);
19   OEThrow.Info("%6s %8s %8s %8s %8s", "Index", "Element", "-dE/dx",
20                "-dE/dy", "-dE/dz");
21 }
22
23 void PrintForces(const OEMolBase &mol, const float *forces)
24 {
25   OEIter<OEAtomBase> atom;
26   for(atom=mol.GetAtoms(); atom; ++atom)
27   {
28     OEThrow.Info("%6d %8s %8.2f %8.2f %8.2f",
29            atom->GetIdx(), OEGetAtomicSymbol(atom->GetAtomicNum()),
30            forces[atom->GetIdx()], forces[atom->GetIdx()+1],
31            forces[atom->GetIdx()+2]);
32   }
33 }
34
35 int main(int argc, char *argv[])
36 {
37   if (argc!=2)
38     OEThrow.Usage("%s <molfile>", argv[0]);
39
40   oemolistream ifs;
41   if (!ifs.open(argv[1]))
42     OEThrow.Fatal("Could not open %s for reading", argv[1]);
43
44   float epsin = 1.0f;
45   OEGraphMol mol;
46   OEZap zap;
47   zap.SetInnerDielectric(epsin);
48
49   while (OEReadMolecule(ifs,mol))
50   {
51     PrintHeader(mol.GetTitle());
52     float *forces = new float[mol.GetMaxAtomIdx()*3];
53     OEAssignBondiVdWRadii(mol);
54     OEMMFFAtomTypes(mol);
55     OEMMFF94PartialCharges(mol);
56     zap.SetMolecule(mol);
57     zap.CalcForces(forces);
58     PrintForces(mol, forces);
59     delete [] forces;
60   }
61   return 0;
62 }

Listing:3.10 Calculating Solvation Forces