Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Patch/ambig stereo #826

Merged
merged 2 commits into from
Feb 8, 2022
Merged

Patch/ambig stereo #826

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3365c79
Logging and additional checks where perceiving 2D stereochemistry fro…
johnmay Feb 7, 2022
305079f
Additional strict mode check for stereochemistry, checking the wedge/…
johnmay Feb 7, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
98 changes: 98 additions & 0 deletions base/standard/src/main/java/org/openscience/cdk/geometry/GeometryUtil.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Comparator;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
Expand Down Expand Up @@ -1853,4 +1854,101 @@ public static double[] shiftReactionVertical(IReaction reaction, double[] bounds
}
}

private static class PolarAtomComparator implements Comparator<IAtom> {
private final Point2d center;

public PolarAtomComparator(Point2d center) {
this.center = center;
if (center == null)
throw new NullPointerException("Center point cannot be null!");
}

@Override
public int compare(IAtom aAtom, IAtom bAtom) {
Point2d a = aAtom.getPoint2d();
Point2d b = bAtom.getPoint2d();
if (a == null || b == null)
throw new IllegalArgumentException("Missing 2D coordinated!");

double deltaXa = a.x - center.x;
double deltaXb = b.x - center.x;

if (deltaXa >= 0 && deltaXb < 0) return -1;
if (deltaXa < 0 && deltaXb >= 0) return +1;

double deltaYa = a.y - center.y;
double deltaYb = b.y - center.y;

if (deltaXa == 0 && deltaXb == 0) {
if (deltaYa >= 0 || deltaYb >= 0)
return a.y > b.y ? -1 : +1;
return b.y > a.y ? -1 : +1;
}

// compute the cross product of vectors (center -> a) x (center -> b)
double det = (deltaXa) * (deltaYb) - (deltaXb) * (deltaYa);
if (det < 0) return -1;
if (det > 0) return +1;

// points a and b are on the same line from the center
// check which point is closer to the center
double d1 = (deltaXa) * (deltaXa) + (deltaYa) * (deltaYa);
double d2 = (deltaXb) * (deltaXb) + (deltaYb) * (deltaYb);
return Double.compare(d1, d2);
}
}

private static final class PolarBondComparator implements Comparator<IBond> {

private final IAtom centerAtom;
private final PolarAtomComparator polarAtomCompare;

private PolarBondComparator(IAtom atom) {
this.centerAtom = atom;
if (centerAtom == null)
throw new NullPointerException("Central atom cannot be null!");
this.polarAtomCompare = new PolarAtomComparator(centerAtom.getPoint2d());
}

/**
* Returns the polar ordering of bonds a and b relative to some central
* point.
* @param a one bond
* @param b another bond
* @return -1 a is less than b, 0 a == b, +1 a is greater than b
* @see <a href="http://stackoverflow.com/a/6989383">Sort points in clockwise order, ciamej</a>
*/
@Override
public int compare(IBond a, IBond b) {
IAtom aAtom = a.getOther(centerAtom);
IAtom bAtom = b.getOther(centerAtom);
if (aAtom == null || bAtom == null)
throw new IllegalArgumentException("Bond is not connected to central atom!");
return polarAtomCompare.compare(aAtom, bAtom);
}
}

/**
* Sort the bonds using polar/radial coords around a central atom in 2D.
* @param central the central atom
*/
public static Comparator<IBond> polarBondComparator(final IAtom central) {
return new PolarBondComparator(central);
}

/**
* Sort the atoms using polar/radial coords around a central point in 2D.
* @param central the central atom
*/
public static Comparator<IAtom> polarAtomComparator(final Point2d central) {
return new PolarAtomComparator(central);
}

/**
* Sort the atoms using polar/radial coords around a central atom in 2D.
* @param central the central atom
*/
public static Comparator<IAtom> polarAtomComparator(final IAtom central) {
return polarAtomComparator(central.getPoint2d());
}
}
116 changes: 99 additions & 17 deletions base/standard/src/main/java/org/openscience/cdk/stereo/StereoElementFactory.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@

package org.openscience.cdk.stereo;

import org.openscience.cdk.geometry.GeometryUtil;
import org.openscience.cdk.graph.Cycles;
import org.openscience.cdk.graph.GraphUtil;
import org.openscience.cdk.interfaces.IAtom;
Expand All @@ -32,6 +33,8 @@
import org.openscience.cdk.interfaces.IDoubleBondStereochemistry;
import org.openscience.cdk.interfaces.IStereoElement;
import org.openscience.cdk.interfaces.ITetrahedralChirality;
import org.openscience.cdk.tools.ILoggingTool;
import org.openscience.cdk.tools.LoggingToolFactory;

import javax.vecmath.Point2d;
import javax.vecmath.Point3d;
Expand All @@ -46,6 +49,8 @@
import static org.openscience.cdk.graph.GraphUtil.EdgeToBondMap;
import static org.openscience.cdk.interfaces.IBond.Stereo.DOWN;
import static org.openscience.cdk.interfaces.IBond.Stereo.DOWN_INVERTED;
import static org.openscience.cdk.interfaces.IBond.Stereo.UP;
import static org.openscience.cdk.interfaces.IBond.Stereo.UP_INVERTED;
import static org.openscience.cdk.interfaces.IDoubleBondStereochemistry.Conformation;
import static org.openscience.cdk.interfaces.ITetrahedralChirality.Stereo;

