Skip to content

Instantly share code, notes, and snippets.

@ghutchis
Created August 4, 2023 20:55
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ghutchis/de5359f411efa7cb0803bae654617641 to your computer and use it in GitHub Desktop.
Save ghutchis/de5359f411efa7cb0803bae654617641 to your computer and use it in GitHub Desktop.
Import cclib files to RDKit molecules
#!/usr/bin/env python
import sys
import os
import cclib
from cclib.parser import ccopen
from rdkit import Chem
try:
from rdkit.Chem import rdDetermineBonds
except ImportError:
exit('This script requires RDKit version 2022.09 or later')
# repeat through all the files on the command-line
# we can change this to use the glob module as well
# e.g., find all the files in a set of folders
for argument in sys.argv[1:]:
filename, extension = os.path.splitext(argument)
try:
parser = ccopen(argument)
data = parser.parse()
except:
print("Error importing %s" % argument)
continue
# okay, now let's build the RDKit molecule
mol = Chem.RWMol()
# we'll need space for the 3D coordinates
conformer = Chem.Conformer(len(data.atomnos))
for atom in data.atomnos:
mol.AddAtom(Chem.Atom(int(atom)))
for i in range(len(data.atomcoords[-1])):
conformer.SetAtomPosition(i, data.atomcoords[-1][i])
mol.AddConformer(conformer)
# now the new bond connectivity and bond orders
bonded_mol = Chem.Mol(mol)
rdDetermineBonds.DetermineConnectivity(bonded_mol)
rdDetermineBonds.DetermineBondOrders(bonded_mol, data.charge)
# and finally the stereochemistry
Chem.AssignStereochemistryFrom3D(bonded_mol)
# print the SMILES
smi = Chem.MolToSmiles(bonded_mol)
escape = f"Charge: {data.charge} Spin: {data.mult}"
print(f"From file {argument} I found SMILES {smi}")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment