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.

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

Additional stereo wedge checks and corner cases. We need to check the… #839

Merged
merged 1 commit into from
Feb 12, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ public abstract class StereoElementFactory {
/** A bond map for fast access to bond labels between two atom indices. */
protected final EdgeToBondMap bondMap;


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

protected boolean strict;
Expand Down Expand Up @@ -460,14 +460,14 @@ abstract ExtendedCisTrans createExtendedCisTrans(List<IBond> bonds,

/**
* Indicate that stereochemistry drawn as a certain projection should be
* interpreted.
* interpreted.
*
* <pre>{@code
* StereoElementFactory factory =
* StereoElementFactory factory =
* StereoElementFactory.using2DCoordinates(container)
* .interpretProjections(Projection.Fischer, Projection.Haworth);
* }</pre>
*
*
* @param projections types of projection
* @return self
* @see org.openscience.cdk.stereo.Projection
Expand Down Expand Up @@ -661,19 +661,69 @@ private boolean isOkay(int a, int 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) {
bonds.sort(GeometryUtil.polarBondComparator(focus));
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;
double angle = getMaxSweep(focus, bonds);
double delta = angle - Math.PI;
double threshold = 0.01;
// largest angle between 2 neighbours is > 180 => wedges should
// alternate
if (delta > threshold) {
int ref = 0;
for (IBond bond : bonds) {
int curr = elevationOf(focus, bond);
if (!isOkay(ref, curr)) {
logger.error("Invalid wedge pattern, up/down bonds should be mixed when there is an acute angle!");
return false;
} else
ref = curr;
ref = -ref;
}
}
// larges angle between 2 neighbours is < 180 => all wedges
// should be the same
else if (delta < -threshold) {
int ref = 0;
for (IBond bond : bonds) {
int curr = elevationOf(focus, bond);
if (!isOkay(ref, curr)) {
logger.error("Invalid wedge pattern, up/down bonds should be same when there is not an acute angle!");
return false;
} else
ref = curr;
}
} else {
// 180-degrees, check where the wedge is
if (bonds.size() != 3)
throw new IllegalArgumentException("3 bonds only!");
Vector2d v1 = toUnitVector(focus, bonds.get(0).getOther(focus));
Vector2d v2 = toUnitVector(focus, bonds.get(1).getOther(focus));
Vector2d v3 = toUnitVector(focus, bonds.get(2).getOther(focus));
String ambiuousStereoMesg = "Ambiguous stereochemistry - 3 neighbours and two";
if (Math.abs(signedAngle(v1,v2) - Math.PI) < threshold) {
if (elevationOf(focus, bonds.get(0)) == 0 &&
elevationOf(focus, bonds.get(1)) == 0 &&
elevationOf(focus, bonds.get(2)) != 0) {
logger.error(ambiuousStereoMesg);
return false;
}
} else if (Math.abs(signedAngle(v2,v3) - Math.PI) < threshold) {
if (elevationOf(focus, bonds.get(0)) != 0 &&
elevationOf(focus, bonds.get(1)) == 0 &&
elevationOf(focus, bonds.get(2)) == 0) {
logger.error(ambiuousStereoMesg);
return false;
}
} else if (Math.abs(signedAngle(v3,v1) - Math.PI) < threshold) {
if (elevationOf(focus, bonds.get(0)) == 0 &&
elevationOf(focus, bonds.get(1)) != 0 &&
elevationOf(focus, bonds.get(2)) == 0) {
logger.error(ambiuousStereoMesg);
return false;
}
}
}
} else { // n == 4
bonds.sort(GeometryUtil.polarBondComparator(focus));
int ref = 0;
for (IBond bond : bonds) {
int curr = elevationOf(focus, bond);
Expand All @@ -692,6 +742,26 @@ private boolean verifyWedgePattern(IAtom focus, int n, List<IBond> bonds) {
return true;
}

private double signedAngle(Vector2d a, Vector2d b) {
double angle = Math.atan2(a.x*b.y-a.y*b.x, a.x*b.x+a.y*b.y);
return angle >= 0 ? (2*Math.PI)-angle : -angle;
}

private double max(double a, double b, double c) {
return Math.max(a, Math.max(b, c));
}

private double getMaxSweep(IAtom atom, List<IBond> bonds) {
if (bonds.size() != 3)
throw new IllegalArgumentException("3 bonds only!");
Vector2d v1 = toUnitVector(atom, bonds.get(0).getOther(atom));
Vector2d v2 = toUnitVector(atom, bonds.get(1).getOther(atom));
Vector2d v3 = toUnitVector(atom, bonds.get(2).getOther(atom));
return max(signedAngle(v1, v2),
signedAngle(v2, v3),
signedAngle(v3, v1));
}

private static boolean isWedged(IBond bond) {
switch (bond.getStereo()) {
case UP:
Expand Down Expand Up @@ -1049,6 +1119,12 @@ private Point2d toUnitVector(Point2d from, Point2d to) {
return new Point2d(v2d);
}

private Vector2d toUnitVector(IAtom from, IAtom to) {
if (from.equals(to))
return new Vector2d(0, 0);
return new Vector2d(toUnitVector(from.getPoint2d(), to.getPoint2d()));
}

/**
* Compute the signed volume of the tetrahedron from the planar points
* and elevations.
Expand Down Expand Up @@ -1292,7 +1368,7 @@ ExtendedTetrahedral createExtendedTetrahedral(int v, Stereocenters stereocenters
IAtom focus = container.getAtom(v);

if (hasUnspecifiedParity(focus)) return null;

if (container.getConnectedBondsCount(focus) != 2)
return null;

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
import org.openscience.cdk.interfaces.IStereoElement;
import org.openscience.cdk.interfaces.ITetrahedralChirality;
import org.openscience.cdk.io.MDLV2000Reader;
import org.openscience.cdk.io.MDLV2000Writer;
import org.openscience.cdk.silent.Atom;
import org.openscience.cdk.silent.AtomContainer;
import org.openscience.cdk.silent.SilentChemObjectBuilder;
Expand Down Expand Up @@ -626,11 +627,10 @@ public void badWedgePatternWithThreeNeighbors() throws CDKException {
m.addBond(1, 0, IBond.Order.SINGLE, IBond.Stereo.DOWN); // CH-OH
m.addBond(1, 2, IBond.Order.SINGLE, IBond.Stereo.UP); // CH-CH2
m.addBond(1, 3, IBond.Order.SINGLE); // CH-Cl

StereoElementFactory factory = StereoElementFactory.using2DCoordinates(m);
Assert.assertNotNull(factory.createTetrahedral(m.getAtom(1), Stereocenters.of(m)));
factory.withStrictMode();
Assert.assertNull(factory.createTetrahedral(m.getAtom(1), Stereocenters.of(m)));
Assert.assertNull(factory.createTetrahedral(m.getAtom(0), Stereocenters.of(m)));
}

@Test
Expand All @@ -650,10 +650,83 @@ public void okWedgePatternWithThreeNeighbors() throws CDKException {
Assert.assertNotNull(factory.createTetrahedral(m.getAtom(1), Stereocenters.of(m)));
}

// see http://efficientbits.blogspot.com/2019/09/rules-for-interpreting-updown-wedge.html
@Test
public void okWedgePatternWithThreeNeighbors2() throws Exception {
IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
try (InputStream in = getClass().getResourceAsStream("wedge_okay_d3.mol");
MDLV2000Reader mdlr = new MDLV2000Reader(in)) {
mdlr.getSetting("AddStereoElements").setSetting("false");
IAtomContainer mol = mdlr.read(builder.newAtomContainer());
int numStereo = 0;
for (IStereoElement<?,?> se : mol.stereoElements())
numStereo++;
assertEquals(0, numStereo);
StereoElementFactory stereoFactory = StereoElementFactory.using2DCoordinates(mol);
Assert.assertEquals(1, stereoFactory.createAll().size());
stereoFactory.withStrictMode();
Assert.assertEquals(1, stereoFactory.createAll().size());
}
}

@Test
public void badWedgePatternWithThreeNeighbors2() throws Exception {
IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
try (InputStream in = getClass().getResourceAsStream("wedge_bad_d3.mol");
MDLV2000Reader mdlr = new MDLV2000Reader(in)) {
mdlr.getSetting("AddStereoElements").setSetting("false");
IAtomContainer mol = mdlr.read(builder.newAtomContainer());
int numStereo = 0;
for (IStereoElement<?,?> se : mol.stereoElements())
numStereo++;
assertEquals(0, numStereo);
StereoElementFactory stereoFactory = StereoElementFactory.using2DCoordinates(mol);
Assert.assertEquals(1, stereoFactory.createAll().size());
stereoFactory.withStrictMode();
Assert.assertEquals(0, stereoFactory.createAll().size());
}
}

@Test
public void badWedgePatternWithThreeNeighbors180() throws Exception {
IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
try (InputStream in = getClass().getResourceAsStream("wedge_180_bad_d3.mol");
MDLV2000Reader mdlr = new MDLV2000Reader(in)) {
mdlr.getSetting("AddStereoElements").setSetting("false");
IAtomContainer mol = mdlr.read(builder.newAtomContainer());
int numStereo = 0;
for (IStereoElement<?,?> se : mol.stereoElements())
numStereo++;
assertEquals(0, numStereo);
StereoElementFactory stereoFactory = StereoElementFactory.using2DCoordinates(mol);
Assert.assertEquals(1, stereoFactory.createAll().size());
stereoFactory.withStrictMode();
Assert.assertEquals(0, stereoFactory.createAll().size());
}
}

@Test
public void okWedgePatternWithThreeNeighbors180() throws Exception {
IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
try (InputStream in = getClass().getResourceAsStream("wedge_180_okay_d3.mol");
MDLV2000Reader mdlr = new MDLV2000Reader(in)) {
mdlr.getSetting("AddStereoElements").setSetting("false");
IAtomContainer mol = mdlr.read(builder.newAtomContainer());
int numStereo = 0;
for (IStereoElement<?,?> se : mol.stereoElements())
numStereo++;
assertEquals(0, numStereo);
StereoElementFactory stereoFactory = StereoElementFactory.using2DCoordinates(mol);
Assert.assertEquals(1, stereoFactory.createAll().size());
stereoFactory.withStrictMode();
Assert.assertEquals(1, stereoFactory.createAll().size());
}
}

@Test
public void badWedgePatternWithFourNeighbors() throws Exception {
IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
try (InputStream in = getClass().getResourceAsStream("bad-wedges-4nbors.mol");
try (InputStream in = getClass().getResourceAsStream("wedge_bad_d4.mol");
MDLV2000Reader mdlr = new MDLV2000Reader(in)) {
mdlr.getSetting("AddStereoElements").setSetting("false");
IAtomContainer mol = mdlr.read(builder.newAtomContainer());
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@

ChemDraw02112215352D

4 3 0 0 0 0 0 0 0 0999 V2000
-0.8250 0.4094 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 0.4115 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.8250 0.4135 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
0.0021 -0.4135 0.0000 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
2 4 1 1
M END
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@

ChemDraw02112216072D

4 3 0 0 1 0 0 0 0 0999 V2000
-0.8250 0.4094 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.0000 0.4115 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.8250 0.4135 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
0.0021 -0.4135 0.0000 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 1
2 4 1 0
M END
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@

ChemDraw02112215332D

4 3 0 0 1 0 0 0 0 0999 V2000
-0.7145 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 0.4125 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.7145 0.0000 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 -0.4125 0.0000 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 1
2 4 1 1
M END
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@

ChemDraw02112215322D

4 3 0 0 1 0 0 0 0 0999 V2000
-0.7145 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 0.4125 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.7145 0.0000 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 -0.4125 0.0000 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 6
2 4 1 1
M END