Re: [Rdkit-discuss] Setting custom coordinates for new atoms

[Paolo T. <pao...@gm...>]

Hi Illimar,

The problem here is that you are adding atoms to your RWMol, but the 
Conformer object that you are trying to access with the newly added atom 
indices is still tied to the original molecule which does not have the 
newly added atoms, hence your indices are out of bounds. If you 
reference the Conformer object associated to the RWMol this won't 
happen. I have slightly modifed your code below:

from  rdkit  import  Chem
from  rdkit.Chem  import  AllChem,  rdFreeSASA,  RDConfig,  \
                        rdmolops,  rdMolTransforms,  rdGeometry

def  block_apolar_N_atoms(mol):
     """
Method:
The unit vector is found from two neighbouring atoms.
This unit vector is scaled to the appropiate distance
of the two new dummy atoms from the N-atom (2.4 Å)
The two atoms are placed at the N-atom coord + and -
the scaled unit vector
"""
     #places two C-atoms 2.4 Å above and below aromatic N-atoms.
     mw  =  Chem.RWMol(mol)
     conf  =  mw.GetConformer()
     for  atom  in  mol.GetAtoms():
         idx  =  0
         if  atom.GetSymbol()  ==  "N"  and  atom.GetIsAromatic():
             idx  =  atom.GetIdx()
             print(idx)
             print("Found aromatic N")
             atom_3d_point  =  conf.GetAtomPosition(idx)
             n  =  0
             vector_list  =  []
             #need two neighbors to find cross product -> unit vector
             while  n  <  2:
                 for  nbr  in  atom.GetNeighbors():
                     nbridx  =  nbr.GetIdx()
                     nbr_3d_point  =  conf.GetAtomPosition(nbridx)
                     vector  =  nbr_3d_point  -  atom_3d_point
                     print("vector",  vector,  list(vector))
                     vector_list.append(vector)
                     n  +=  1
             #the unit vector is perpendicular to the two vectors
             unit_vector  =  rdGeometry.Point3D.CrossProduct(
                 vector_list[0],  vector_list[1])
             #this is usually ~1.72 Å for a N in an aromatic six-membered ring
             length_unit_vector  =  rdGeometry.Point3D.Length(unit_vector)
             print("length:",  length_unit_vector,  "Å")
             #let's figure out how much the unit vector needs to be scaled
             scalar  =  2.4  /  length_unit_vector
             #let's scale the vector
             scaled_vector  =[i  *  scalar  for  i  in  list(unit_vector)]
             #make a Point3D object. Easier to handle in RDKit
             scaled_vector_3d  =  rdGeometry.Point3D(
                 scaled_vector[0],  scaled_vector[1],  scaled_vector[2])
             A1  =  atom_3d_point  +  scaled_vector_3d  # New atom coordinate
             A2  =  atom_3d_point  -  scaled_vector_3d  # New atom coordinate
             print("A1",  list(A1),  A1,"A2",  list(A2),  A2)
             A1_idx  =  mw.AddAtom(Chem.Atom(6))
             print("A1_ix",  A1_idx)
             print(type(A1_idx))
             c  =  conf.SetAtomPosition(A1_idx,  A1)
             A2_idx  =  mw.AddAtom(Chem.Atom(6))
             print("A2_idx",  A2_idx)

prot  =  Chem.MolFromPDBFile("C:/data/paolo/support/illimar/3etr.pdb")

block_apolar_N_atoms(prot)

Cheers,
p.

