C++ MatchVectType类代码示例

double AlignMolecule(ROMol &prbMol, const ROMol &refMol, int prbCid = -1,
                     int refCid = -1, python::object atomMap = python::list(),
                     python::object weights = python::list(),
                     bool reflect = false, unsigned int maxIters = 50) {
    MatchVectType *aMap = _translateAtomMap(atomMap);
    unsigned int nAtms;
    if (aMap) {
        nAtms = aMap->size();
    } else {
        nAtms = prbMol.getNumAtoms();
    }
    RDNumeric::DoubleVector *wtsVec = _translateWeights(weights);
    if (wtsVec) {
        if (wtsVec->size() != nAtms) {
            throw_value_error("Incorrect number of weights specified");
        }
    }

    double rmsd;
    {
        NOGIL gil;
        rmsd = MolAlign::alignMol(prbMol, refMol, prbCid, refCid, aMap, wtsVec,
                                  reflect, maxIters);
    }
    if (aMap) {
        delete aMap;
    }
    if (wtsVec) {
        delete wtsVec;
    }
    return rmsd;
}

void testMMFFO3AConstraintsAndLocalOnly() {
  std::string rdbase = getenv("RDBASE");
  std::string sdf = rdbase + "/Code/GraphMol/MolAlign/test_data/ref_e2.sdf";
  SDMolSupplier supplier(sdf, true, false);
  const int refNum = 23;
  const int prbNum = 32;
  ROMol *refMol = supplier[refNum];
  ROMol *prbMol = supplier[prbNum];
  unsigned int refNAtoms = refMol->getNumAtoms();
  std::vector<double> refLogpContribs(refNAtoms);
  std::vector<double> refMRContribs(refNAtoms);
  std::vector<unsigned int> refAtomTypes(refNAtoms);
  std::vector<std::string> refAtomTypeLabels(refNAtoms);
  Descriptors::getCrippenAtomContribs(*refMol, refLogpContribs,
    refMRContribs, true, &refAtomTypes, &refAtomTypeLabels);
  unsigned int prbNAtoms = prbMol->getNumAtoms();
  std::vector<double> prbLogpContribs(prbNAtoms);
  std::vector<double> prbMRContribs(prbNAtoms);
  std::vector<unsigned int> prbAtomTypes(prbNAtoms);
  std::vector<std::string> prbAtomTypeLabels(prbNAtoms);
  Descriptors::getCrippenAtomContribs(*prbMol, prbLogpContribs,
    prbMRContribs, true, &prbAtomTypes, &prbAtomTypeLabels);
  RWMol *patt = SmartsToMol("S");
  MatchVectType matchVect;
  TEST_ASSERT(SubstructMatch(*refMol, (ROMol &)*patt, matchVect));
  delete patt;
  unsigned int refSIdx = matchVect[0].second;
  matchVect.clear();
  patt = SmartsToMol("O");
  TEST_ASSERT(SubstructMatch(*prbMol, (ROMol &)*patt, matchVect));
  delete patt;
  unsigned int prbOIdx = matchVect[0].second;
  std::vector<double> distOS(2);
  distOS[0] = 2.7;
  distOS[1] = 0.4;
  std::vector<double> weights(2);
  weights[0] = 0.1;
  weights[1] = 100.0;
  for (unsigned int i = 0; i < 2; ++i) {
    MatchVectType constraintMap;
    constraintMap.push_back(std::make_pair(prbOIdx, refSIdx));
    RDNumeric::DoubleVector constraintWeights(1);
    constraintWeights[0] = weights[i];
    MolAlign::O3A *o3a = new MolAlign::O3A(*prbMol, *refMol,
      &prbLogpContribs, &refLogpContribs, MolAlign::O3A::CRIPPEN,
      -1, -1, false, 50, 0, &constraintMap, &constraintWeights);
    TEST_ASSERT(o3a);
    o3a->align();
    delete o3a;
    o3a = new MolAlign::O3A(*prbMol, *refMol, &prbLogpContribs,
      &refLogpContribs, MolAlign::O3A::CRIPPEN, -1, -1,
      false, 50, MolAlign::O3_LOCAL_ONLY);
    TEST_ASSERT(o3a);
    o3a->align();
    delete o3a;
    double d = (prbMol->getConformer().getAtomPos(prbOIdx)
      - refMol->getConformer().getAtomPos(refSIdx)).length();
    TEST_ASSERT(feq(d, distOS[i], 0.1));
  }
}

