Hi David,
Thanks for the clarification. One more computationally-expensive and naive
approach would be to generate multiple conformers and check them against a
known reference structure.
Here's a quick take on the approach that will be a bit slower. It assumes that
mol is a probe molecule that has no conformers, ref_mol is a mol that has the
desired geometry, and that you already know what the core scaffold is:
mol=
Chem.MolFromSmiles('[17OH]C1=C2C=C(N)C=C1CC3=C([17O-])C(CC4=C([17OH])C(CC5=C(N)C(C2)=CC(S(=O)([O-])=O)=C5)=CC(S(=O)([O-])=O)=C4)=CC(N)=C3')
add_conformer_match(mol, ref_mol,
'C1(C2)=CC(CC3=CC(CC4=CC=CC(CC5=CC=CC2=C5)=C4)=CC=C3)=CC=C1')
def add_conformer_match(probe_mol, # assume no conformers yet
reference_mol, # assume the reference has a single
correct conf
substructure, # the conserved scaffold
threshold=0.5, # define an rmsd threshold
num_confs=50,
num_threads=4):
# create a consistent scaffold to match atom_ids
scaffold = Chem.MolFromSmiles(substructure)
probe_indices = probe_mol.GetSubstructMatch(scaffold)
ref_indices = reference_mol.GetSubstructMatch(scaffold)
atom_map = list(zip(probe_indices, ref_indices))
# generate diverse conformers for new molecule
probe_mol_h = Chem.AddHs(probe_mol)
conf_ids = AllChem.EmbedMultipleConfs(probe_mol_h,
numConfs=num_confs,
numThreads=num_threads,
pruneRmsThresh=0.5)
probe_mol_confs = Chem.RemoveHs(probe_mol_h)
# iterate over conformers to see if they match
rms_list = []
for conf_id in conf_ids:
rmsd = AllChem.AlignMol(probe_mol_confs,
reference_mol,
prbCid=conf_id,
atomMap=atom_map)
rms_list.append(rmsd)
# check to see if you find a reasonable match
if min(rms_list) > threshold:
print(f'No conformer found within RMS threshold of {threshold}')
return None
else:
# add the lowest_rms conformer to the original mol
lowest_rms = rms_list.index(min(rms_list))
probe_mol.AddConformer(probe_mol_confs.GetConformer(lowest_rms))
return probe_mol # return the original object with added conformer
Kangway
________________________________
From: David Turnbull <[email protected]>
Sent: Thursday, July 23, 2020 7:58 AM
To: [email protected] <[email protected]>;
Chuang, Kangway <[email protected]>
Subject: Re: Conformer generation
That is the structure I want, however I found that it doesn't give that
structure every time (sometimes it inverts the rings).
Get Outlook for
Android<https://urldefense.proofpoint.com/v2/url?u=https-3A__aka.ms_ghei36&d=DwMGaQ&c=iORugZls2LlYyCAZRB3XLg&r=Z0E5F87lf3GPcsIl1f2OYQw4iwqHJfffu3dwlNgH2Zs&m=zKUYd_Jfhr2WYR_OBhn5whmIq1pRnUMMWKZXet-DsRg&s=_KGfUMbYuqPwR3m2g7VO79jxh9F-XP_TRtO9itOLZ40&e=>
________________________________
From: Chuang, Kangway <[email protected]>
Sent: Thursday, July 23, 2020 8:55:01 AM
To: David Turnbull <[email protected]>;
[email protected] <[email protected]>
Subject: Re: Conformer generation
[△EXTERNAL]
Hi David,
Do you have a specific example of the bowl conformation you're looking for
(e.g. do you have an image of the desired conformation vs what you are seeing)?
Running your current code I get the following conformer (shown in two different
views).
Kangway
________________________________
From: David Turnbull <[email protected]>
Sent: Thursday, July 23, 2020 6:49 AM
To: [email protected] <[email protected]>
Subject: [Rdkit-discuss] Conformer generation
Hi all,
I am trying to generate structures of calixarenes in a set shape, I am trying
to use constrain distances but struggling. The lowest energy conformer is not
what I want as I want it in the bowl shape. I labelled 4 oxygens with O17 (for
finding purposes) and set the distance of them as that should set the geometry
but it isn’t working. Any help would be great. Code below.
mol=
Chem.MolFromSmiles('[17OH]C1=C2C=C([Y])C=C1CC3=C([17O-])C(CC4=C([17OH])C(CC5=C([17OH])C(C2)=CC(S(=O)([O-])=O)=C5)=CC(S(=O)([O-])=O)=C4)=CC([Y])=C3')
y= Chem.MolFromSmiles('OC1=CC=C([Y])C=C1')
rxn = AllChem.ReactionFromSmarts("[Y][*:1].[Y][*:2]>>[*:1][*:2]")
results=rxn.RunReactants([mol,y])
for products in results:
for mol in products:
m2=mol
results2=rxn.RunReactants([m2,y])
for products in results2:
for mol in products:
m3=mol
m4=Chem.MolToSmiles(m3)
x=Chem.MolFromSmiles(m4)
index={}
k=0
for atom in x.GetAtoms():
if atom.GetSymbol() == 'O':
if atom.GetIsotope() == 17:
index[k]=atom.GetIdx()
k+=1
m5=Chem.AddHs(x)
AllChem.EmbedMolecule(m5, useRandomCoords=True)
ff=AllChem.UFFGetMoleculeForceField(m5)
ff.UFFAddDistanceConstraint( index[0], index[1], False, 3.5, 4.5, 99.0)
ff.UFFAddDistanceConstraint( index[0], index[2], False, 3.5, 4.5, 99.0)
ff.UFFAddDistanceConstraint( index[0], index[3], False, 3.5, 4.5, 99.0)
m6=ff.Minimize(maxIts=50)
David
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
