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