# Computing force fields using Open Babel

I am trying to use Open Babel (OB) to compute universal force field (UFF) given an xyz file. I read each xyz file and turn it into a pybel molecule. Here is the code I have so far

def xyz_to_smiles(fname: str) -> str:
mol = next(pybel.readfile("xyz", fname))
smi = mol.write(format="smi")

return smi.split()[0].strip(), mol

def get_qm9_OB(target_index=np.arange(0, 133884)):

smiles = []
molecule = []
for i in target_index:
mol_idx = f"{i+1:06d}"

tmp_smiles, tmp_molecule = xyz_to_smiles('../dsgdb9nsd.xyz/dsgdb9nsd_' + mol_idx + '.xyz')

smiles.append(tmp_smiles)
molecule.append(tmp_molecule)

return molecule, smiles

molecules, smiles = get_qm9_OB(target_index=np.arange(100,200))

for i in range(len(smiles)):
mol = molecules[i].OBMol

ff = ob.OBForceField.FindForceField('UFF')
ff.Setup(mol)

ff.GetAtomTypes(mol)
ff.GetCoordinates(mol)
for at in ob.OBMolAtomIter(mol):
type = at.GetData("FFAtomType") # an OBPairData object
print(type.GetValue())


Now, I am trying to compute bond stretching interactions. My questions are:

1. How can I get the UFF constants for bond stretching? (IOW, is there an OB equivalent for RDKit function GetUFFBondStretchParams?)
2. Is there a function to calculate bond stretching interactions for each bond? (IOW, is there an OB implementation of either harmonic oscillator or the Morse function?)
3. How can I get coordinates from OBAtom object? I found GetCoordinate() but the output is

<Swig Object of type 'double *' at 0x1a1eb4a9f0>


Alternatively, GetVector() outputs

<openbabel.vector3; proxy of <Swig Object of type 'OpenBabel::vector3 *' at 0x1a1fbddbd0> >


How can I get the actual values?

P.S. I posted the problem on OB mailing list, with no answers.

• I honestly don't understand your question, and your code seems incomplete; I don't see how you are actually reading in an xyz file. I'm not sure how well this question is actually suited here, as it seems to be rather a programming issue. Maybe it is more on topic at another site, e.g. scientific computing it even stack overflow. – Martin - マーチン Oct 26 '19 at 12:13
• @Martin-マーチン Oh, I forgot to include xyz_to_smiles. Now you can see how it reads the file. Let me know what part you don't understand so that I can clarify. – Blade Oct 26 '19 at 18:15
• I did not see your message - not sure why others didn't respond, though. – Geoff Hutchison Oct 30 '19 at 19:51

1. While people have asked, at the moment, there's no API call to get parameters for various force fields. It's possible to hack the C++ code, of course. (It hasn't been clear to me what such an API call should look like - feel free to suggest something.)

2. There isn't an implementation because a harmonic oscillator and Morse potential are fairly easy to implement yourself. (Incidentally, despite mention of Morse in the UFF paper, the parameterization involves a harmonic oscillator, not a Morse potential for bonds.)

3. IIRC there's a SWIG convenience function to cast from vector3 to Python vectors, but you see that there are multiple ways to do this.

1. Coordinates can be accessed using

[OBAtom.GetX(),OBAtom.GetY(),OBAtom.GetZ()]


look at this Github discussion.

If you have a better answer please don't hesitate to answer the question.