[Rdkit-discuss] Alignment of alkane chains

[Patel, Lara Anne via Rdkit-discuss]

Dear Rdkit Community,

I am currently trying to figure out how to get two alkane chains of the exact 
same length to align a top one another.  Hypothetically, it should be possible 
to consistently get the molecules to align onto one another each time.  Thus 
far my code seems to be giving me a number of faulty alignments and 
occasionally has given me conformations that do not make sense.  I have tried 
several different variations of the code using the Chem.rdMolAlign.GetO3A() and 
Chem.rdMolAlign.AlignMol()functions that have consistently gotten the same bad 
alignments.  An example of this is:

[cid:image001.png@01D66904.D2C9BB30]

The red indicates the configuration of the reference molecule that all other 
hexane chains (cyan) are aligned onto.  Any help that could be given on how I 
might avoid the bad alignments would be much appreciated.
A sample of the code to reproduce the alignment of hexane onto hexane is:

import sys, os
from rdkit.Chem import rdMolAlign, rdForceFieldHelpers, AllChem
from rdkit import Chem
import numpy as np

if __name__ == '__main__' :

    ## Smile string for the input file.

    smile_Ref = "CCCCCC"

    smile_Targ = "CCCCCC"

    ## Opening output file.

    FP_SCORE = "RMD-Score.out"

    F_SCORE  = open(FP_SCORE,"w")

    ## Generating the reference molecule from the smile strings,

    ## adding Hs for proper embedding in 3D space, and embedding

    ## the molecule in 3D with a randomized orientation.

    MOL_Ref  = Chem.MolFromSmiles(smile_Ref)

    MOL_Ref = Chem.AddHs(MOL_Ref)

    AllChem.EmbedMolecule(MOL_Ref , useRandomCoords=True, randomSeed=-1)

    rdForceFieldHelpers.UFFOptimizeMolecule(MOL_Ref)

    Chem.MolToPDBFile(MOL_Ref,"newMol-round{0:0>4}.pdb".format(0))

    ## We do the alignment several times to see if poor alignments arise.

    for i in range(1,100):

        ## Generating the target molecule from the smile strings,

        ## adding Hs for proper embedding in 3D space, and embedding

        ## the molecule in 3D with a randomized orientation.

        MOL_Targ = Chem.MolFromSmiles(smile_Targ)

        MOL_Targ = Chem.AddHs(MOL_Targ)

        AllChem.EmbedMolecule(MOL_Targ, useRandomCoords=True, randomSeed=-1)

        rdForceFieldHelpers.UFFOptimizeMolecule(MOL_Targ)

        ############################################################

        ############################################################

        ## Aligning the Target molecule onto the reference molecule.

        ############################################################

        ############################################################

        pyO3A = rdMolAlign.GetO3A(MOL_Targ, MOL_Ref)

        score = pyO3A.Align()

        F_SCORE.write("{0}\n".format(np.around(score,decimals=2)))

        print("RMD Score: ",score)

        ## Printing the molecule to a pdb file.

        Chem.MolToPDBFile(MOL_Targ,"newMol-round{0:0>4}.pdb".format(i))

    ## Concatenating the pdb files into a single file and removing the 
individual prints.

    COMMAND = "cat newMol*.pdb > all.pdb"

    os.system(COMMAND)

    COMMAND = "rm newMol*pdb"

    os.system(COMMAND)

    ## Generate the histogram of the RMDS scores

    FP_SCORE_H = "RMSD-hist.out"

    COMMAND = "sort -n {0} | uniq -c > {1}".format(FP_SCORE, FP_SCORE_H)

    os.system(COMMAND)




_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss