create molecules using RDKit

molecule_generator.py

from rdkit import Chem
from rdkit.Chem import AllChem

parent_smiles = 'c1ccccc1'
rxn_smarts_list = ['[cX3;H1:1]>>[*:1]F','[cX3;H1:1]>>[*:1]O']
molecule_substitutions = {}
molecules = []
mol = Chem.MolFromSmiles(parent_smiles)
molecule_substitutions[0] = [mol]
molecules.append(mol)

substitutions = 6

for i in range(1,substitutions+1):
    molecule_substitutions[i] = []
    smiles_list = []
    for mol in molecule_substitutions[i-1]:
        for rxn_smarts in rxn_smarts_list:
            rxn = AllChem.ReactionFromSmarts(rxn_smarts)
            new_mols = rxn.RunReactants((mol,))

            for new_mol in new_mols:
                new_smiles = Chem.MolToSmiles(new_mol[0],isomericSmiles=True)
                if new_smiles not in smiles_list:
                    smiles_list.append(new_smiles)
                    molecule_substitutions[i].append(Chem.MolFromSmiles(new_smiles))
                    
    molecules += molecule_substitutions[i]

# Write 3D coordinates in sdf format
folder = "/Users/jan/Desktop/"
for i, mol in enumerate(molecules):
    mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol)
    file = folder+"comp"+str(i)+".sdf"
    writer = Chem.SDWriter(file)
    writer.write(mol)