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 #!/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))

Listing:3.10 Calculating Solvation Forces