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 #!/usr/bin/env python
2 #############################################################################
3 # Copyright (C) 2006, 2007, 2008 OpenEye Scientific Software, Inc.
4 #############################################################################
5 import os, sys
6 from openeye.oechem import *
7 from openeye.oezap import *
8
9 def PrintHeader(title):
10 OEThrow.Info("\nForce Components for %s, in kT/Angstrom" % title)
11 OEThrow.Info("%6s %8s %8s %8s %8s" % ("Index", "Element", "-dE/dx",
12 "-dE/dy", "-dE/dz"))
13
14 def PrintForces(mol, forces):
15 for atom in mol.GetAtoms():
16 OEThrow.Info("%6d %8s %8.2f %8.2f %8.2f" % (atom.GetIdx(),
17 OEGetAtomicSymbol(atom.GetAtomicNum()), forces[atom.GetIdx()],
18 forces[atom.GetIdx()+1], forces[atom.GetIdx()+2]))
19
20 def main(argv = [__name__]):
21 if len(argv) != 2:
22 OEThrow.Usage("%s <molfile>" % argv[0])
23
24 epsin = 1.0
25 zap = OEZap()
26 zap.SetInnerDielectric(epsin);
27
28 mol = OEGraphMol()
29 ifs = oemolistream()
30 if not ifs.open(argv[1]):
31 OEThrow.Fatal("Unable to open %s for reading" % argv[1])
32 while OEReadMolecule(ifs, mol):
33 PrintHeader(mol.GetTitle())
34 forces = OEFloatArray(mol.GetMaxAtomIdx()*3)
35 OEAssignBondiVdWRadii(mol)
36 OEMMFFAtomTypes(mol)
37 OEMMFF94PartialCharges(mol)
38 zap.SetMolecule(mol)
39 zap.CalcForces(forces)
40 PrintForces(mol, forces)
41
42 return 0
43
44 if __name__ == "__main__":
45 sys.exit(main(sys.argv))