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 }