Computing force fields using Open Babel

Question

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.

Answer 1

3. 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.

Answer 2

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.