Rules for Interpreting Up/Down Wedge Bonds

Yesterday I was reminded of an old anecdote about a maths professor. A professor was lecturing an auditorium and writing down a proof. They proclaim at one part "and obviously x infers y". A student raises their hand and asks, "is it obvious?". The professor then studies the equation for 30 minutes until stating "Yes, it is obvious" and continued on with the proof.

chemists and SMILES experts, is the bridge carbon point towards us, or away from us? FC(F)(F)c1ccc(nc1)O[C@H]1C[C@H](C[C@H]12)C([2H])([2H])N2C(=O)c1cccc(F)c1-c1=Nccc=N1 https://t.co/6WwEvUSbkl (cc @jwmay @cdsouthan) pic.twitter.com/93q8TpZA58

— Egon Willighⓐgen (@egonwillighagen) September 1, 2019

Egon asked is the bridge in the following compound up or down? I replied obviously it's pointing down. Since it might not be obvious to everyone I thought I'd explain the three rules to easily assign up/down wedges to other bonds around a stereocentre in 2D.

The insight actually comes from the algorithm used to assign up/down wedges to a depiction. Since only a handful of people have ever had to write such an algorithm I'm not sure how common knowledge it is (i.e. is it actually obvious?).

Rule 1 (D₄)

A tetrahedral centre with four bonds must have alternating up/down bonds. Therefore no mater what the angle if one bond is labelled as up, we know the bond opposite it must also be up, and the two either side must be down (inverse of up).

Rule 2 (D₃)

For three bonds, when bonds are spaced evenly (i.e. all angles < 180 degrees) then all bonds are the same direction, all up, or all down.

Rule 3 (D₃)

For three bonds, when an angle > 180 exists then think about it like the D₄ case with one neighbour missing. The "outside" bonds are opposite direction to the "middle" one.

Exceptions

Like all rules there are exceptions, a well know ambiguity is when we have three bonds and the angle is exactly 180 degrees:

The problem here is you could move the central atom slightly to apply either Rule 2 or Rule 3. Some chemistry toolkits will refuse to read this others will side on the more likely interpretation (Rule 3).

A bigger and perhaps more common issue is mixing up/down wedges with perspective projection. More precisely MDL (and then SYMYX, Accelrys, now BIOVIA) had something known as the "triangle rule". The idea was if you were looking at a molecule in 3D the lengths of the bond would indicate which way round you were looking at it. They imposed this concept on 2D interpretation.

In practice what the this means is these two structures are read as different enantiomers (by BIOVIA) depending on whether the H is inside or outside the "triangle":

You're unlikely to encounter such cases except when a projection is involved. For example for the bridged system pictured below, perspective has been used and we may end up with the H within the "triangle". Note it's the point stored in the file for the atom not the actual "H" glyph that maters.

This isn't to say projections are bad, only that mixing perspective with up/down wedges can be problematic.

Creating Chemical Structure Animations

I've just got back from the Eighth Joint Sheffield Conference on Chemoinfomatics where I presented about the technical details of subgraph isomorphism algorithms. It was a great conference (as usual) with good science, interesting posters and lots of fun. Noel was live tweeting the whole thing so check out the #Shef2019 hashtag if you want to catch up.

To help explain the algorithms in my talk I created some animations (as videos) that showed the backtracking search procedures. After the talk several delegates asked how I created these so I thought I do a quick blog post on it (and also to remind me in future how to do it again). I'd done something similar before to demonstrate the Sayle and Delany algorithm for tautomer enumeration. Here's the PDF (and PPTX with the videos) of my Sheffield talk.

CDK + ffmpeg

The general idea is to generate a bunch of PNGs and then stick them together with the ffmpeg command line tool. In older blog posts I used to generate a GIF but it turns out mp4 compresses better with higher quality.

Step 1: Generating the PNGs with CDK

The example code below loads the SMILES for indole and then loops around highlighting one atom. Other than some normal params we also set the symbol visibility. By default the depiction will add in symbols for highlighted carbons, we can turn this off by overriding the parameter.

String dest = "/tmp/example1";
IChemObjectBuilder bldr   = SilentChemObjectBuilder.getInstance();
SmilesParser       smipar = new SmilesParser(bldr);
IAtomContainer     mol    = smipar.parseSmiles("[nH]1ccc2c1cccc2");
DepictionGenerator dg = new DepictionGenerator().withZoom(5)
                                                .withSize(1280, 720)
                                                .withParam(StandardGenerator.Visibility.class,
                                                           SymbolVisibility.iupacRecommendationsWithoutTerminalCarbon())
                                                .withOuterGlowHighlight();
new File(dest).mkdirs();
for (int i=0; i<100; i++) {
    IAtom atm = mol.getAtom(i % mol.getAtomCount());
    dg.withHighlight(Collections.singleton(atm), Color.RED)
      .depict(mol)
      .writeTo(String.format("%s/frame.%03d.png", dest, i));
}