void testDativeMatch() {
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "    Test dative-bond matching" << std::endl;
  {
    std::string smi = "[Cu]->[Fe]";
    ROMol *mol = SmilesToMol(smi);
    TEST_ASSERT(mol);

    // make sure a self-match works
    MatchVectType match;
    TEST_ASSERT(SubstructMatch(*mol, *mol, match));
    TEST_ASSERT(match.size() == mol->getNumAtoms());

    {  // reverse the order and make sure that works
      std::string sma = "[Fe]<-[Cu]";
      ROMol *qmol = SmilesToMol(sma);
      TEST_ASSERT(qmol);
      MatchVectType match;
      TEST_ASSERT(SubstructMatch(*mol, *qmol, match));
      TEST_ASSERT(match.size() == qmol->getNumAtoms());
      delete qmol;
    }
    {  // reverse the direction and make sure that does not work.
      std::string sma = "[Fe]->[Cu]";
      ROMol *qmol = SmilesToMol(sma);
      TEST_ASSERT(qmol);
      MatchVectType match;
      TEST_ASSERT(!SubstructMatch(*mol, *qmol, match));
      delete qmol;
    }

    delete mol;
  }
  BOOST_LOG(rdErrorLog) << "  done" << std::endl;
}

PyObject *getMolAlignTransform(const ROMol &prbMol, const ROMol &refMol,
                               int prbCid = -1, int refCid = -1,
                               python::object atomMap = python::list(),
                               python::object weights = python::list(),
                               bool reflect = false,
                               unsigned int maxIters = 50) {
    MatchVectType *aMap = _translateAtomMap(atomMap);
    unsigned int nAtms;
    if (aMap) {
        nAtms = aMap->size();
    } else {
        nAtms = prbMol.getNumAtoms();
    }
    RDNumeric::DoubleVector *wtsVec = _translateWeights(weights);
    if (wtsVec) {
        if (wtsVec->size() != nAtms) {
            throw_value_error("Incorrect number of weights specified");
        }
    }
    RDGeom::Transform3D trans;
    double rmsd;
    {
        NOGIL gil;
        rmsd = MolAlign::getAlignmentTransform(
                   prbMol, refMol, trans, prbCid, refCid, aMap, wtsVec, reflect, maxIters);
    }
    if (aMap) {
        delete aMap;
    }
    if (wtsVec) {
        delete wtsVec;
    }

    return generateRmsdTransPyTuple(rmsd, trans);
}

