Using OpenBabel's gen3d in a Python notebook

Question

I want to preface this by saying I am really new to computational chemistry and this kind of programming in general, but I am having real trouble figuring out how to use OpenBabel's gen3d function in a Python notebook. I want to use it to make a single force field optimized conformer from a SMILES string but I don't really know how to call it on my SMILES string. Any help or links would be greatly appreciated! Thank you!

Answer 1

Here is a simple python recipe in addition to the CLI approach mentioned in the other answer, for more details see the docs:

from openbabel import pybel as pb # pybel is syntax sugar for openbabel

smiles = '[O-]C(=O)[C@H]([NH3+])Cc1ccc(O)cc1' # your smiles
mol = pb.readstring('smiles', smiles) # first argument is a string format

mol.addh() # add Hs for 3D
mol.make3D() # generate 3D coordinates

However, my experience shows that open babel is not very reliable for 3D generation. I recommend looking through the RDKit library, which is in my opinion more powerful and handy tool if you need to work with molecules:

from rdkit import Chem # main tools
from rdkit.Chem import AllChem # additional tools, including 3D

smiles = '[O-]C(=O)[C@H]([NH3+])Cc1ccc(O)cc1' # your smiles
mol = Chem.MolFromSmiles(smiles) # initialize molecule

mol = Chem.AddHs(mol) # adding explicit Hs for 3D generation
cid = AllChem.EmbedMolecule(mol) # returns the id of the generated conformer,
                                 # and -1 if no conformers were generated

AllChem.MMFFOptimizeMolecule(mol) # optimize molecule with MMFF94
Chem.MolToMolFile(mol, 'some_amino_acid.mol') # save to MOL-file

For multiple conformers use AllChem.EmbedMultipleConformers and AllChem.MMFFOptimizeMoleculeConfs, for more details see the docs, they are written really nicely in contrast to OpenBabel.

Answer 2

From the openbabel documentation :

Open Babel will not generate coordinates unless asked, so while a conversion from SMILES to SDF will generate a valid SDF file, the resulting file will not contain coordinates. To generate coordinates, use either the --gen3d or the --gen2d option.

The command line option --conformer allows performing conformer searches using a range of different algorithms and options.

To generate 3D coordinates and conformer searching you can do the following:

In a text editor add the SMILES, one per line, and save the file with .smi extension.

For example add the following in a .smi file:

CC
CCN
CO

After that you can use the --gen3d and --conformer as follows:

obabel infile.smi -O outfile.sdf --gen3d --conformer --nconf 20 --weighted

nconf is the number of conformers.

The output file will look like this:

OpenBabel03192322333D

  8  7  0  0  0  0  0  0  0  0999 V2000
    0.9292   -0.0246   -0.0664 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.4251   -1.1907    0.7537 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.5744    0.9217    0.3526 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.5744   -0.0990   -1.0987 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.0230   -0.0076   -0.0784 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.7799   -1.1164    1.7859 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6687   -1.2078    0.7657 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.7799   -2.1371    0.3346 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1  0  0  0  0
  1  5  1  0  0  0  0
  2  6  1  0  0  0  0
  2  7  1  0  0  0  0
  2  8  1  0  0  0  0
M  END
$$$$

 OpenBabel03192322333D

 10  9  0  0  0  0  0  0  0  0999 V2000
    1.0339    0.0918    0.0081 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.5341    0.6110   -1.3285 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.0092   -0.2403   -2.4112 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.6799   -0.9278    0.1942 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6696    0.7287    0.8207 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.1284    0.0873    0.0482 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.5606    0.6263   -1.3329 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.8849    1.6362   -1.4846 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.7017    0.1435   -3.3037 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.0281   -0.2323   -2.4323 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  4  1  0  0  0  0
  1  5  1  0  0  0  0
  1  6  1  0  0  0  0
  2  3  1  0  0  0  0
  2  7  1  0  0  0  0
  2  8  1  0  0  0  0
  3  9  1  0  0  0  0
  3 10  1  0  0  0  0
M  END
$$$$

 OpenBabel03192322333D

  6  5  0  0  0  0  0  0  0  0999 V2000
    1.0511    0.0027    0.0995 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.5766    0.8295    1.1462 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.6845    0.3908   -0.8536 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6815   -1.0157    0.2421 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.1439    0.0081    0.0964 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.9177    0.4607    1.9788 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1  0  0  0  0
  1  5  1  0  0  0  0
  2  6  1  0  0  0  0
M  END
$$$$

Answer 3

The command-line --gen3d option can do a variety of things, including coordinate generation, some force field cleanup, and a quick conformer search.

As mentioned in another answer, Pybel has make3D() which generates coordinates followed by a quick cleanup using MMFF94 or UFF (if MMFF94 doesn't work).

On the other hand, the other answer suggests (correctly) to perform a conformer search, but not how to do that in Python...

from openbabel import pybel as pb # pybel is syntax sugar for openbabel
from openbabel import openbabel as ob

smiles = '[O-]C(=O)[C@H]([NH3+])Cc1ccc(O)cc1' # your smiles
mol = pb.readstring('smiles', smiles) # first argument is a string format

mol.make3D() # generate 3D coordinates - handles adding H atoms
# also handles a quick forcefield cleanup

# now perform conformer searching
ff = pb._forcefields["mmff94"]
success = ff.Setup(mol.OBMol)
if not success:
    ff = pb._forcefields["uff"]
    success = ff.Setup(mol.OBMol)
    if not success:
        sys.exit("Cannot set up forcefield")

# feel free to tweak these to your balance of time / quality
ff.ConjugateGradients(100, 1.0e-3)
ff.FastRotorSearch(true) # permute central bonds
ff.WeightedRotorSearch(100, 25) # 100 cycles, each with 25 forcefield ops
# final optimization
ff.ConjugateGradients(250, 1.0e-4)
# update the coordinates
ff.GetCoordinates(mol.OBMol)

# or pick your favorite format or filename
mol.write("sdf", "optimized.sdf", overwrite=True)