Skip to content

Commit

Permalink
Reading of 0D stereochemistry from the CFG=1 and CFG=2 values in MDLV…
Browse files Browse the repository at this point in the history 
…3000 to be more consistent with V2000. This is a corner case and I may change to off by default.
  • Loading branch information
@johnmay @egonw
johnmay authored and egonw committed Feb 7, 2022
1 parent 3f55fb3 commit 23cec7e
Show file tree
Hide file tree
Showing 4 changed files with 150 additions and 80 deletions.
71 changes: 39 additions & 32 deletions storage/ctab/src/main/java/org/openscience/cdk/io/MDLV2000Reader.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
import org.openscience.cdk.interfaces.IChemObject;
import org.openscience.cdk.interfaces.IChemObjectBuilder;
import org.openscience.cdk.interfaces.IChemSequence;
import org.openscience.cdk.interfaces.IElement;
import org.openscience.cdk.interfaces.IIsotope;
import org.openscience.cdk.interfaces.IPseudoAtom;
import org.openscience.cdk.interfaces.ISingleElectron;
Expand Down Expand Up @@ -457,39 +458,11 @@ private IAtomContainer readAtomContainer(IAtomContainer molecule) throws CDKExce

// create 0D stereochemistry
if (optStereoPerc.isSet() && optStereo0d.isSet()) {
Parities:
for (Map.Entry<IAtom, Integer> e : parities.entrySet()) {
int parity = e.getValue();
if (parity != 1 && parity != 2)
continue; // 3=unspec
int idx = 0;
IAtom focus = e.getKey();
IAtom[] carriers = new IAtom[4];
int hidx = -1;
for (IAtom nbr : outputContainer.getConnectedAtomsList(focus)) {
if (idx == 4)
continue Parities; // too many neighbors
if (nbr.getAtomicNumber() == 1) {
if (hidx >= 0)
continue Parities;
hidx = idx;
}
carriers[idx++] = nbr;
}
// to few neighbors, or already have a hydrogen defined
if (idx < 3 || idx < 4 && hidx >= 0)
continue;
if (idx == 3)
carriers[idx++] = focus;

if (idx == 4) {
Stereo winding = parity == 1 ? Stereo.CLOCKWISE : Stereo.ANTI_CLOCKWISE;
// H is always at back, even if explicit! At least this seems to be the case.
// we adjust the winding as needed
if (hidx == 0 || hidx == 2)
winding = winding.invert();
outputContainer.addStereoElement(new TetrahedralChirality(focus, carriers, winding));
}
IStereoElement<IAtom,IAtom> stereoElement
= createStereo0d(outputContainer, e.getKey(), e.getValue());
if (stereoElement != null)
molecule.addStereoElement(stereoElement);
}
}

Expand Down Expand Up @@ -560,6 +533,40 @@ private IAtomContainer readAtomContainer(IAtomContainer molecule) throws CDKExce
return outputContainer;
}

static IStereoElement<IAtom, IAtom> createStereo0d(IAtomContainer mol, IAtom focus, int parity) {
if (parity != 1 && parity != 2)
return null; // 3=unspec
int numNbrs = 0;
IAtom[] carriers = new IAtom[4];
int idxOfHyd = -1;
for (IAtom nbr : mol.getConnectedAtomsList(focus)) {
if (numNbrs == 4)
return null; // too many neighbors
if (nbr.getAtomicNumber() == IElement.H) {
if (idxOfHyd >= 0)
return null; // too many hydrogens
idxOfHyd = numNbrs;
}
carriers[numNbrs++] = nbr;
}
// incorrect number of neighbours?
if (numNbrs < 3 || numNbrs < 4 && idxOfHyd >= 0)
return null;
// implicit neighbour (H or lone-pair)
if (numNbrs == 3)
carriers[numNbrs++] = focus;
if (numNbrs != 4)
return null;

Stereo winding = parity == 1 ? Stereo.CLOCKWISE : Stereo.ANTI_CLOCKWISE;
// H is always at back, even if explicit! At least this seems to be the case.
// we adjust the winding as needed which is when the explict H is in slot
// 0 or 2 (odd number of swaps to get to index 3)
if (idxOfHyd == 0 || idxOfHyd == 2)
winding = winding.invert();
return new TetrahedralChirality(focus, carriers, winding);
}