void test4() {
  std::cout << " ----------------- Test 4" << std::endl;
  MatchVectType matchV;
  std::vector<MatchVectType> matches;
  int n;
  bool updateLabel = true;
  bool takeOwnership = true;

  RWMol *m, *q1, *q2;
  auto *a6 = new Atom(6);
  auto *a8 = new Atom(8);
  m = new RWMol();
  m->addAtom(a6);
  m->addAtom(a6);
  m->addAtom(a8);
  m->addAtom(a6);
  m->addAtom(a6);
  m->addBond(1, 0, Bond::SINGLE);
  m->addBond(1, 2, Bond::SINGLE);
  m->addBond(1, 3, Bond::SINGLE);
  m->addBond(2, 4, Bond::SINGLE);

  // this will be the recursive query
  q1 = new RWMol();
  q1->addAtom(new QueryAtom(6), updateLabel, takeOwnership);
  q1->addAtom(new QueryAtom(8), updateLabel, takeOwnership);
  q1->addBond(0, 1, Bond::UNSPECIFIED);

  // here's the main query
  q2 = new RWMol();
  auto *qA = new QueryAtom(6);
  auto *rsq = new RecursiveStructureQuery(q1);
  qA->expandQuery(rsq, Queries::COMPOSITE_AND);
  // std::cout << "post expand: " << qA->getQuery() << std::endl;
  q2->addAtom(qA, true, true);
  // std::cout << "mol: " << q2->getAtomWithIdx(0)->getQuery() << std::endl;
  q2->addAtom(new QueryAtom(6), true, true);
  q2->addBond(0, 1, Bond::UNSPECIFIED);

  bool found = SubstructMatch(*m, *q2, matchV);
  CHECK_INVARIANT(found, "");
  CHECK_INVARIANT(matchV.size() == 2, "");
  TEST_ASSERT(matchV[0].first == 0);
  TEST_ASSERT(matchV[0].second == 1);
  TEST_ASSERT(matchV[1].first == 1);
  TEST_ASSERT(matchV[1].second == 0 || matchV[1].second == 3);
  n = SubstructMatch(*m, *q2, matches, true);
  TEST_ASSERT(n == 2);
  TEST_ASSERT(matches.size() == (size_t)n);
  TEST_ASSERT(matches[0].size() == 2);
  TEST_ASSERT(matches[1].size() == 2);
  TEST_ASSERT(matches[0][0].second == matches[1][0].second);
  TEST_ASSERT(matches[0][1].second != matches[1][1].second);
  delete m;
  delete a6;
  delete a8;
  delete q2;

  std::cout << "Done\n" << std::endl;
}

