Hi Tim,
I had a look at the 1DWD example and I detail in the following the
analysis I did as it may be useful to other users.
The conformer generation function has the option to print the
experimental torsion preferences that were used in the generation:
printExpTorsionAngles=True
For 1DWD, it gives the following output:
> ref =
Chem.MolFromSmiles('NC(=[NH2+])c1ccc(C[C@@H](NC(=O)CNS(=O)(=O)c2ccc3ccccc3c2)C(=O)N2CCCCC2)cc1’)
> mol1 =
Chem.MolFromPDBFile(RDConfig.RDBaseDir+'/rdkit/Chem/test_data/1DWD_ligand.pdb')
> mol1 = AllChem.AssignBondOrdersFromTemplate(ref, mol1)
> mol1 = Chem.AddHs(mol1)
> numConfs = 100
> cids = AllChem.EmbedMultipleConfs(mol1, numConfs=numConfs,
useExpTorsionAnglePrefs=True, useBasicKnowledge=True,
printExpTorsionAngles=True)
O=[S:1](=O)[NX3H1:2]!@;-[CX4H2:3][!#1:4]: 0 13 14 15, (17.9, -13.3,
9.2, 4.7, -2.3, -0.9)
O=[C:1][NX3H1:2]!@;-[CX4H1:3][H:4]: 15 17 18 48, (5, 0, 0, 0, 0, 0)
[O:1]=[CX3:2]!@;-[NX3H0:3][!#1:4]: 20 19 31 32, (0, 8, 0, 0, 0, 0)
[O:1]=[CX3:2]!@;-[NX3H1:3][!#1:4]: 16 15 17 18, (100, 0, 0, 0, 0, 0)
[c:1][S:2](=O)(=O)!@;-[NX3H1:3][C:4]: 4 0 13 14, (0, 16, 5, 7, 0, 0)
[aH1:1][c:2]([aH1])!@;-[SX4:3][!#1:4]: 3 4 0 13, (0, 1.5, 0, 0, 0, 0)
N[C:2](=[O:1])!@;-[CH2:3][N:4]: 16 15 14 13, (0, 10, 0, 0, 0, 0)
[!#1:1][CX4:2]!@;-[CX4:3][!#1:4]: 17 18 21 22, (0, 0, 7, 0, 0, 0)
[NH1:1][CX4:2]!@;-[CX3:3]=[O:4]: 17 18 19 20, (0, 2.1, -0.1, -1, -0.4,
-0.1)
[cH1:1][c:2]([cH1])!@;-[CX4H2:3][CX4:4]: 23 22 21 18, (0, 3.6, 0, 0,
0, 0)
[a:1][c:2]!@;-[C:3](=[N:4]): 25 27 28 30, (0, 5, 0, 0, 0, 0)
[*:1][X3,X2:2]=[X3,X2:3][*:4]: 27 28 30 57, (0, 100, 0, 0, 0, 0)
The output is structured as follows: First comes the SMARTS pattern,
then the indices of the four atoms involved and finally the six force
constants K_1, K_2, K_3, K_4, K_5, K_6 for the torsion potential (see
Eq. (2) in JCIM, 55, 2562, 2015).
The position of the non-zero force constant tells you the multiplicity
i of the cosine. E.g. for the torsion between atoms 18 and 21 the
third force constant is 7.0 and all others are zero, thus the torsion
potential for this bond has a multiplicity i = 3 (i.e. three maxima).
The information about the phase shift is not exposed, but we can
change that if people find it useful.
Currently to get the full information about the torsion potential, you
need to search for the SMARTS pattern in
Code/GraphMol/ForceFieldHelpers/CrystalFF/TorsionPreferences.cpp
There you find the line: "[!#1:1][CX4:2]!@;-[CX4:3][!#1:4] 1 0.0 1 0.0
1 7.0 1 0.0 1 0.0 1 0.0\n"
Here we have twelve parameters, always two for a given multiplicity:
first the phase shift (can be 1 or -1) and second the force constant.
For this particular torsion pattern we have a multiplicity of 3, a
phase shift of cos(d) = 1 and a force constant K_3 = 7.0, which
corresponds to three peaks at 60°, 180° and 300° in the range [0°,
360°] (or -60°, 60° and 180° in the range [-180°, 180°]).
The two bonds connecting the benzamidine ring to the rest of the
molecule are between atoms 22/21 and 21/18.
I calculated the signed dihedral angle for these two bonds, here is
the code for the second one:
> angles = []
> for cid in cids:
> conf = mol1.GetConformer(id=cid)
> # torsion of atoms 17 18 21 22
> p1 = conf.GetAtomPosition(17)
> p2 = conf.GetAtomPosition(18)
> p3 = conf.GetAtomPosition(21)
> p4 = conf.GetAtomPosition(22)
> a = Geometry.ComputeSignedDihedralAngle(p1, p2, p3, p4)/math.pi *
180.0
> angles.append(a)
For 100 confs, this gives the following distribution:
So, the torsion potential is sampled as it is supposed to do.
For the other bond between atmos 22/21, we have the following line in
TorsionPreferences.cpp:
"[cH1:1][c:2]([cH1])!@;-[CX4H2:3][CX4:4] 1 0.0 1 3.6 1 0.0 1 0.0 1 0.0
1 0.0\n”
This means, we have a multiplicity of 2, a phase shift cos(d) = 1 and
a force constant of 3.6, which corresponds to two peaks at -120° and
120° in the range [-180°, 180°].
The histogram looks like:
Again, the potential is largely sampled as it’s supposed to be.
Now, let’s have a look at the two bonds connecting atom 18 to the
other two branches.
For the bond between atoms 18/17, we should have a single peak at 0°,
which we get most of the time
The bond between atoms 18/19 has a multi-term potential, which is a
bit more difficult to interpret from the numbers alone.
"[NH1:1][CX4:2]!@;-[CX3:3]=[O:4] 1 0.0 -1 2.1 -1 -0.1 -1 -1.0 1 -0.4 1
-0.1\n"
It has three peaks at -150°, 0° and 150° in the range [-180°, 180°].
The peak at 0° is broad.
The histogram looks like:
For the 100 conformers, the peaks at -150° and 150° are not sampled. I
therefore generated 1000 conformers, but the sampling does not really
improve:
The algorithm of ETKDG takes a randomly generated conformer from
standard distance geometry and minimizes it with the ET and K terms.
For each bond with an exp. torsion term, the torsional angle will be
minimized to the closest minima. In other words, ETKDG relies on
distance geometry for the sampling. Nevertheless, I will have a look
at this particular pattern in the next version of ETKDG to see if we
can improve it.
Best,
Sereina
On 22 Jun 2016, at 10:41, Tim Dudgeon <[email protected]
<mailto:[email protected]>> wrote:
This topic (https://sourceforge.net/p/rdkit/mailman/message/35173301/)
discussed using conformer generation as input into Open3DAlign.
One thing I noticed is that the conformer generation (using the
useExpTorsionAnglePrefs=True and
useBasicKnowledge=True options) does not generate conformers that align
well for this example. The input is based on the 1DWD_ligand.pdb
structure in the RDKit distro. What I find is the the rotation of the
benzamidine ring is never in the right place for alignment (the other
two ring systems align well). This is when generating up to 10,000
conformers.
Does this suggest that the conformer generation does not sample
conformational space very effectively? Are there options for
improving this?
Thanks
Tim
------------------------------------------------------------------------------
Attend Shape: An AT&T Tech Expo July 15-16. Meet us at AT&T Park in San
Francisco, CA to explore cutting-edge tech and listen to tech luminaries
present their vision of the future. This family event has something for
everyone, including kids. Get more information and register today.
http://sdm.link/attshape
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
