I have run into a particular problem while trying to automate my input file generation for ORCA calculations. Apparently, I either do not fully understand how SMILES strings work, or I am demanding more than they can really deliver. What I generally do is to use a python script to generate xyz files from SMILES strings to get the initial coordinates of my molecules. This usually works just fine, but metal-organic complexes seem to pose a particular challenge.

As an example, I will show you the SMILES [Zn--]([NH3+])([NH3+])([NH3+])N1C=C(Cl)N=C1, which corresponds to the following molecule:

initial coords xyz from RDKit EmbedMolecule

Right off the bat: There is something weird about the SMILES. I wasn't able to generate the right molecule with the SMILES [Zn--]([NH3+])([NH3+])([NH3+])N1C=C(Cl)N=C1, which would show the correct distribution of charges and numbers of hydrogens on the ammonia ligands. So the only way to get it to work is to use cationic "ammonium" side groups/ligands and compensate their overall charge using negative charges on the Zinc atom.

I've looked for SMILES other metal-organic complexes, which I could use as a reference, to construct a proper SMILES for my mols, but it seems that I am not the only one facing this issue. Take a look at the page for Tetraamminecopper(2+) sulfate on pubchem. The SMILES on there is incorrect and no geometry can be generated on the website. This indicates to me that SMILES strings are, at least regular SMILES strings, not well adapted to depict metal-organic complexes.

This isn't that much of an issue in itself, since I have found a workaround to generate my initial coordinates using RDKit anyways. The real problem arises when I try to do an initial MMFF94 Energy Minimization of the initial coordinates: How is MMFF94 supposed to know which atom types to assign to this molecule and perform it's energy minimization? Short answer: It obviously can't and so the energy minimization isn't performed.

What I have discovered though is that Avogadro can indeed perform an energy minimization using UFF if you load in the xyz file of the initial coordinates. So it seems that Avogadro has some kind of fancy way to assign atom types based on xyz file info alone and perform a proper energy minimization using UFF:

Optimized coords using UFF in Avogadro

Now here is where I got excited: So there is a way to properly preoptimize my mols using UFF and guess what? You can use UFF preoptimization in RDKit:

from rdkit import Chem
from rdkit.Chem import AllChem

# Load SMILES as mol obj
mol = Chem.MolFromSmiles('[Zn--]([NH3+])([NH3+])([NH3+])N1C=C(Cl)N=C1')  
mol = AllChem.AddHs(mol) # make sure to add explicit hydrogens

# Generate 3D coordinates
xyz = Chem.AllChem.MolToXYZBlock(mol) # initial coords

with open('unoptimized.xyz','w') as outf:

# Perform UFF optimization
ff = AllChem.UFFGetMoleculeForceField(mol)

xyz = Chem.AllChem.MolToXYZBlock(mol)
with open('B.xyz','w') as outf:

Here is the problem though: I am still using a fundamentally incorrect SMILES. And while the xyz information does not remember the incorrect charges of the Zinc and the Nitrogens, the RDKit mol object certainly does. And so I cannot use that mol object to preoptimize my mols.

But there could be a way out if I just create a new mol object only from the xyz information from the initial coords and then perform the minimization. It worked in Avogadro, so why not in RDKit, right? Well, unfortunately no...

I added this code fragment to my code after the initial coords generation:

mol = Chem.MolFromXYZBlock(xyz)

which overwrites the mol object with one only generated from xyz information. But when I run the UFF minimization on that mol, the UFF typer of RDKit doesn't let me.

UFFTYPER: Unrecognized charge state for atom: 0

Keep in mind: The charges are no longer provided by the incorrect SMILES, but must be dynamically assigned by RDKit's UFF typer. And since the charge of the Zinc atom is not recognized by the Universal Force Field, the geometry optimization doesn't do anything.

