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 <paolo.tosco.m...@gmail.com>
Sent: Thursday, August 22, 2019 11:47:19 AM
To: Illimar Hugo Rekand; rdkit-discuss@lists.sourceforge.net
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
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
_______________________________________________
Rdkit-discuss mailing list
Rdkit-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
