I would like to identify the atom that used to participate in a bond before the molecule was modified by MoleculeModify[mol, {"DeleteBond", {idx1, idx2}}].

According to the official documentation of MoleculeModify:

When removing or replacing an atom, the number of explicit hydrogen atoms may be adjusted to maintain proper valence. Disconnected hydrogen atoms will be removed.

This suggests to me that I should be able to find the formerly participating atom in the connected component by finding the current bond that features the highest atomic index (mol // AtomList // Length).

My questions are:

  1. Is the new hydrogen atom, added to balance out the valences, always guaranteed to have the highest index in the list of atoms of the connected component? (In other words, is such an approach robust?)

  2. Is there a better way to do this?


As per Jason's request, here is some code to better illustrate the question.

mol = Molecule[

getIdxOfAtomInDeletedBond[mol_] := mol // (
        Bond[{OrderlessPatternSequence[Length @ AtomList @ #, _]}, _] &
    } /* Through /* Apply[Cases] /*
    {First, #[[1, 1]] & /* Min} /* Through

mol //

    (* Original molecule *)
    MoleculePlot[#, {Bond[{30, 31}]}] &,

    (* Randomly selected component after the deletion of the bond *)
    MoleculeModify[#, {"DeleteBond", {30, 31}}] & /*
    WLBugFix`connectedMolComps /* RandomChoice /*
    (MoleculePlot[#, getIdxOfAtomInDeletedBond[#]] &)

} // Through

Please note that WLBugFix``connectedMolComps is the temporary bug fix for ConnectedMoleculeComponents as per this question.

Here is a sample output (sample because there is a random selection from the connected components):

Sample output

In the code above, the indices {30, 31} are provided, but in the actual run of the application, the indices are determined by a number of different probabilistic factors.

In the general case, we cannot expect the same range of atomic indices to be available in the selected component after the deletion of the bond. This is illustrated here by the fact that the molecule on the right cannot possibly have atoms with indices 30 or 31.

However, I would still like to know which atom on the right used to participate in the deleted bond in the parent molecule. It is shown in blue here, and computed by the means detailed in the code.

I just wanted to know if this is a robust way of calculating this, or if there is indeed a better way.

  • 1
    $\begingroup$ Can you give a concrete example that demonstrates the issue? $\endgroup$ – Jason B. Dec 27 '20 at 12:40
  • $\begingroup$ @JasonB. I have updated the question. Thank you for looking at it. $\endgroup$ – Shredderroy Dec 27 '20 at 17:59
  • 1
    $\begingroup$ Thank you - an example is always helpful. One point you are free to ignore - I find your function getIdxOfAtomInDeletedBond a bit difficult to read. The combination of right-composition infix operators with postfix and Through is very unfamiliar to me. I personally find something using modules like this to be more digestible when reading code I didn't write. Lastly, you use Length @ AtomList @ mol where AtomCount is more direct. $\endgroup$ – Jason B. Dec 27 '20 at 20:18

Using the index of the last hydrogen atom to find the index of the atom in the broken bond is not a robust method of keeping track of the atoms before and after a modification.

Hydrogen atoms are often instantiated as needed, especially in molecules created from SMILES strings with implicit hydrogens.

A quick and easy way to keep track of the atoms is to label them before breaking the bond. Ideally we could set some MetaInformation tag on the atoms, and query for it after breaking the bond and getting the connected components. Per-atom meta information isn't supported yet (request it here), so we can just use the mass number to label the atoms.

ClearAll @ mol;
mol = Molecule @ "O=C1C(Cc2ccc(O)cc2)N2C(=O)CCN(C(=O)Nc3ccccc3)C2CN1Cc1cccc2ccccc12";
bondToBreak = {30, 31};
taggedMol = MoleculeModify[mol, {"SetMassNumber", Thread[bondToBreak -> 12]}];
newMols = connectedMoleculeComponents @ MoleculeModify[taggedMol, {"DeleteBond", bondToBreak}];
(MoleculePlot[#1, Atom["MassNumber" -> Except[None]]]&) /@ newMols

enter image description here

Now you can use the pattern Atom["MassNumber" -> Except[None]] to find the labeled atoms after the transformation.

  • $\begingroup$ Thank you very much. I am going to fill out the form to request per-atom meta-information. $\endgroup$ – Shredderroy Dec 27 '20 at 20:57
  • $\begingroup$ I also used the opportunity to request a feature that is at the root of all the gymnastics in this question. $\endgroup$ – Shredderroy Dec 27 '20 at 21:38
  • $\begingroup$ @Shredderroy that’s great. I wondered if this was an XY problem. $\endgroup$ – Jason B. Dec 28 '20 at 0:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.