On 26/08/2019 02:13, Illimar Hugo Rekand wrote:
> Hello, everyone and thanks for the help.
>
>
> I tried out implementing the functionality in my code, but after running the following function:
>
>
> #!/usr/bin/env python3.7
>
> from rdkit import Chem
> from rdkit.Chem import AllChem, rdFreeSASA, RDConfig, rdmolops, rdMolTransforms, rdGeometry
>
> def block_apolar_N_atoms(mol, conf): #places two C-atoms 2.4 Å above and below aromatic N-atoms.
>      """
>      Method:
>      The unit vector is found from two neighbouring atoms. This unit vector is scaled to the appropiate distance of the two new dummy atoms from the N-atom (2.4 Å)
>      The two atoms are placed at the N-atom coord + and - the scaled unit vector
>      """
>      mw = Chem.RWMol(mol)
>      for atom in mol.GetAtoms():
>          idx = 0
>          if atom.GetSymbol() == "N" and atom.GetIsAromatic():
>              idx = atom.GetIdx()
>              print(idx)
>              print("Found aromatic N")
>              atom_3d_point = conf.GetAtomPosition(idx)
>              n = 0
>              vector_list = []
>              while n < 2: #need two neighbors to find cross product -> unit vector
>                  for nbr in atom.GetNeighbors():
>                      nbridx = nbr.GetIdx()
>                      nbr_3d_point = conf.GetAtomPosition(nbridx)
>                      vector = nbr_3d_point - atom_3d_point
>                      print("vector", vector, list(vector))
>                      vector_list.append(vector)
>                      n += 1
>              unit_vector = rdGeometry.Point3D.CrossProduct(vector_list[0], vector_list[1]) #the unit vector is perpendicular to the two vectors
>              length_unit_vector = rdGeometry.Point3D.Length(unit_vector) #this is usually ~1.72 Å for a N in an aromatic six-membered ring
>              print("length:", length_unit_vector, "Å")
>              scalar = 2.4 / length_unit_vector #let's figure out how much the unit vector needs to be scaled
>              scaled_vector =[i * scalar for i in list(unit_vector)] #let's scale the vector
>              scaled_vector_3d = rdGeometry.Point3D(scaled_vector[0], scaled_vector[1], scaled_vector[2]) #make a Point3D object. Easier to handle in RDKit
>              A1 = atom_3d_point + scaled_vector_3d # New atom coordinate
>              A2 = atom_3d_point - scaled_vector_3d # New atom coordinate
>              print("A1", list(A1), A1,"A2", list(A2), A2)
>              A1_idx = mw.AddAtom(Chem.Atom(6))
>              print("A1_ix", A1_idx)
>              print(type(A1_idx))
>              c = conf.SetAtomPosition(A1_idx, A1)
>              A2_idx = mw.AddAtom(Chem.Atom(6))
>              print("A2_idx", A2_idx)
>
>
> prot = Chem.MolFromPDBFile("./3etr.pdb")
> protconf = prot.GetConformer()
>
> block_apolar_N_atoms(prot, protconf)
>
>
>
> I get the following error:
>
>
> RuntimeError: Pre-condition Violation
>
>      Violation occurred on line 51 in file Code/GraphMol/Conformer.cpp
>      Failed Expression: dp_mol->getNumAtoms() == d_positions.size()
>      RDKIT: 2019.09.1dev1
>      BOOST: 1_67
>
> Hoping to hear from you soon!
>
>
> Illimar Rekand
> Ph.D. candidate,
> Brenk-lab, Haug-lab
> Department of Biomedicine
> Department of Chemistry
> University of Bergen
>
>
> ________________________________
> From: Paolo Tosco <pao...@gm...>
> Sent: Thursday, August 22, 2019 11:47:19 AM
> To: Illimar Hugo Rekand; rdk...@li...
> Subject: Re: [Rdkit-discuss] Setting custom coordinates for new atoms
>
> Hi Illimar,
>
> AddAtom() will return the index i of the added atom, then you can call
> SetAtomPosition on that index on the molecule conformer and pass a
> Point3D with the desired coordinates:
>
> conf.SetAtomPosition(i, Point3D(x, y, z))
>
> Cheers,
> p.
>
> On 08/22/19 09:24, Illimar Hugo Rekand wrote:
>> Hello, everyone
>>
>>
>> I'm wondering whether there is a way to set custom coordinates to an atom in a conformer?
>>
>> In particular I'm interested in using the AddAtom() function in the RWMol class to place a new dummy atom in a PDB-file.
>>
>>
>> Illimar Rekand
>> Ph.D. candidate,
>> Brenk-lab, Haug-lab
>> Department of Biomedicine
>> Department of Chemistry
>> University of Bergen
>>
>>
>>
>> _______________________________________________
>> Rdkit-discuss mailing list
>> Rdk...@li...
>> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
>
>
> _______________________________________________
> Rdkit-discuss mailing list
> Rdk...@li...
> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss