1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ – Martin - マーチン Oct 26 '19 at 12:13
  • 1
    $\begingroup$ @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. $\endgroup$ – Blade Oct 26 '19 at 18:15
  • 1
    $\begingroup$ I did not see your message - not sure why others didn't respond, though. $\endgroup$ – Geoff Hutchison Oct 30 '19 at 19:51
2
$\begingroup$
  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.

$\endgroup$
3
$\begingroup$
  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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.