void testSubstructMatchDMAP() {
  BOOST_LOG(rdInfoLog) << "-----------------------\n"
                       << "testSubstructMatchDMAP" << std::endl;
  RWMol *mol = SmilesToMol("C(C)Nc1cc[nH+]cc1");
  RWMol *query = SmartsToMol("[#7+]");
  ResonanceMolSupplier *resMolSuppl;
  unsigned int n;
  MatchVectType p;

  resMolSuppl = new ResonanceMolSupplier((ROMol &)*mol);
  std::vector<MatchVectType> matchVect;
  n = SubstructMatch(*mol, *query, matchVect, false, true, false, false);
  TEST_ASSERT((n == 1) && (matchVect.size() == 1));
  p = matchVect[0];
  TEST_ASSERT((p.size() == 1) && (p[0].second == 6));
  matchVect.clear();
  n = SubstructMatch(*resMolSuppl, *query, matchVect, false, true, false,
                     false);
  TEST_ASSERT((n == 2) && (matchVect.size() == 2));
  std::vector<unsigned int> v;
  p = matchVect[0];
  TEST_ASSERT(p.size() == 1);
  v.push_back(p[0].second);
  p = matchVect[1];
  TEST_ASSERT(p.size() == 1);
  v.push_back(p[0].second);
  std::sort(v.begin(), v.end());
  TEST_ASSERT(v[0] == 2);
  TEST_ASSERT(v[1] == 6);
  delete mol;
  delete query;
  delete resMolSuppl;
}

 double getAlignmentTransform(const ROMol &prbMol, const ROMol &refMol,
                                      RDGeom::Transform3D &trans,
                                      int prbCid, int refCid, 
                                      const MatchVectType *atomMap, 
                                      const RDNumeric::DoubleVector *weights, bool reflect,
                                      unsigned int maxIterations) {
   RDGeom::Point3DConstPtrVect refPoints, prbPoints;
   const Conformer &prbCnf = prbMol.getConformer(prbCid);
   const Conformer &refCnf = refMol.getConformer(refCid);
   if (atomMap == 0) {
     // we have to figure out the mapping between the two molecule
     MatchVectType match;
     if (SubstructMatch(refMol, prbMol, match)) {
       MatchVectType::const_iterator mi;
       for (mi = match.begin(); mi != match.end(); mi++) {
         prbPoints.push_back(&prbCnf.getAtomPos(mi->first));
         refPoints.push_back(&refCnf.getAtomPos(mi->second));
       }
     } else {
       throw MolAlignException("No sub-structure match found between the probe and query mol");
     } 
   } else {
     MatchVectType::const_iterator mi;
     for (mi = atomMap->begin(); mi != atomMap->end(); mi++) {
       prbPoints.push_back(&prbCnf.getAtomPos(mi->first));
       refPoints.push_back(&refCnf.getAtomPos(mi->second));
     }
   }
   double ssr = RDNumeric::Alignments::AlignPoints(refPoints, prbPoints, 
                                                   trans, weights, reflect, maxIterations);
   ssr /= (prbPoints.size());
   return sqrt(ssr);
 }

  ROMol *prepareMol(const ROMol &mol, const FragCatParams *fparams,
		    MatchVectType &aToFmap) {
    PRECONDITION(fparams,"");
    
    // get a mapping of the functional groups onto the molecule
    INT_VECT fgBonds;
    MatchVectType aidToFid = findFuncGroupsOnMol(mol, fparams, fgBonds);
    
    // get the core piece of molecule (i.e. the part of the molecule 
    // without the functional groups). This basically the part of the molecule
    // that does not contain the function group bonds given by "fgBonds"
    INT_VECT cBonds;
    int bid, nbds = mol.getNumBonds();

    for (bid = 0; bid < nbds; bid++) {
      if (std::find(fgBonds.begin(), fgBonds.end(), bid) == fgBonds.end()) {
	cBonds.push_back(bid);
      }
    }

    INT_MAP_INT aIdxMap; // a map from atom id in mol to the new atoms id in coreMol
    
    ROMol *coreMol = Subgraphs::pathToSubmol(mol, cBonds, false, aIdxMap);
    
    // now map the functional groups on mol to coreMol using aIdxMap
    MatchVectType::iterator mati;
    
    int newID;
    for (mati = aidToFid.begin(); mati != aidToFid.end(); mati++) {
      newID = aIdxMap[mati->first];
      aToFmap.push_back(std::pair<int, int>(newID, mati->second));
    }
    
    return coreMol;
  }

  //*************************************************************************************
  //
  //  Adds 2D coordinates to a molecule using the Depict.dll
  //
  //  ARGUMENTS:
  //    mol:          the molecule to be altered
  //    tempFilename: (OPTIONAL) the name of the temporary file 
  //
  //  RETURNS:
  //   1 on success, 0 otherwise
  //
  //  Here's the process by which this works (it's kind of contorted):
  //   1) convert the mol to SMILES
  //   2) use the DLL to convert the SMILES to a mol file (on disk)
  //   3) parse the mol file into a temporary molecule
  //   4) do a substructure match from the old molecule to the
  //      temp one (which may have a different atom numbering or additional
  //      atoms added).
  //   5) Update the positions of the atoms on the old molecule.
  //   6) Free the temp molecule.
  //
  //  NOTES:
  //   - *FIX:* at the moment we're not doing anything to clear up the
  //     temp file created in this process.  It'll always have the same
  //     name unless the user explicitly asks that we do something different.
  //   - To use the DLL, it's essential that the COMBICHEM_ROOT and COMBICHEM_RELEASE
  //     environment variables be set.  If this isn't done, this whole process
  //     will fail.
  //   - See the notes above about failures when opening the DLL.
  //
  //*************************************************************************************
  int Add2DCoordsToMolDLL(ROMol &mol,std::string tempFilename){
    std::string smi=MolToSmiles(mol,true);
    int tmp = SmilesToMolFileDLL(smi,tempFilename);
    int res = -1;
    if(tmp){
      // build another mol from that mol file:
      RWMol *tmpMol = MolFileToMol(tempFilename,false);
      // match it up with the starting mol:
      //  (We need to do this because the depict.dll conversion
      //   to a mol file may have added Hs)
      MatchVectType matchVect;
      bool hasMatch=SubstructMatch(tmpMol,&mol,matchVect);
      if(hasMatch){
        const Conformer &conf = tmpMol->getCoformer(0);
        Coformer *nconf = new Coformer(mol.getNumAtoms());
	for(MatchVectType::const_iterator mvi=matchVect.begin();
	    mvi!=matchVect.end();mvi++){
	  
          nconf->setAtomPos(conf.getAtomPos(mvi->first));
	}
        confId = (int)mol.addConformer(nconf, true);
      }
      delete tmpMol;
    }
    return res;
  }