Step 2: Stick them all together

$ ffmpeg -r 12/1 -start_number 0 -i /tmp/example1/frame.%03d.png \
  -c:v libx264 -r 30 -pix_fmt yuv420p example1.mp4

The key parameters here are the output name (last arg) and input name template (-i /tmp/example1/frame.%03d.png). I zero padded the numbers so they will be ordered correctly when alphabetically sorted and match this in the argument. You don't need to zero-pad but it means they should then alphabetical in your OS file system which is handy. The first -r is used to set the input frame rate to 12/1 (12 per 1 second).

Okay so how did it turn out....

Download: example1.mp4

Not bad but the inter-frame alignment is off due to the different size caused by the moving outer-glow. In future I might have an anchor attribute that would lock the depiction in place relative to some atom but that's quite specialised. We can fix the shifting by highlighting all other atoms to the same as the background. This is possible with the high-level APIs but it's easy enough just to set the field on the atom:

for (int i=0; i<100; i++) {
    IAtom target = mol.getAtom(i % mol.getAtomCount());
    // reset all
    for (IAtom a : mol.atoms())
        a.setProperty(StandardGenerator.HIGHLIGHT_COLOR,
                      Color.WHITE);
    target.setProperty(StandardGenerator.HIGHLIGHT_COLOR,
                       Color.RED);
    dg.depict(mol)
      .writeTo(String.format("%s/frame.%03d.png", dest, i));
}

Download: example2.mp4

Add/Remove Atoms and Bonds?

Because of the alignment issue we need to use some tricks to add/remove atoms and bonds. Essentially you draw the whole thing and then hide the parts you don't want by setting them to the background color. For example:

int frameId = 0;
for (IAtom atom : mol.atoms())
    hide(atom);
for (IBond bond : mol.bonds())
    hide(bond);
// add bonds
for (int i = 0; i < mol.getBondCount(); i++) {
    IBond bnd = mol.getBond(i);
    show(bnd);
    show(bnd.getBegin());
    show(bnd.getEnd());
    dg.depict(mol)
      .writeTo(String.format("%s/frame.%03d.png", dest, ++frameId));
}
// hold for 12 frames
for (int i = 0; i < 12; i++)
    dg.depict(mol)
      .writeTo(String.format("%s/frame.%03d.png", dest, ++frameId));
// remove bonds
for (int i = 0; i < mol.getBondCount(); i++) {
    IBond bnd = mol.getBond(mol.getBondCount()-i-1);
    hide(bnd);
    if (!visible(bnd.getBegin()))
        hide(bnd.getBegin());
    if (!visible(bnd.getEnd()))
        hide(bnd.getEnd());
    dg.depict(mol)
      .writeTo(String.format("%s/frame.%03d.png", dest, ++frameId));
}

Where show/hide are:

private static void hide(IChemObject x) {
    x.setProperty(StandardGenerator.HIGHLIGHT_COLOR,
                  Color.WHITE);
}

private static void show(IChemObject x) {
    x.setProperty(StandardGenerator.HIGHLIGHT_COLOR,
                  Color.BLACK);
}

Producing:

Download: example3.mp4

That's just about it, the CDK depiction is quite configurable but not really intended to be a general purpose drawing/animation tool. However using some tricks you can work around some quirks and get some nice results. If you make animations let me know as I'd love to see them!

CDK Depict on HTTPS

Just a quick post to say CDK Depict is now using HTTPS https://simolecule.com/cdkdepict/depict.html. The main reason for this was Blogger stopped allowing image links to HTTP resources. In general browsers are being more fussy about non HTTPS content.

I used LetsEncrypt that turned out to be very easy to configure with TomCat.

Step 1

Install the certbot utility and use it generate a certificate.

$ sudo certbot certonly

Step 2

Configure TomCat 8+ connectors. This used to be more complex on older TomCat servers with the need to generate a separate keystore. Editing $CATALINA_HOME/confg/server.xml we configure the base connected, redirectPort is changed from 8443 to 443 (and 8080 to 80).

<Connector port="80" protocol="HTTP/1.1"
           connectionTimeout="20000"
           redirectPort="443" />

We also configure SSL connector, using port 443, change to NIO based protocol (the default requires extra native library) org.apache.coyote.http2.Http2Protocol, and set the file paths to the .pem files generated by certbot.