private boolean is3Dfile(String program) {
return program.length() >= 22 && program.substring(20, 22).equals("3D");
}
Expand Down
104 changes: 56 additions & 48 deletions storage/ctab/src/main/java/org/openscience/cdk/io/MDLV3000Reader.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,60 +231,67 @@ private void finalizeMol(ReadState state) {
}
}

if (!isQuery) {
if (!isQuery)
finalizeStereochemistry(state, readData);
}

/* Add stereo information */
boolean is3d = true;
for (IAtom atom : readData.atoms()) {
if (atom.getPoint3d() == null)
is3d = false;
}
if (optStereoPerc.isSet()) {
if (is3d) { // has 3D coordinates
readData.setStereoElements(StereoElementFactory.using3DCoordinates(readData)
.createAll());
} else { // has 2D coordinates (set as 2D coordinates)
readData.setStereoElements(StereoElementFactory.using2DCoordinates(readData)
.createAll());
private void finalizeStereochemistry(ReadState state, IAtomContainer readData) {
if (optStereoPerc.isSet()) {

if (state.dimensions == 3) { // has 3D coordinates
readData.setStereoElements(StereoElementFactory.using3DCoordinates(readData)
.createAll());
} else if (state.dimensions == 2) { // has 2D coordinates (set as 2D coordinates)
readData.setStereoElements(StereoElementFactory.using2DCoordinates(readData)
.createAll());
} else if (state.dimensions == 0 && optStereo0d.isSet()) {
// technically if a molecule is 2D/3D and has the CFG=1 or CFG=2
// specified this gives us hints information but it's safer to
// just use the coordinates or wedge bonds
for (Map.Entry<IAtom, Integer> e : state.stereo0d.entrySet()) {
IStereoElement<IAtom,IAtom> stereoElement
= MDLV2000Reader.createStereo0d(state.mol, e.getKey(), e.getValue());
if (stereoElement != null)
state.mol.addStereoElement(stereoElement);
}
}

if (state.stereoflags != null && !state.stereoflags.isEmpty()) {

// work out the next available group, if we have &1, &2, etc then we choose &3
// this is only needed if
int defaultRacGrp = 0;
if (!state.chiral) {
int max = 0;
for (Integer val : state.stereoflags.values()) {
if ((val & IStereoElement.GRP_TYPE_MASK) == IStereoElement.GRP_RAC) {
int num = val >>> IStereoElement.GRP_NUM_SHIFT;
if (num > max)
max = num;
}
if (state.stereoflags != null && !state.stereoflags.isEmpty()) {

// work out the next available group, if we have &1, &2, etc then we choose &3
// this is only needed if
int defaultRacGrp = 0;
if (!state.chiral) {
int max = 0;
for (Integer val : state.stereoflags.values()) {
if ((val & IStereoElement.GRP_TYPE_MASK) == IStereoElement.GRP_RAC) {
int num = val >>> IStereoElement.GRP_NUM_SHIFT;
if (num > max)
max = num;
}
defaultRacGrp = IStereoElement.GRP_RAC | (((max + 1) << IStereoElement.GRP_NUM_SHIFT));
}
defaultRacGrp = IStereoElement.GRP_RAC | (((max + 1) << IStereoElement.GRP_NUM_SHIFT));
}

for (IStereoElement<?, ?> se : readData.stereoElements()) {
if (se.getConfigClass() != IStereoElement.TH)
continue;
IAtom focus = (IAtom) se.getFocus();
int idx = readData.indexOf(focus);
if (idx < 0)
continue;
Integer grpinfo = state.stereoflags.get(idx);
if (grpinfo != null)
se.setGroupInfo(grpinfo);
else if (!state.chiral)
se.setGroupInfo(defaultRacGrp);
}
} else if (!state.chiral) {
// chiral flag not set which means this molecule is this stereoisomer "and" the enantiomer, mark all
// Tetrahedral stereo as AND1 (&1)
for (IStereoElement<?, ?> se : readData.stereoElements()) {
if (se.getConfigClass() == IStereoElement.TH) {
se.setGroupInfo(IStereoElement.GRP_RAC1);
}
for (IStereoElement<?, ?> se : readData.stereoElements()) {
if (se.getConfigClass() != IStereoElement.TH)
continue;
IAtom focus = (IAtom) se.getFocus();
int idx = readData.indexOf(focus);
if (idx < 0)
continue;
Integer grpinfo = state.stereoflags.get(idx);
if (grpinfo != null)
se.setGroupInfo(grpinfo);
else if (!state.chiral)
se.setGroupInfo(defaultRacGrp);
}
} else if (!state.chiral) {
// chiral flag not set which means this molecule is this stereoisomer "and" the enantiomer, mark all
// Tetrahedral stereo as AND1 (&1)
for (IStereoElement<?, ?> se : readData.stereoElements()) {
if (se.getConfigClass() == IStereoElement.TH) {
se.setGroupInfo(IStereoElement.GRP_RAC1);
}
}
}
Expand All @@ -311,6 +318,7 @@ private void finalizeDimensions(ReadState state) {
// check the global header
if (dimensions == 0)
dimensions = state.dimensions;
state.dimensions = dimensions;

if (dimensions == 0) {
// remove all coords we set
Expand Down
27 changes: 27 additions & 0 deletions storage/ctab/src/test/java/org/openscience/cdk/io/MDLV3000ReaderTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,4 +401,31 @@ public void testStereoRel1() throws Exception {
}
}
}

@Test
public void testStereo0d() throws Exception {
IChemObjectBuilder bldr = SilentChemObjectBuilder.getInstance();
try (InputStream in = getClass().getResourceAsStream("stereo0d.mdl3");
MDLV3000Reader mdlr = new MDLV3000Reader(in)) {
IAtomContainer mol = mdlr.read(bldr.newAtomContainer());
int numTetrahedrals = 0;
for (IStereoElement<?,?> se : mol.stereoElements()) {
if (se.getConfigClass() == IStereoElement.TH)
numTetrahedrals++;
}
Assert.assertEquals(2, numTetrahedrals);
}

try (InputStream in = getClass().getResourceAsStream("stereo0d.mdl3");
MDLV3000Reader mdlr = new MDLV3000Reader(in)) {
mdlr.getSetting("AddStereo0d").setSetting("false");
IAtomContainer mol = mdlr.read(bldr.newAtomContainer());
int numTetrahedrals = 0;
for (IStereoElement<?,?> se : mol.stereoElements()) {
if (se.getConfigClass() == IStereoElement.TH)
numTetrahedrals++;
}
Assert.assertEquals(0, numTetrahedrals);
}
}
}
28 changes: 28 additions & 0 deletions storage/ctab/src/test/resources/org/openscience/cdk/io/stereo0d.mdl3
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@

ChemDraw0207221100

0 0 0 0 0 0 V3000
M V30 BEGIN CTAB
M V30 COUNTS 8 8 0 0 0
M V30 BEGIN ATOM
M V30 1 C 0.0 0.0 0.0 0
M V30 2 N 0.0 0.0 0.0 0
M V30 3 C 0.0 0.0 0.0 0
M V30 4 C 0.0 0.0 0.0 0
M V30 5 C 0.0 0.0 0.0 0 CFG=1
M V30 6 C 0.0 0.0 0.0 0 CFG=1
M V30 7 O 0.0 0.0 0.0 0
M V30 8 O 0.0 0.0 0.0 0
M V30 END ATOM
M V30 BEGIN BOND
M V30 1 2 1 2
M V30 2 1 2 3
M V30 3 1 3 4
M V30 4 1 4 5
M V30 5 1 5 6
M V30 6 1 6 1
M V30 7 1 6 7
M V30 8 1 5 8
M V30 END BOND
M V30 END CTAB
M END

0 comments on commit 23cec7e

Please sign in to comment.