bool StructCheckTautomer::applyTautomer(unsigned it) {
  if (Options.FromTautomer.size() <= it || Options.ToTautomer.size() <= it) {
    if (Options.Verbose)
      BOOST_LOG(rdInfoLog) << "ERROR: incorrect Tautomer index it=" << it
                           << "\n";
    return false;
  }
  const ROMol &fromTautomer = *Options.FromTautomer[it];
  const ROMol &toTautomer = *Options.ToTautomer[it];
  if (toTautomer.getNumAtoms() != fromTautomer.getNumAtoms()) {
    if (Options.Verbose)
      BOOST_LOG(rdInfoLog) << "ERROR: incorrect data toTautomer.getNumAtoms() "
                              "!= fromTautomer.getNumAtoms()\n";
    // incorrect data
    // throw(.....);
    return false;
  }
  const unsigned nta = toTautomer.getNumAtoms();
  const unsigned ntb = toTautomer.getNumBonds();
  MatchVectType match;  // The format is (queryAtomIdx, molAtomIdx)
  std::vector<unsigned> atomIdxMap(
      Mol.getNumAtoms());  // matched tau atom indeces

  if (!SubstructMatch(Mol, *Options.FromTautomer[it],
                      match))  // SSMatch(mp, from_tautomer, SINGLE_MATCH);
    return false;
  if (Options.Verbose)
    BOOST_LOG(rdInfoLog) << "found match for from_tautomer with " << nta
                         << " atoms\n";
  // init
  for (unsigned i = 0; i < Mol.getNumAtoms(); i++) atomIdxMap[i] = -1;
  for (MatchVectType::const_iterator mit = match.begin(); mit != match.end();
       mit++) {
    unsigned tai = mit->first;   // From and To Tautomer Atom index
    unsigned mai = mit->second;  // Mol Atom index
    atomIdxMap[mai] = tai;
  }
  // scan for completely mapped bonds and replace bond order with mapped bond
  // from to_tautomer
  for (RDKit::BondIterator_ bond = Mol.beginBonds(); bond != Mol.endBonds();
       bond++) {
    unsigned ti = atomIdxMap[(*bond)->getBeginAtomIdx()];
    unsigned tj = atomIdxMap[(*bond)->getEndAtomIdx()];
    if (-1 == ti || -1 == tj) continue;
    const Bond *tb = toTautomer.getBondBetweenAtoms(ti, tj);
    if (tb && (*bond)->getBondType() != tb->getBondType()) {
      (*bond)->setBondType(tb->getBondType());
    }
  }
  // apply charge/radical fixes if any
  for (unsigned i = 0; i < match.size(); i++) {
    Atom &atom = *Mol.getAtomWithIdx(match[i].second);
    const Atom &ta = *toTautomer.getAtomWithIdx(match[i].first);
    atom.setFormalCharge(ta.getFormalCharge());
    atom.setNumRadicalElectrons(ta.getNumRadicalElectrons());
  }

  return true;
}