<Connector port="443" 
           protocol="org.apache.coyote.http11.Http11NioProtocol"
           maxThreads="150" SSLEnabled="true" >
  <UpgradeProtocol className="org.apache.coyote.http2.Http2Protocol" />
  <SSLHostConfig>
    <Certificate certificateKeyFile="/etc/letsencrypt/live/www.simolecule.com/privkey.pem"
                 certificateFile="/etc/letsencrypt/live/www.simolecule.com/cert.pem"
                 certificateChainFile="/etc/letsencrypt/live/www.simolecule.com/chain.pem"
                 type="RSA" />
  </SSLHostConfig>
</Connector>       

Step 3 (optional)

If a client tries to visit the HTTP site we want to redirect them to HTTPS. To do this we edit $CATALINA_HOME/confg/web.xml adding this section to the end of the <web-app> block

<security-constraint>
  <web-resource-collection>
    <web-resource-name>Entire Application</web-resource-name>
    <url-pattern>/*</url-pattern>
  </web-resource-collection>
  <user-data-constraint>
    <transport-guarantee>CONFIDENTIAL</transport-guarantee>
  </user-data-constraint>
</security-constraint>

Bit packing for fast atom type assignment

Many cheminformatics algorithms perform some form of atom typing as their first step. The atom types you need and the granularity depend on the algorithm. At the core of all atom typing lies some form of decision tree, usually manifesting as a complex if-else cascade.

Using an elegant technique Roger showed me a few years ago, you can replace this if-else cascade with a table/switch lookup. This turns out to be very clean, easy to extend, and efficient. The technique relies on first computing a value that captures the bond environment around an atom and packing this in to a single value.

Here's what it looks like:

Algorithm 1

int btypes = atom.getImplicitHydrogenCount();
for (IBond bond : atom.bonds()) {
  switch (bond.getOrder()) {
    case SINGLE: btypes += 0x0001; break;
    case DOUBLE: btypes += 0x0010; break;
    case TRIPLE: btypes += 0x0100; break;
    default:     btypes += 0x1000; break;
  }
}

After the value btypes has been calculated (once for each atom) it contains the count of single, double, triple, and other bonds. We can inspect each of these counts individually my masking of the relevant bits or the entire value, for example:

Algorithm 2

switch (btypes) {
 case 0x004: // 4 single bonds, e.g. Sp3 carbon
  break;
 case 0x012: // 1 double bond, 2 single bonds e.g. Sp2 carbon
  break;
 case 0x020: // 2 double bonds e.g. Sp cumulated carbon
  break;
 case 0x101: // 1 triple bond, 1 single bond e.g. Sp triple bonded carbon
  break;
}

Using a nibble allows us to store numbers up to 16 (2⁴) - more than enough for any sane chemistry. In the example above I shoved default bonds under an 'other' category but of course you could extend it to handle quadruple bonds and even additional properties of the bonds:

Algorithm 3

int btypes = atom.getImplicitHydrogenCount();
for (IBond bond : atom.bonds()) {
  switch (bond.getOrder()) {
    case SINGLE: 
      if (bond.isAromatic())
        btypes += 0x010001; 
      else
        btypes += 0x000001;
      break;
    case DOUBLE: 
      if (bond.isAromatic())
        btypes += 0x010010; 
      else if (isOxygen(bond.getOther(atom)))
        btypes += 0x100010; // dbs to oxygens
      else
        btypes += 0x000010;
      break;
    case TRIPLE: 
       btypes += 0x000100;
       break;
  }
}

RDKit Reaction SMARTS

Update 19/09/22

I believe I was wrong and Daylight did allow SMIRKS to be run backwards and that the "direction" parameter controls it ("personal communication" from Roger).. The key insight is the documentation notes SMARTS patterns should only appear on atoms that don't have bond changes around them. From the man page:

"Atomic expressions may be any valid atomic SMARTS expression for nodes where the bonding (connectivity & bond order) doesn't change.1 Otherwise, the atomic expressions must be valid SMILES."

Whether this make sense or not for writing portable SMIRKS is questionable.

This is supplementary to my original grumble on ambiguous naming. I still believe RDKit should just call it SMIRKS or at least something which didn't already mean something else (e.g. "RdSmirks"). I do agree there is a lot of overlap between SMARTS and SMIRKS but conflating these terms is problematic. Shortly after I originally wrote this original post, SMIRKS Native Open Force Field (SMIRNOFF) was published. In that work SMARTS is used, but they called it SMIRKS. There is some wiggling attempted by the authors that they use "SMIRKS features" of atom maps. However it is incorrect that these were a SMIRKS feature as in what Daylight called Reaction SMARTS, see the table "Examples of Reaction SMARTS" (SMARTS Theory Page). Of course SMIRNOFF was just too good a pun to change, but it again creates more confusion.

I recently ran some comparisons/benchmarks on different SMIRKS implementations, everyone's semantics are consistently inconsistent and the community would benefit by an effort to standardise. Perhaps then I can convince RDKit that they really do have a (good) SMIRKS implementation and it would be less confusing to just call it that.

---