Expand Down Expand Up @@ -94,10 +99,15 @@ public abstract class StereoElementFactory {

protected final Set<Projection> projections = EnumSet.noneOf(Projection.class);

protected boolean strict;

protected final ILoggingTool logger
= LoggingToolFactory.createLoggingTool(StereoElementFactory.class);

/**
* Verify if created stereochemistry are actually stereo-centres.
*/
private boolean check = false;
protected boolean check = false;

/**
* Internal constructor.
Expand Down Expand Up @@ -473,6 +483,21 @@ public StereoElementFactory checkSymmetry(boolean check) {
return this;
}

/**
* Enables stricter stereochemistry checking, specifically tetrahedral
* centres may not be created from inverse down wedges (i.e. Daylight
* style depictions). This also sets that all stereocentres are tested
* for asymmetry.
*
* @return stereo element factory.
*/
public StereoElementFactory withStrictMode()
{
this.check = true;
this.strict = true;
return this;
}

/**
* Create a stereo element factory for creating stereo elements using 2D
* coordinates and depiction labels (up/down, wedge/hatch).
Expand Down Expand Up @@ -544,60 +569,82 @@ ITetrahedralChirality createTetrahedral(int v, Stereocenters stereocenters) {

if (hasUnspecifiedParity(focus)) return null;

if (check) {
if (!stereocenters.isStereocenter(v))
return null;
}

IAtom[] neighbors = new IAtom[4];
int[] elevation = new int[4];

neighbors[3] = focus;
neighbors[3] = focus; // implicit neighbour if needed

boolean nonplanar = false;
int n = 0;

List<IBond> bonds = new ArrayList<>();
for (int w : graph[v]) {
IBond bond = bondMap.get(v, w);
bonds.add(bond);

// wavy bond
if (isUnspecified(bond)) return null;

neighbors[n] = container.getAtom(w);
elevation[n] = elevationOf(focus, bond);

if (elevation[n] != 0) nonplanar = true;
if (elevation[n] != 0)
nonplanar = true;

n++;
}

// too few/many neighbors
if (n < 3 || n > 4) return null;
// too few neighbors
if (n < 3) return null;

// TODO: verify valid wedge/hatch configurations using similar procedure
// to NonPlanarBonds in the cdk-sdg package.
if (strict) {
if (!verifyWedgePattern(focus, n, bonds))
return null;
}

// no up/down bonds present - check for inverted down/hatch
if (!nonplanar) {
if (!nonplanar && !strict) {
int[] ws = graph[v];
for (int i = 0; i < ws.length; i++) {
int w = ws[i];
IBond bond = bondMap.get(v, w);

// we have already previously checked whether 'v' is at the
// 'point' and so these must be inverse (fat-end @
// stereocenter) ala Daylight
// 'point' and so these must be inverse (fat-end of hatched
// wedge is a stereocenter) ala Daylight
if (bond.getStereo() == DOWN || bond.getStereo() == DOWN_INVERTED) {

// we stick to the 'point' end convention but can
// interpret if the bond isn't connected to another
// stereocenter - otherwise it's ambiguous!
if (stereocenters.isStereocenter(w)) continue;
if (stereocenters.stereocenterType(w) != Stereocenters.Stereocenter.Non) {
stereocenters.checkSymmetry();
if (stereocenters.isStereocenter(w)) {
logger.error("Ambiguous down wedge bond between atom indexes ", w, " and ", v);
return null;
}
}

logger.warn("Inverse wedge bond used for stereo at atom index ", v);
elevation[i] = -1;
nonplanar = true;
}
// stereo at the "fat-end" of the bold wedge?
else if (bond.getStereo() == UP || bond.getStereo() == UP_INVERTED) {
logger.warn("Ignoring inverted up wedge bond connected to atom idx=", w);
return null;
}
}

// still no bonds to use
if (!nonplanar) return null;
}

// still no bonds to use
if (!nonplanar) return null;

int parity = parity(focus, neighbors, elevation);

if (parity == 0) return null;
Expand All @@ -607,6 +654,44 @@ ITetrahedralChirality createTetrahedral(int v, Stereocenters stereocenters) {
return new TetrahedralChirality(focus, neighbors, winding);
}

private boolean isOkay(int a, int b) {
return a == 0 || b == 0 || a == b;
}

// check some obvious stereo chemistry errors, see the InChI
// technical manual "Definition of 2D drawing correctness"
private boolean verifyWedgePattern(IAtom focus, int n, List<IBond> bonds) {
if (n == 3) {
// no sort needed
int ref = 0;
for (IBond bond : bonds) {
int curr = elevationOf(focus, bond);
if (!isOkay(ref, curr)) {
logger.error("Badly drawn stereochemistry with 3 neighbours, up/down bonds should not be mixed!");
return false;
} else
ref = curr;
}
} else { // n == 4
bonds.sort(GeometryUtil.polarBondComparator(focus));
int ref = 0;
for (IBond bond : bonds) {
int curr = elevationOf(focus, bond);
if (curr != 0) {
if (ref != 0 && ref != curr) {
logger.error("Badly drawn stereochemistry with 4 neighbours, up/down bonds should alternate!");
return false;
}
else {
ref = curr;
}
}
ref = -ref; // flip for next check
}
}
return true;
}

private static boolean isWedged(IBond bond) {
switch (bond.getStereo()) {
case UP:
Expand Down Expand Up @@ -1082,9 +1167,6 @@ ITetrahedralChirality createTetrahedral(int v, Stereocenters stereocenters) {
// too few/many neighbors
if (n < 3 || n > 4) return null;

// TODO: verify valid wedge/hatch configurations using similar procedure
// to NonPlanarBonds in the cdk-sdg package

int parity = parity(neighbors);

Stereo winding = parity > 0 ? Stereo.ANTI_CLOCKWISE : Stereo.CLOCKWISE;
Expand Down