void test2(){
  std::cout << " ----------------- Test 2" << std::endl;
  MatchVectType matchV;
  std::vector< MatchVectType > matches;
  unsigned int n;

  RWMol *m,*q1;
  m = new RWMol();
  m->addAtom(new Atom(6));
  m->addAtom(new Atom(6));
  m->addAtom(new Atom(8));
  m->addBond(0,1,Bond::SINGLE);
  m->addBond(1,2,Bond::SINGLE);

  q1 = new RWMol();
  q1->addAtom(new QueryAtom(6));
  q1->addAtom(new QueryAtom(8));
  q1->addBond(0,1,Bond::SINGLE);

  n = SubstructMatch(*m,*q1,matchV);
  TEST_ASSERT(n);
  TEST_ASSERT(matchV.size()==2);
  TEST_ASSERT(matchV[0].first==0);
  TEST_ASSERT(matchV[0].second==1);
  TEST_ASSERT(matchV[1].first==1);
  TEST_ASSERT(matchV[1].second==2);

  n = SubstructMatch(*m,*q1,matches,false);
  CHECK_INVARIANT(n==1,"");
  CHECK_INVARIANT(matches.size()==n,"");
  CHECK_INVARIANT(matches[0].size()==2,"");
  n = SubstructMatch(*m,*q1,matches,true);
  CHECK_INVARIANT(n==1,"");
  CHECK_INVARIANT(matches.size()==n,"");
  CHECK_INVARIANT(matches[0].size()==2,"");

  CHECK_INVARIANT(SubstructMatch(*m,*q1,matchV),"");
  CHECK_INVARIANT(matchV.size()==2,"");

  delete m;
  m = new RWMol();
  m->addAtom(new Atom(6));
  m->addAtom(new Atom(6));
  m->addAtom(new Atom(8));
  m->addBond(0,1,Bond::SINGLE);
  m->addBond(1,2,Bond::DOUBLE);

  matches.clear();
  n = SubstructMatch(*m,*q1,matches,false);
  CHECK_INVARIANT(n==0,"");
  CHECK_INVARIANT(matches.size()==n,"");
  n = SubstructMatch(*m,*q1,matches,true);
  CHECK_INVARIANT(n==0,"");
  CHECK_INVARIANT(matches.size()==n,"");
  CHECK_INVARIANT(!SubstructMatch(*m,*q1,matchV),"");

  std::cout << "Done\n" << std::endl;

}

void testRecursiveSerialNumbers(){
  std::cout << " ----------------- Testing serial numbers on recursive queries" << std::endl;
  MatchVectType matchV;
  std::vector< MatchVectType > matches;
  int n;

  RWMol *m,*q1,*q2;
  Atom *a6 = new Atom(6);
  Atom *a8 = new Atom(8);
  m = new RWMol();
  m->addAtom(a6);
  m->addAtom(a6);
  m->addAtom(a8);
  m->addAtom(a6);
  m->addAtom(a6);
  m->addBond(1,0,Bond::SINGLE);
  m->addBond(1,2,Bond::SINGLE);
  m->addBond(1,3,Bond::SINGLE);
  m->addBond(2,4,Bond::SINGLE);

  {
    // this will be the recursive query
    q1 = new RWMol();
    q1->addAtom(new QueryAtom(6),true);
    q1->addAtom(new QueryAtom(8),true);
    q1->addBond(0,1,Bond::UNSPECIFIED);

    // here's the main query
    q2 = new RWMol();
    QueryAtom *qA = new QueryAtom(6);
    RecursiveStructureQuery *rsq = new RecursiveStructureQuery(new RWMol(*q1),1);
    qA->expandQuery(rsq,Queries::COMPOSITE_AND);
    //std::cout << "post expand: " << qA->getQuery() << std::endl;
    q2->addAtom(qA,true,true);
    //std::cout << "mol: " << q2->getAtomWithIdx(0)->getQuery() << std::endl;
    q2->addAtom(new QueryAtom(8),true,true);
    q2->addBond(0,1,Bond::UNSPECIFIED);

    qA = new QueryAtom(6);
    rsq = new RecursiveStructureQuery(new RWMol(*q1),1);
    qA->expandQuery(rsq,Queries::COMPOSITE_AND);
    q2->addAtom(qA,true,true);
    q2->addBond(1,2,Bond::UNSPECIFIED);

    bool found = SubstructMatch(*m,*q2,matchV);
    CHECK_INVARIANT(found,"");
    CHECK_INVARIANT(matchV.size()==3,"");
    n = SubstructMatch(*m,*q2,matches,true);
    TEST_ASSERT(n==1);
    TEST_ASSERT(matches.size()==1);
    TEST_ASSERT(matches[0].size()==3);

    delete q1;
    delete q2;
  }
  delete m;
  std::cout << "Done\n" << std::endl;
}

void test4() {
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "FeatureDef functionality testing." << std::endl;

  ROMol *testMol;
  MatchVectType mv;
  std::string inText;
  int res;
  MolChemicalFeatureDef::CollectionType featureDefs;
  MolChemicalFeatureDef::CollectionType::const_iterator featDefIt;
  MolChemicalFeatureDef::CollectionType::value_type featDef;

  inText =
      "DefineFeature HDonor1 [N,O;!H0]\n"
      "  Family HBondDonor\n"
      "  Weights 1.0\n"
      "EndFeature\n"
      "DefineFeature HAcceptor1 [N,O]\n"
      "  Family HAcceptor\n"
      "  Weights 1.0\n"
      "EndFeature\n";
  res = parseFeatureData(inText, featureDefs);
  TEST_ASSERT(!res);
  TEST_ASSERT(featureDefs.size() == 2);
  featDefIt = featureDefs.begin();
  featDef = *featDefIt;
  TEST_ASSERT(featDef->getSmarts() == "[N,O;!H0]");
  TEST_ASSERT(featDef->getPattern());
  TEST_ASSERT(featDef->getPattern()->getNumAtoms() == 1);
  TEST_ASSERT(featDef->getNumWeights() == 1);

  featDef = *(++featureDefs.begin());
  TEST_ASSERT(featDef->getSmarts() == "[N,O]");
  TEST_ASSERT(featDef->getPattern());
  TEST_ASSERT(featDef->getPattern()->getNumAtoms() == 1);
  TEST_ASSERT(featDef->getNumWeights() == 1);

  testMol = SmilesToMol("COCN");
  TEST_ASSERT(testMol);
  featDef = *featureDefs.begin();
  TEST_ASSERT(SubstructMatch(*testMol, *featDef->getPattern(), mv));
  TEST_ASSERT(mv.size() == 1);
  TEST_ASSERT(mv[0].first == 0);
  TEST_ASSERT(mv[0].second == 3);
  featDef = *(++featureDefs.begin());
  TEST_ASSERT(SubstructMatch(*testMol, *featDef->getPattern(), mv));
  TEST_ASSERT(mv.size() == 1);
  TEST_ASSERT(mv[0].first == 0);
  TEST_ASSERT(mv[0].second == 1 || mv[0].second == 3);
  delete testMol;

  BOOST_LOG(rdErrorLog) << "  done" << std::endl;
}

