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 <[email protected]>
Sent: Thursday, August 22, 2019 11:47:19 AM
To: Illimar Hugo Rekand; [email protected]
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
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
