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:
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:
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
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: