Dear Greg,
Patch attached (to be applied in the GraphMol/ChemReactions folder).
Some explanation:
DaylightParser.cpp:
The changes introduce two new values (4,5) for an atom's molInversion
flag which specify absolute stereochemistry (i.e. for cases 8&9: no
stereochemistry specified on the reactant atom, but stereochemistry in
the product). Exceptions need to be added here for the cases where the
stereochemistry isn't properly defined.
Reaction.cpp:
Added 'invertedBonding' function which tests whether the bond order of
mapped atoms around a tetrahedral stereocenter are changed. This is
derived from existing code.
Added a chiralTemplateAtomsToCheck vector which records atoms in the
template which are chiral in the reactant or product, to be checked
later, once the product is fully built. Removed the existing checking
of template atoms as it didn't have enough information to work while
the product was being built.
Added a loop over those chiralTemplateAtomsToCheck to do the work
(relying upon the molInversion flag). Note that this loop is a copy
and hack of the loop above it (chiralAtomsToCheck) and I noticed as I
made the patch that the comments are not to the new loop.
I've changed the sense of some of the tests from !=CHI_UNSPECIFIED to
==CHI_TETRAHEDRAL_XX; there was a CHI_OTHER in some of the tests -
not sure if this needs to be supported.
Best wishes,
Richard
On Wed, Aug 17, 2011 at 6:32 AM, Greg Landrum <[email protected]> wrote:
> Dear Richard,
>
> On Tue, Aug 16, 2011 at 10:35 PM, Richard Cooper
> <[email protected]> wrote:
> <snip>
>>
>> I've made a patch and some tests against 2011_03_2 to implement all
>> these cases - but best to discuss here first whether the spec above is
>> right.
>
> <snip>
>
> Thanks for doing a patch!
>
> It's going to take me a while to work my way through this. Can you
> please send the patch so I have something to play with as I'm doing
> so?
>
> Best,
> -greg
>
diff -Naur DaylightParser.orig.cpp DaylightParser.cpp
--- DaylightParser.orig.cpp 2010-08-21 01:16:02.000000000 +0100
+++ DaylightParser.cpp 2011-08-16 21:46:34.927215097 +0100
@@ -47,7 +47,57 @@
prodIt!=rxn->endProductTemplates();++prodIt){
for(ROMol::AtomIterator prodAtomIt=(*prodIt)->beginAtoms();
prodAtomIt!=(*prodIt)->endAtoms();++prodAtomIt){
- if((*prodAtomIt)->getChiralTag()!=Atom::CHI_UNSPECIFIED &&
+
+ if ( (*prodAtomIt)->hasProp("molAtomMapNumber") ) {
+ int mapNum;
+ (*prodAtomIt)->getProp("molAtomMapNumber",mapNum);
+
+ for(MOL_SPTR_VECT::const_iterator
reactIt=rxn->beginReactantTemplates();
+ reactIt!=rxn->endReactantTemplates();++reactIt){
+
+ for(ROMol::AtomIterator reactAtomIt=(*reactIt)->beginAtoms();
+ reactAtomIt!=(*reactIt)->endAtoms();++reactAtomIt){
+
+ if((*reactAtomIt)->hasProp("molAtomMapNumber")) {
+ int reactMapNum;
+ (*reactAtomIt)->getProp("molAtomMapNumber",reactMapNum);
+ if(reactMapNum==mapNum){
+ // finally, in the bowels of the nesting, we get to some
actual
+ // work:
+ if(
(*reactAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CW ) {
+ if(
(*prodAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CW ) {
+ (*prodAtomIt)->setProp("molInversionFlag",2);
// Keep
+ } else if (
(*prodAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CCW ) {
+ (*prodAtomIt)->setProp("molInversionFlag",1);
// Invert
+ } else {
+ (*prodAtomIt)->setProp("molInversionFlag",3);
// Racemize
+ }
+ } else if
((*reactAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CCW ) {
+ if(
(*prodAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CW ) {
+ (*prodAtomIt)->setProp("molInversionFlag",1);
// Invert
+ } else if (
(*prodAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CCW ) {
+ (*prodAtomIt)->setProp("molInversionFlag",2);
// Keep
+ } else {
+ (*prodAtomIt)->setProp("molInversionFlag",3);
// Racemize
+ }
+ } else {
+ if(
(*prodAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CW ) {
+ (*prodAtomIt)->setProp("molInversionFlag",4);
// Add CW
+ } else if (
(*prodAtomIt)->getChiralTag()==Atom::CHI_TETRAHEDRAL_CCW ) {
+ (*prodAtomIt)->setProp("molInversionFlag",5);
// Add CCW
+ } else {
+ (*prodAtomIt)->setProp("molInversionFlag",2);
// Implicit keep
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+
+/* if((*prodAtomIt)->getChiralTag()!=Atom::CHI_UNSPECIFIED &&
(*prodAtomIt)->getChiralTag()!=Atom::CHI_OTHER &&
(*prodAtomIt)->hasProp("molAtomMapNumber")) {
int mapNum;
@@ -81,6 +131,8 @@
}
}
}
+*/
+
}
}
}
diff -Naur Reaction.orig.cpp Reaction.cpp
--- Reaction.orig.cpp 2010-12-20 05:00:23.000000000 +0000
+++ Reaction.cpp 2011-08-16 21:51:45.180753567 +0100
@@ -232,6 +232,27 @@
return RWMOL_SPTR(res);
} // end of initProduct()
+ bool invertedBonding(const Atom *reactantAtom, RWMOL_SPTR product,
std::map<unsigned int,unsigned int> reactProdAtomMap) {
+ Atom *productAtom =
product->getAtomWithIdx(reactProdAtomMap[reactantAtom->getIdx()]);
+
+ if ( productAtom->getDegree() < 3 ) {
+ return false;
+ }
+ INT_LIST newOrder;
+ ROMol::OEDGE_ITER beg,end;
+ boost::tie(beg,end) =
reactantAtom->getOwningMol().getAtomBonds(reactantAtom);
+ while(beg!=end){
+ const BOND_SPTR reactantBond=reactantAtom->getOwningMol()[*beg];
+ unsigned int
oAtomIdx=reactantBond->getOtherAtomIdx(reactantAtom->getIdx());
+ const Bond *productBond;
+ productBond=product->getBondBetweenAtoms(productAtom->getIdx(),
+
reactProdAtomMap[oAtomIdx]);
+ newOrder.push_back(productBond->getIdx());
+ ++beg;
+ }
+ return ( productAtom->getPerturbationOrder(newOrder)%2 != 0 ) ;
+ }
+
void addReactantAtomsAndBonds(const ChemicalReaction *rxn,
RWMOL_SPTR product,const ROMOL_SPTR
reactantSptr,
const MatchVectType &match,
@@ -249,6 +270,7 @@
std::map<unsigned int,unsigned int> reactProdAtomMap; // this maps atom
indices from reactant->product
std::vector<int> prodReactAtomMap(product->getNumAtoms(),-1); // this
maps atom indices from product->reactant
std::vector<const Atom *> chiralAtomsToCheck;
+ std::vector<const Atom *> chiralTemplateAtomsToCheck;
for(unsigned int i=0;i<match.size();i++){
const Atom
*templateAtom=reactantTemplate->getAtomWithIdx(match[i].first);
if(templateAtom->hasProp("molAtomMapNumber")){
@@ -336,34 +358,35 @@
// --------- --------- --------- --------- --------- ---------
// While we're here, set the stereochemistry
// FIX: this should be free-standing, not in this function.
- if(reactantAtom->getChiralTag()!=Atom::CHI_UNSPECIFIED &&
- reactantAtom->getChiralTag()!=Atom::CHI_OTHER &&
- productAtom->hasProp("molInversionFlag")){
- int flagVal;
- productAtom->getProp("molInversionFlag",flagVal);
- productAtom->setChiralTag(reactantAtom->getChiralTag());
- switch(flagVal){
- case 0:
- // FIX: should we clear the chirality or leave it alone? for now
we leave it alone
- //productAtom->setChiralTag(Atom::ChiralType::CHI_UNSPECIFIED);
- break;
- case 1:
- // inversion
- if(reactantAtom->getChiralTag()!=Atom::CHI_TETRAHEDRAL_CW &&
- reactantAtom->getChiralTag()!=Atom::CHI_TETRAHEDRAL_CCW){
- BOOST_LOG(rdWarningLog) << "unsupported chiral type on
reactant atom ignored\n";
- } else {
- productAtom->invertChirality();
- }
- break;
- case 2:
- // retention: do nothing
- break;
- default:
- BOOST_LOG(rdWarningLog) << "unrecognized chiral
inversion/retention flag on product atom ignored\n";
- }
+ // NB - RIC: Moved stereochemistry checks to after building of
product, otherwise sometimes required bonds are not present.
+
+ // Chiral atoms in the product template have possibly different
behaviour to those outside the template.
+ // 1a. [C@:1]>>[C@:1] keep stereochemistry of reactant atom
+ // 1b. [C@@:1]>>[C@@:1] keep stereochemistry of reactant atom
+ // 1c. [C:1]>>[C:1] implilcitly keep stereochemistry of reactant
atom (unless bond degree changes)
+ // 2. [C@:1]>>[C@@:1] invert stereochemistry
+ // 3. [C@:1]>>[C:1] explicitly remove stereochemistry (racemization)
+ // 4. [C:1]>>[C@:1] explicitly add or replace stereochemistry
(chiral catalyst) - requires 4 connected atoms (or 3+H).
+
+ if ((reactantAtom->getChiralTag()!=Atom::CHI_UNSPECIFIED
)||(productAtom->getChiralTag()!=Atom::CHI_UNSPECIFIED )) {
+ chiralTemplateAtomsToCheck.push_back(reactantAtom);
//deals later with cases 1a&b, and 2 (excludes case 3&4)
}
-
+
+ //If the reaction transformed atom is specified chiral, but the
reactant atom is unspecified, then set the product atom to unspecified.
+ // /e.g./
+ // >>> rxn = AllChem.ReactionFromSmarts( '[O:2][C@:1]>>[Cl:2][C@:1]')
# replace OH with Cl. Keep chirality of C, if present.
+ // >>> print
Chem.MolToSmiles(rxn.RunReactants((Chem.MolFromSmiles('C[C@](O)(Br)I'),))[0][0],True)
+ // C[C@](Cl)(Br)I
+ //
+ // >>> rxn = AllChem.ReactionFromSmarts( '[O:2][C@:1]>>[Cl:2][C@:1]')
# replace OH with Cl. Keep chirality of C, if present.
+ // >>> print
Chem.MolToSmiles(rxn.RunReactants((Chem.MolFromSmiles('CC(O)(Br)I'),))[0][0],True)
+ // CC(Cl)(Br)I
+ //Could cause problems with reactions that form chiral products?
+ if(productAtom->getChiralTag()!=Atom::CHI_UNSPECIFIED &&
+ (reactantAtom->getChiralTag()==Atom::CHI_UNSPECIFIED ||
reactantAtom->getChiralTag()==Atom::CHI_UNSPECIFIED)){
+ productAtom->setChiralTag(reactantAtom->getChiralTag());
+ }
+
// now traverse:
std::list< const Atom * > atomStack;
atomStack.push_back(reactantAtom);
@@ -491,6 +514,90 @@
} // end of loop over chiralAtomsToCheck
// ---------- ---------- ---------- ---------- ---------- ----------
+ // now we need to loop over atoms from the reactants that were chiral
and
+ // directly involved in the reaction in order to make sure their
chirality has
+ // been correctly set
+ for(std::vector<const Atom *>::const_iterator
atomIt=chiralTemplateAtomsToCheck.begin();
+ atomIt!=chiralTemplateAtomsToCheck.end();++atomIt){
+ const Atom *reactantAtom=*atomIt;
+ Atom
*productAtom=product->getAtomWithIdx(reactProdAtomMap[reactantAtom->getIdx()]);
+// CHECK_INVARIANT(reactantAtom->getChiralTag()!=Atom::CHI_UNSPECIFIED,
+// "missing atom chirality.");
+//
CHECK_INVARIANT(reactantAtom->getChiralTag()==productAtom->getChiralTag(),
+// "invalid product chirality.");
+
+ if( reactantAtom->getOwningMol().getAtomDegree(reactantAtom) !=
+ product->getAtomDegree(productAtom) ){
+ // If the number of bonds to the atom has changed in the course of
the
+ // reaction we're lost, so remove chirality.
+ // A word of explanation here: the atoms in the
chiralTemplateAtomsToCheck set are
+ // not explicitly mapped atoms of the reaction, so we really have
no idea what
+ // to do with this case. At the moment I'm not even really sure how
this
+ // could happen, but better safe than sorry.
+ productAtom->setChiralTag(Atom::CHI_UNSPECIFIED);
+ } else if( ( reactantAtom->getChiralTag()==Atom::CHI_TETRAHEDRAL_CW ||
+ reactantAtom->getChiralTag()==Atom::CHI_TETRAHEDRAL_CCW )
) {
+ int flagVal=2;
+ if ( productAtom->hasProp("molInversionFlag") ) {
+ productAtom->getProp("molInversionFlag",flagVal);
+ }
+ productAtom->setChiralTag(reactantAtom->getChiralTag());
+ switch(flagVal){
+ case 0:
+ // FIX: should we clear the chirality or leave it alone? for now
we leave it alone
+ //productAtom->setChiralTag(Atom::ChiralType::CHI_UNSPECIFIED);
+ break;
+ case 1:
+ // inversion
+ if ( ! invertedBonding(reactantAtom, product, reactProdAtomMap)
) {
+ productAtom->invertChirality();
+ }
+ break;
+ case 2:
+ {
+ if ( invertedBonding(reactantAtom, product, reactProdAtomMap) )
{
+ productAtom->invertChirality();
+ }
+ }
+ break;
+ case 3: //Racemize
+ {
+ productAtom->setChiralTag(Atom::CHI_UNSPECIFIED);
+ }
+ break;
+ case 4: //Force CW
+ productAtom->setChiralTag(Atom::CHI_TETRAHEDRAL_CW);
+ break;
+ case 5: //Force CCW
+ productAtom->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
+ break;
+ default:
+ BOOST_LOG(rdWarningLog) << "unrecognized chiral
inversion/retention flag on product atom ignored\n";
+ }
+ } else if( reactantAtom->hasProp("_CIPCode") ) { // Reactant has no
assigned chirality, but it chiral
+ int flagVal=0;
+ if ( productAtom->hasProp("molInversionFlag") ) {
+ productAtom->getProp("molInversionFlag",flagVal);
+ }
+ switch(flagVal){
+ case 0:
+ //No relevant flag set. Ignore.
+ break;
+ case 4: //Force CW
+ productAtom->setChiralTag(Atom::CHI_TETRAHEDRAL_CW);
+ break;
+ case 5: //Force CCW
+ productAtom->setChiralTag(Atom::CHI_TETRAHEDRAL_CCW);
+ break;
+ default:
+ BOOST_LOG(rdWarningLog) << "unrecognized chiral
inversion/retention flag on product atom ignored\n";
+ }
+ }
+ } // end of loop over chiralTemplateAtomsToCheck
+
+
+
+ // ---------- ---------- ---------- ---------- ---------- ----------
// finally we may need to set the coordinates in the product conformer:
if(productConf){
productConf->resize(product->getNumAtoms());
@@ -568,13 +675,15 @@
return productMols;
}
+
// find the matches for each reactant:
VectVectMatchVectType matchesByReactant;
if(!ReactionUtils::getReactantMatches(reactants,this->m_reactantTemplates,matchesByReactant)){
// some reactants didn't find a match, return an empty product list:
return productMols;
}
-
+
+
// -------------------------------------------------------
// we now have matches for each reactant, so we can start creating
products:
------------------------------------------------------------------------------
Get a FREE DOWNLOAD! and learn more about uberSVN rich system,
user administration capabilities and model configuration. Take
the hassle out of deploying and managing Subversion and the
tools developers use with it. http://p.sf.net/sfu/wandisco-d2d-2
_______________________________________________
Rdkit-discuss mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/rdkit-discuss