void test5QueryRoot() {
  std::cout << " ----------------- Test 5 QueryRoot" << std::endl;
  MatchVectType matchV;
  std::vector<MatchVectType> matches;
  int n;
  bool updateLabel = true;
  bool takeOwnership = true;

  RWMol *m, *q1, *q2;
  auto *a6 = new Atom(6);
  auto *a8 = new Atom(8);
  // CC(OC)C
  m = new RWMol();
  m->addAtom(a6);
  m->addAtom(a6);
  m->addAtom(a8);
  m->addAtom(a6);
  m->addAtom(a6);
  m->addBond(0, 1, Bond::SINGLE);
  m->addBond(1, 2, Bond::SINGLE);
  m->addBond(1, 4, Bond::SINGLE);
  m->addBond(2, 3, Bond::SINGLE);

  // this will be the recursive query
  q1 = new RWMol();
  q1->addAtom(new QueryAtom(8), updateLabel, takeOwnership);
  q1->addAtom(new QueryAtom(6), updateLabel, takeOwnership);
  q1->addBond(0, 1, Bond::UNSPECIFIED);
  q1->setProp(common_properties::_queryRootAtom, 1);

  // here's the main query
  q2 = new RWMol();
  auto *qA = new QueryAtom();
  auto *rsq = new RecursiveStructureQuery(q1);
  qA->setQuery(rsq);
  q2->addAtom(qA, true, true);
  q2->addAtom(new QueryAtom(6), true, true);
  q2->addBond(0, 1, Bond::UNSPECIFIED);

  bool found = SubstructMatch(*m, *q2, matchV);
  CHECK_INVARIANT(found, "");
  CHECK_INVARIANT(matchV.size() == 2, "");
  n = SubstructMatch(*m, *q2, matches, true);
  CHECK_INVARIANT(n == 2, "");
  CHECK_INVARIANT(matches[0].size() == 2, "");
  delete m;
  delete a6;
  delete a8;
  delete q2;

  std::cout << "Done\n" << std::endl;
}

void testCrippenO3AConstraints() {
  ROMol *m = SmilesToMol("n1ccc(cc1)-c1ccccc1");
  TEST_ASSERT(m);
  ROMol *m1 = MolOps::addHs(*m);
  delete m;
  TEST_ASSERT(m1);
  DGeomHelpers::EmbedMolecule(*m1);
  MMFF::sanitizeMMFFMol((RWMol &)(*m1));
  MMFF::MMFFMolProperties mp(*m1);
  TEST_ASSERT(mp.isValid());
  ForceFields::ForceField *field = MMFF::constructForceField(*m1, &mp);
  field->initialize();
  field->minimize();
  RWMol *patt = SmartsToMol("nccc-cccc");
  MatchVectType matchVect;
  TEST_ASSERT(SubstructMatch(*m1, (ROMol &)*patt, matchVect));
  unsigned int nIdx = matchVect[0].second;
  unsigned int cIdx = matchVect[matchVect.size() - 1].second;
  MolTransforms::setDihedralDeg(m1->getConformer(), matchVect[2].second,
    matchVect[3].second, matchVect[4].second, matchVect[5].second, 0.0);
  ROMol m2(*m1);
  MolAlign::randomTransform(m2);
  ROMol m3(m2);
  unsigned int prbNAtoms = m2.getNumAtoms();
  std::vector<double> prbLogpContribs(prbNAtoms);
  std::vector<double> prbMRContribs(prbNAtoms);
  std::vector<unsigned int> prbAtomTypes(prbNAtoms);
  std::vector<std::string> prbAtomTypeLabels(prbNAtoms);
  Descriptors::getCrippenAtomContribs(m2, prbLogpContribs,
    prbMRContribs, true, &prbAtomTypes, &prbAtomTypeLabels);
  MolAlign::O3A *o3a = new MolAlign::O3A(m2, *m1,
    &prbLogpContribs, &prbLogpContribs, MolAlign::O3A::CRIPPEN);
  TEST_ASSERT(o3a);
  o3a->align();
  delete o3a;
  double d = (m2.getConformer().getAtomPos(cIdx)
    - m1->getConformer().getAtomPos(cIdx)).length();
  TEST_ASSERT(feq(d, 0.0, 1));
  MatchVectType constraintMap;
  constraintMap.push_back(std::make_pair(cIdx, nIdx));
  o3a = new MolAlign::O3A(m3, *m1, &prbLogpContribs, &prbLogpContribs,
    MolAlign::O3A::CRIPPEN, -1, -1, false, 50, 0, &constraintMap);
  TEST_ASSERT(o3a);
  o3a->align();
  delete o3a;
  d = (m3.getConformer().getAtomPos(cIdx)
    - m1->getConformer().getAtomPos(cIdx)).length();
  TEST_ASSERT(feq(d, 7.0, 1.0));
  delete m1;
}