So it would seem that while the RDKit UFF typer is incapable of assigning any meaningful atom types for my mol, the UFF typer of Avogadro can perform some kind of eldritch magic and get it done regardless. Which means a lot of manual optimization for me.

I have tried to switch out the Zinc for other elements (Sn, Si, Pb bc tetragonal) so that the UFF typer would at least do some of the work for me, but no luck.

So my question is this: Is there a way to get RDKit to recognize such mols and preoptimize them anyway? It is clearly possible with Avogadro and it is almost certainly a UFF Typer assignment issue. Maybe there is a better SMILES that I can use? Perhaps the charges of the Zinc can be edited manually somehow? I have some ideas of what I can try, but not how to test them...

EDIT: I messed up the SMILES in the pics. The correct SMILES is: [Zn--]([NH3+])([NH3+])([NH3+])N1C=C(Cl)N=C1


1 Answer 1


Okay, there's a bunch to answer here, but I'll do my best.

Take a look at the page for Tetraamminecopper(2+) sulfate on pubchem. The SMILES on there is incorrect and no geometry can be generated on the website.

Yes, that's because PubChem uses MMFF94 to generate 3D geometries so not all molecules on PubChem have them.

You ask about MMFF94 minimization - not all force fields support metals. MMFF94 supports a few types of isolated ions, but not in molecules.

So it seems that Avogadro has some kind of fancy way to assign atom types based on xyz file info alone and perform a proper energy minimization using UFF:

Avogadro uses Open Babel's implementation for UFF, which I wrote. It's intended to be very forgiving of weird charges and coordination environments. A strict UFF implementation (e.g., RDKit) will insist the Zn3+2 parameters only apply to $\ce{Zn^{2+}}$ ions.

I'm less familiar with feeding metals to RDKit's UFF implementation - your better bet is to ask there.

But you can feed whatever you want into Open Babel as well:

obabel -:'[Zn--]([NH3+])([NH3+])([NH3+])N1C=C(Cl)N=C1' -O zn.xyz --gen3d

Or in Python:

from openbabel import pybel

# assuming a list of SMILES
smiles_list = [ '[Zn--]([NH3+])([NH3+])([NH3+])N1C=C(Cl)N=C1' ]

for smiles in smiles_list:

    mol = pybel.readstring("smi", smiles)

    # make it 3D

    # if possible (e.g., a ligand) use MMFF94
    #  .. otherwise UFF
    ff = pybel._forcefields["mmff94"]
    success = ff.Setup(mol.OBMol)
    if not success:
        ff = pybel._forcefields["uff"]
        success = ff.Setup(mol.OBMol)
        if not success:
           sys.exit("Cannot set up forcefield")

    # optimize a bit and do a quick conformer search
    ff.ConjugateGradients(100, 1.0e-3)
    ff.WeightedRotorSearch(100, 25)
    ff.ConjugateGradients(250, 1.0e-4)

    # update the coordinates

    # pick your favorite format
    # here, I'm suggesting SDF 
    # if you want to retain formal charges
    mol.write("sdf", "test.sdf", overwrite=True)
  • 1
    $\begingroup$ Ah yes...the python lib for openbabel...I have tried and failed to install it multiple times since some wheels refuse to be built. Just tried it again. Same problem I had a year ago...and the year before that. I'd have never gone for RDKit if obabel wasn't so difficult to handle in installation alone. $\endgroup$
    – J.Doe
    Commented Nov 2, 2023 at 13:31
  • $\begingroup$ One other comment is that I'd suggest optimizing beyond UFF, for example GFN-FF or GFN-2 since those methods handle metals much better than UFF. $\endgroup$ Commented Nov 2, 2023 at 13:31
  • $\begingroup$ @J.Doe conda install -c condo-forge openbabel $\endgroup$ Commented Nov 2, 2023 at 13:32
  • $\begingroup$ Also, if you have obabel installed, you can call it from Python with os.system or subprocess $\endgroup$ Commented Nov 2, 2023 at 14:00

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .