Dear Andrew,
On Tue, Jan 24, 2012 at 8:08 AM, Andrew Dalke <[email protected]> wrote:
> I wrote some code to split a molecule up by its ring systems to make
> fragments. The code is after the signature of this email.
>
> Structure editing was more difficult than I expected. I had a bond that I
> wanted to break, and to create a dummy "*" atom for each new fragment, to
> indicate the type of connection which was broken; "[1*]" means the other atom
> was aliphatic and "[2*]" means the other atom was aromatic.
>
> For example, in "c1ccccc1O" I want to see "c1ccccc1[1*].[2*]O"
>
> The difficulty was that I needed to keep track of the expected atom index of
> newly added atoms. That is:
>
> # Keep track of the size, because that will be the newly added atom's index
> n = mol.GetNumAtoms()
> emol = Chem.EditableMol(mol)
>
> atom1 = mol.GetAtomWithIdx(atom1_idx)
> atom2 = mol.GetAtomWithIdx(atom2_idx)
> bond_type = mol.GetBondBetweenAtoms(atom1_idx, atom2_idx).GetBondType()
>
> emol.RemoveBond(atom1_idx, atom2_idx)
> add_dummy_atom(emol, n, atom1, atom2, bond_type)
> n += 1
> add_dummy_atom(emol, n, atom2, atom1, bond_type)
> n += 1
>
> where
>
> def add_dummy_atom(emol, new_idx, atom1, atom2, bond_type):
> new_atom = Chem.Atom(0)
> if atom2.GetIsAromatic():
> new_atom.SetMass(2.0)
> else:
> new_atom.SetMass(1.0)
> emol.AddAtom(new_atom)
> emol.AddBond(atom1.GetIdx(), new_idx, bond_type)
>
>
> Do I need to keep track of the expected atom index myself, or is there a
> nicer way to do this?
At the moment you need to keep track yourself. I agree that this is
not particularly nice. It would be much better if AddAtom() (and
AddBond()) were to return the index of the atom/bond just added.
Fortunately this is an easy one to remedy: I will check in a change
this morning so that you can do this:
In [2]: m1 = Chem.MolFromSmiles('CC')
In [3]: em = Chem.EditableMol(m1)
In [4]: em.AddAtom(Chem.Atom(8))
Out[4]: 2
In [5]: em.AddBond(0,2)
Out[5]: 2
> My thought was that I could get the atom index after it has been added to the
> EditableMolecule, as in:
>
> def add_dummy_atom(emol, atom1, atom2, bond_type):
> new_atom = Chem.Atom(0)
> if atom2.GetIsAromatic():
> new_atom.SetMass(2.0)
> else:
> new_atom.SetMass(1.0)
> emol.AddAtom(new_atom)
> new_idx = new_atom.GetIdx() # <----- This always returns 0
> emol.AddBond(atom1.GetIdx(), new_idx, bond_type)
>
> That does not work. The newly added atom's GetIdx() always returns 0.
That's expected behavior. EditableMol.AddAtom copies the atom you give it.
-greg
------------------------------------------------------------------------------
Keep Your Developer Skills Current with LearnDevNow!
The most comprehensive online learning library for Microsoft developers
is just $99.99! Visual Studio, SharePoint, SQL - plus HTML5, CSS3, MVC3,
Metro Style Apps, more. Free future releases when you subscribe now!
http://p.sf.net/sfu/learndevnow-d2d
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
