Created
August 4, 2023 20:55
-
-
Save ghutchis/de5359f411efa7cb0803bae654617641 to your computer and use it in GitHub Desktop.
Import cclib files to RDKit molecules
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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