Do not allow the minimizer to flip constrained fragments

[ptosco]

Hello,

working with the RDKit interface to CoordGen I noticed that when I do a constrained depiction in RDKit using rdDepictor.GenerateDepictionMatching2DStructure and my constrained scaffold consists of two ring systems connected by a rotatable bond, CoordGen may flip the two rings with respect to each other in the attempt to escape a steric clash. This is a simple reproducible:

from rdkit import Chem
from rdkit.Chem import rdDepictor
from rdkit.Chem.Draw import IPythonConsole, MolsToGridImage
IPythonConsole.drawOptions.addAtomIndices = True

simple_scaffold = Chem.MolFromMolBlock("""
  MJ201100                      

 17 19  0  0  0  0  0  0  0  0999 V2000
   -2.9134    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.7383    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -4.1508   -0.5134    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.7384   -1.2278    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.9134   -1.2278    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.5008   -0.5133    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.6583    0.9857    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.3257    1.4706    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.9932    0.9857    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6758   -0.5133    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8737    1.2405    0.0000 *   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2634   -1.2277    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -0.4383   -1.2277    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0257   -0.5133    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.4382    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2633    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6758    0.9156    0.0000 *   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0  0  0  0
  2  3  1  0  0  0  0
  3  4  2  0  0  0  0
  4  5  1  0  0  0  0
  5  6  2  0  0  0  0
  6  1  1  0  0  0  0
  1  7  1  0  0  0  0
  6 10  1  0  0  0  0
  7  8  2  0  0  0  0
  2  9  1  0  0  0  0
  8  9  1  0  0  0  0
  7 11  1  0  0  0  0
 13 14  2  0  0  0  0
 14 15  1  0  0  0  0
 15 16  2  0  0  0  0
 16 10  1  0  0  0  0
 12 13  1  0  0  0  0
 10 12  2  0  0  0  0
 16 17  1  0  0  0  0
M  END
""")

simple_mol = Chem.MolFromMolBlock("""
  MJ201100                      

 18 20  0  0  0  0  0  0  0  0999 V2000
   -2.9134    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.7383    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -4.1508   -0.5134    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.7384   -1.2278    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.9134   -1.2278    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.5008   -0.5133    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.6583    0.9857    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.3257    1.4706    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.9932    0.9857    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6758   -0.5133    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8737    1.2405    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2634   -1.2277    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
   -0.4383   -1.2277    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0257   -0.5133    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.4382    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.2633    0.2011    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6758    0.9156    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.4209    1.7002    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0  0  0  0
  2  3  1  0  0  0  0
  3  4  2  0  0  0  0
  4  5  1  0  0  0  0
  5  6  2  0  0  0  0
  6  1  1  0  0  0  0
  1  7  1  0  0  0  0
  6 10  1  0  0  0  0
  7  8  2  0  0  0  0
  2  9  1  0  0  0  0
  8  9  1  0  0  0  0
 13 14  2  0  0  0  0
 14 15  1  0  0  0  0
 15 16  2  0  0  0  0
 16 10  1  0  0  0  0
 12 13  1  0  0  0  0
 10 12  2  0  0  0  0
  7 11  1  0  0  0  0
 16 17  1  0  0  0  0
 17 18  1  0  0  0  0
M  END
""")

rdDepictor.SetPreferCoordGen(True)
rdDepictor.GenerateDepictionMatching2DStructure(simple_mol, simple_scaffold);

MolsToGridImage((simple_mol, simple_scaffold), subImgSize=(500,500))

I think this is undesirable as the resulting 2D layout does not honor the original constraint anymore. I did a bit of investigation and I have seen that while sketcherMinimizer::alignWithParentDirectionConstrained indeed prevents the flip at the CoordgenMinimizer::buildMoleculeFromFragments stage, at the CoordgenMinimizer::avoidClashes stage the flip may still take place in an attempt to escape the clash between the amino group and the ethyl group.

With the change proposed in this PR, that prevents the possibility to flip in CoordgenMinimizer::flipFragments also for constrained fragments and not only for fixed fragments, the above flip does not take place anymore, but if the clash is quite bad the molecule is still allowed to bend to alleviate the clash. This still does not look great as it would be better to rather rotate by 180 degrees the ethyl substituent around the 15-16 bond:

You may have better solutions to this issue than the one proposed in this PR, but I hope what I described makes sense to you. When I tried to switch the RDKit CoordGen interface to use fixed rather than constrained results I got really bad layouts, so I guess that's not an option at the moment.

In passing, I have also introduced a change into the sketcherMinimizerFragment::countHeavyAtoms, that does not seem to do what is written on the tin at the moment :-).

[d-b-w]

Thanks for submitting this, @ptosco . Rachel/Nic/I will take a look at this in the next couple of days.

The CI failure isn't real - we disabled appveyor.

[d-b-w]

@ptosco - This is a great bug report and solid proposed fix. I think we're going to iterate on it slightly, but we'll either accept this or propose an alternate fix in the next week or so. (I think that we can slightly improve on the treatment of that ethyl)

[rachelnwalker]

Thank you for submitting this! It looks like the 16-17 bond isn't flipping around the 15-16 bond because the 16-17 bond is considered a constrained fragment. I tried introducing a fully_constrained field to sketcherMinimizerFragment, which is set to true when all atoms in a given fragment are constrained. If we do something like this in CoordgenMinimizer::flipFragments:

for (auto fragment : fragments) {
    if (!fragment->fixed && !fragment->fully_constrained) {
        ...
    }
}

the 16-17 bond will flip around the 15-16 bond, as shown below.

As @ZontaNicola mentioned, it may sometimes be necessary to flip constrained fragments. Maybe we could specifically penalize flipping fully constrained fragments? Or degrees of freedom that only affect constrained atoms?

[ptosco]

@rachelnwalker Thank you Rachel for looking into this. I have just updated my PR taking into account @ZontaNicola's remark on countHeteroatoms and included a possible fix for your finding: a fragment should not be considered constrained if a single atom is contrained, as that's just the attachment point. With that fix, I get a much better layout:

[d-b-w]

Paolo - We're going to merge https://github.com/schrodinger/coordgenlibs/pull/93 instead of this. It's a riff on this with some of Nic and Rachel's ideas. Once again - I really appreciate the work, and feel silly that I didn't just merge it before we worked on #93.