C++ Conformer类代码示例

void testIssue219(){
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "Testing Issue219: Pickle Failure with Conformations." << std::endl;

  std::string smi="CC";
  ROMol *m1 = SmilesToMol(smi);
  TEST_ASSERT(m1);
  std::string pickle;
  ROMol *m2;
  MolPickler::pickleMol(*m1,pickle);
  m2 = new ROMol();
  MolPickler::molFromPickle(pickle,*m2);
  TEST_ASSERT(m1->getNumAtoms()==m2->getNumAtoms());

  Conformer *conf = new Conformer(2);
  conf->setId(23);
  m1->addConformer(conf);
  MolPickler::pickleMol(*m1,pickle);
  delete m2;
  m2 = new ROMol();
  MolPickler::molFromPickle(pickle,*m2);
  TEST_ASSERT(m1->getNumAtoms()==m2->getNumAtoms());
  TEST_ASSERT(m2->getConformer().getId()==23);
  
  BOOST_LOG(rdErrorLog) << "\tdone" << std::endl;
  
}

static void ConnectTheDots_Medium(RWMol *mol)
{
  int count = mol->getNumAtoms();
  std::vector<ProximityEntry> tmp(count);
  PeriodicTable *table = PeriodicTable::getTable();
  Conformer *conf = &mol->getConformer();
  for (int i=0; i<count; i++) {
    Atom *atom = mol->getAtomWithIdx(i);
    unsigned int elem = atom->getAtomicNum();
    RDGeom::Point3D p = conf->getAtomPos(i);
    ProximityEntry *tmpi = &tmp[i];
    tmpi->x = (float)p.x;
    tmpi->y = (float)p.y;
    tmpi->z = (float)p.z;
    tmpi->r = (float)table->getRcovalent(elem);
    tmpi->atm = i;
  }

  std::stable_sort(tmp.begin(),tmp.end());

  for (int j=0; j<count; j++) {
    ProximityEntry *tmpj = &tmp[j];
    double limit = tmpj->x - MAXDIST;
    for (int k=j-1; k>=0; k--) {
      ProximityEntry *tmpk = &tmp[k];
      if (tmpk->x < limit)
        break;
      if (IsBonded(tmpj,tmpk) &&
          !mol->getBondBetweenAtoms(tmpj->atm,tmpk->atm))
        mol->addBond(tmpj->atm,tmpk->atm,Bond::SINGLE);
    }
  }
}

void computeConfBox(const Conformer &conf, RDGeom::Point3D &leftBottom,
                    RDGeom::Point3D &rightTop, const RDGeom::Transform3D *trans,
                    double padding) {
  double xmin, xmax, ymin, ymax, zmin, zmax;
  xmin = ymin = zmin = 1.e8;
  xmax = ymax = zmax = -1.e8;
  unsigned int i, nAtms = conf.getNumAtoms();
  for (i = 0; i < nAtms; ++i) {
    RDGeom::Point3D loc = conf.getAtomPos(i);
    if (trans) {
      trans->TransformPoint(loc);
    }
    xmax = std::max(xmax, loc.x);
    xmin = std::min(xmin, loc.x);
    ymax = std::max(ymax, loc.y);
    ymin = std::min(ymin, loc.y);
    zmax = std::max(zmax, loc.z);
    zmin = std::min(zmin, loc.z);
  }
  RDGeom::Point3D padPt(padding, padding, padding);
  leftBottom.x = xmin;
  leftBottom.y = ymin;
  leftBottom.z = zmin;
  rightTop.x = xmax;
  rightTop.y = ymax;
  rightTop.z = zmax;
  leftBottom -= padPt;
  rightTop += padPt;
}

void test4Coulomb () {
  std::string path = getenv("RDBASE");
  path += "/Code/GraphMol/MIF/test_data/";

  RWMol mol = *MolFileToMol(path + "HCl.mol", true, false);

  computeGasteigerCharges(mol);

  std::vector<double> charges;
  std::vector<Point3D> pos;
  Conformer conf = mol.getConformer(0);
  for (int i = 0; i < mol.getNumAtoms(); ++i) {
    charges.push_back(mol.getAtomWithIdx (i)->getProp<double> ("_GasteigerCharge"));
    pos.push_back(conf.getAtomPos(i));
  }

  UniformRealValueGrid3D grd = *constructGrid(mol);
  UniformRealValueGrid3D grd2 = *constructGrid(mol);

  Coulomb coul(mol);

  calculateDescriptors<Coulomb>(grd, coul);
  calculateDescriptors<Coulomb>(grd2, Coulomb (charges, pos));

  CHECK_INVARIANT(grd.compareGrids(grd2),
                  "Coulomb: Different constructors do not yield the same descriptor.");
  CHECK_INVARIANT(feq (coul(0.0,0.0,0.0, 1000), 0.0),
                  "Coulomb: Potential between atoms wrong.(should be 0)");
  CHECK_INVARIANT(coul(2.0,0.0,0.0, 1000) < 0,
                  "Coulomb: Potential between positive charges not positive.");
  CHECK_INVARIANT(coul(-2.0,0.0,0.0, 1000) > 0,
                  "Coulomb: Potential between positive and negative charges not negative.");
  CHECK_INVARIANT(feq(coul(0.0,0.0,0.0, 0.1), 0.0),
                  "Coulomb: Small threshold dist does not give 0.");
                  
  calculateDescriptors<Coulomb>(grd, Coulomb(mol, 0, 1.0, true));
  for (unsigned int i = 0; i < grd.getSize(); i++) {
    CHECK_INVARIANT(grd.getVal (i) <= 0.0, "Coulomb: Absolute value field not negative");
  }

Coulomb coul1(mol, 0, -1.0, false, "_GasteigerCharge", 0.0, 0.01);
  CHECK_INVARIANT(coul1(-2.0, 0.0, 0.0, 1000) < 0, "Coulomb: Potential between negative charges not positive.");
  CHECK_INVARIANT(coul1(2.0, 0.0, 0.0, 1000) > 0, "Coulomb: Potential between positive and negative charges not negative.");

  Coulomb coul2 = Coulomb(mol, 0, -.5, false, "_GasteigerCharge", 0.0, 0.01);
  CHECK_INVARIANT(coul1(-2.0, 0.0, 0.0, 1000) < coul2 (-2.0, 0.0, 0.0, 1000),
                  "Coulomb: Higher probecharge does not result in stronger forces.");

  Coulomb coul3(mol, 0, 1.0, false, "_GasteigerCharge", 0.01, 1.0);
  CHECK_INVARIANT(coul3(0.0, 0.0, 0.0, 1000) > coul3(0.1, 0.0, 0.0, 1000),
                  "Coulomb: Softcore interaction wrong.");
  CHECK_INVARIANT(coul3(0.66, 0.0, 0.0, 1000) > coul3(0.68, 0.0, 0.0, 1000),
                  "Coulomb: Softcore interaction wrong.");
  CHECK_INVARIANT(coul3(0.70, 0.0, 0.0, 1000) > coul3(0.68, 0.0, 0.0, 1000),
                  "Coulomb: Softcore interaction wrong.");
  CHECK_INVARIANT(feq(coul3(0.0,0.0,0.0, 0.1), 0.0),
                  "Coulomb: Small threshold dist does not give 0.");

}

void setDihedralRad(Conformer &conf, unsigned int iAtomId, unsigned int jAtomId,
                    unsigned int kAtomId, unsigned int lAtomId, double value) {
  RDGeom::POINT3D_VECT &pos = conf.getPositions();
  URANGE_CHECK(iAtomId, pos.size());
  URANGE_CHECK(jAtomId, pos.size());
  URANGE_CHECK(kAtomId, pos.size());
  URANGE_CHECK(lAtomId, pos.size());
  ROMol &mol = conf.getOwningMol();
  Bond *bondIJ = mol.getBondBetweenAtoms(iAtomId, jAtomId);
  if (!bondIJ) throw ValueErrorException("atoms i and j must be bonded");
  Bond *bondJK = mol.getBondBetweenAtoms(jAtomId, kAtomId);
  if (!bondJK) throw ValueErrorException("atoms j and k must be bonded");
  Bond *bondKL = mol.getBondBetweenAtoms(kAtomId, lAtomId);
  if (!bondKL) throw ValueErrorException("atoms k and l must be bonded");

  if (queryIsBondInRing(bondJK))
    throw ValueErrorException("bond (j,k) must not belong to a ring");
  RDGeom::Point3D rIJ = pos[jAtomId] - pos[iAtomId];
  double rIJSqLength = rIJ.lengthSq();
  if (rIJSqLength <= 1.e-16)
    throw ValueErrorException("atoms i and j have identical 3D coordinates");
  RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
  double rJKSqLength = rJK.lengthSq();
  if (rJKSqLength <= 1.e-16)
    throw ValueErrorException("atoms j and k have identical 3D coordinates");
  RDGeom::Point3D rKL = pos[lAtomId] - pos[kAtomId];
  double rKLSqLength = rKL.lengthSq();
  if (rKLSqLength <= 1.e-16)
    throw ValueErrorException("atoms k and l have identical 3D coordinates");

  RDGeom::Point3D nIJK = rIJ.crossProduct(rJK);
  double nIJKSqLength = nIJK.lengthSq();
  RDGeom::Point3D nJKL = rJK.crossProduct(rKL);
  double nJKLSqLength = nJKL.lengthSq();
  RDGeom::Point3D m = nIJK.crossProduct(rJK);
  // we only need to rotate by delta with respect to the current dihedral value
  value -= -atan2(m.dotProduct(nJKL) / sqrt(nJKLSqLength * m.lengthSq()),
                  nIJK.dotProduct(nJKL) / sqrt(nIJKSqLength * nJKLSqLength));
  // our rotation axis is the (j,k) bond
  RDGeom::Point3D &rotAxisBegin = pos[jAtomId];
  RDGeom::Point3D &rotAxisEnd = pos[kAtomId];
  RDGeom::Point3D rotAxis = rotAxisEnd - rotAxisBegin;
  rotAxis.normalize();
  // get all atoms bonded to k and loop through them
  std::list<unsigned int> alist;
  _toBeMovedIdxList(mol, jAtomId, kAtomId, alist);
  for (std::list<unsigned int>::iterator it = alist.begin(); it != alist.end();
       ++it) {
    // translate atom so that it coincides with the origin of rotation
    pos[*it] -= rotAxisBegin;
    // rotate around our rotation axis
    RDGeom::Transform3D rotByAngle;
    rotByAngle.SetRotation(value, rotAxis);
    rotByAngle.TransformPoint(pos[*it]);
    // translate atom back
    pos[*it] += rotAxisBegin;
  }
}

void test6() {
  BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
  BOOST_LOG(rdErrorLog) << "Feature Location testing." << std::endl;

  ROMol *testMol;
  Conformer *conf;
  std::string inText;
  FeatSPtrList featSPtrs;
  boost::shared_ptr<MolChemicalFeature> featSPtr;

  MolChemicalFeatureFactory *factory;

  MolChemicalFeatureDef::CollectionType::value_type featDef;

  inText =
      "DefineFeature HDonor1 [N,O;!H0]\n"
      "  Family HBondDonor\n"
      "  Weights 1.0\n"
      "EndFeature\n"
      "DefineFeature Carboxyl1 C(=O)[O;H1,-]\n"
      "  Family ZnBinder\n"
      "  Weights 1.0,1.0,1.0\n"
      "EndFeature\n";

  factory = buildFeatureFactory(inText);
  TEST_ASSERT(factory);
  TEST_ASSERT(factory->getNumFeatureDefs() == 2);

  testMol = SmilesToMol("C(=O)O");
  TEST_ASSERT(testMol);
  conf = new Conformer(3);
  testMol->addConformer(conf);
  conf->setAtomPos(0, RDGeom::Point3D(0, 0, 0.0));
  conf->setAtomPos(1, RDGeom::Point3D(1.2, 0, 0.0));
  conf->setAtomPos(2, RDGeom::Point3D(0, 1.5, 0.0));

  featSPtrs = factory->getFeaturesForMol(*testMol);
  TEST_ASSERT(featSPtrs.size() == 2);
  featSPtr = *featSPtrs.begin();
  TEST_ASSERT(featSPtr->getFamily() == "HBondDonor");
  TEST_ASSERT(featSPtr->getType() == "HDonor1");
  TEST_ASSERT(feq(featSPtr->getPos().x, 0.0));
  TEST_ASSERT(feq(featSPtr->getPos().y, 1.5));
  TEST_ASSERT(feq(featSPtr->getPos().z, 0.0));

  featSPtr = *(++featSPtrs.begin());
  TEST_ASSERT(featSPtr->getFamily() == "ZnBinder");
  TEST_ASSERT(featSPtr->getType() == "Carboxyl1");
  TEST_ASSERT(feq(featSPtr->getPos().x, 0.4));
  TEST_ASSERT(feq(featSPtr->getPos().y, 0.5));
  TEST_ASSERT(feq(featSPtr->getPos().z, 0.0));
  delete testMol;

  delete factory;

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

double getDihedralRad(const Conformer &conf, unsigned int iAtomId,
                      unsigned int jAtomId, unsigned int kAtomId,
                      unsigned int lAtomId) {
  const RDGeom::POINT3D_VECT &pos = conf.getPositions();
  URANGE_CHECK(iAtomId, pos.size());
  URANGE_CHECK(jAtomId, pos.size());
  URANGE_CHECK(kAtomId, pos.size());
  URANGE_CHECK(lAtomId, pos.size());
  RDGeom::Point3D rIJ = pos[jAtomId] - pos[iAtomId];
  double rIJSqLength = rIJ.lengthSq();
  if (rIJSqLength <= 1.e-16)
    throw ValueErrorException("atoms i and j have identical 3D coordinates");
  RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
  double rJKSqLength = rJK.lengthSq();
  if (rJKSqLength <= 1.e-16)
    throw ValueErrorException("atoms j and k have identical 3D coordinates");
  RDGeom::Point3D rKL = pos[lAtomId] - pos[kAtomId];
  double rKLSqLength = rKL.lengthSq();
  if (rKLSqLength <= 1.e-16)
    throw ValueErrorException("atoms k and l have identical 3D coordinates");

  RDGeom::Point3D nIJK = rIJ.crossProduct(rJK);
  double nIJKSqLength = nIJK.lengthSq();
  RDGeom::Point3D nJKL = rJK.crossProduct(rKL);
  double nJKLSqLength = nJKL.lengthSq();
  RDGeom::Point3D m = nIJK.crossProduct(rJK);
  // we want a signed dihedral, that's why we use atan2 instead of acos
  return -atan2(m.dotProduct(nJKL) / sqrt(nJKLSqLength * m.lengthSq()),
                nIJK.dotProduct(nJKL) / sqrt(nIJKSqLength * nJKLSqLength));
}

void transformConformer(Conformer &conf, const RDGeom::Transform3D &trans) {
  RDGeom::POINT3D_VECT &positions = conf.getPositions();
  RDGeom::POINT3D_VECT_I pi;
  for (pi = positions.begin(); pi != positions.end(); ++pi) {
    trans.TransformPoint(*pi);
  }
}

RDGeom::Point3D computeCentroid(const Conformer &conf, bool ignoreHs) {
  RDGeom::Point3D res(0.0, 0.0, 0.0);
  const ROMol &mol = conf.getOwningMol();
  ROMol::ConstAtomIterator cai;
  unsigned int nAtms = 0;

  for (cai = mol.beginAtoms(); cai != mol.endAtoms(); cai++) {
    if (((*cai)->getAtomicNum() == 1) && (ignoreHs)) {
      continue;
    }
    res += conf.getAtomPos((*cai)->getIdx());
    nAtms++;
  }
  res /= nAtms;
  return res;
}

static void ConnectTheDots_Large(RWMol *mol)
{
  int HashTable[HASHSIZE];
  memset(HashTable,-1,sizeof(HashTable));

  unsigned int count = mol->getNumAtoms();
  ProximityEntry *tmp = (ProximityEntry*)malloc(count*sizeof(ProximityEntry));
  PeriodicTable *table = PeriodicTable::getTable();
  Conformer *conf = &mol->getConformer();

  for (unsigned int i=0; i<count; i++) {
    Atom *atom = mol->getAtomWithIdx(i);
    unsigned int elem = atom->getAtomicNum();
    RDGeom::Point3D p = conf->getAtomPos(i);
    ProximityEntry *tmpi = tmp+i;
    tmpi->x = (float)p.x;
    tmpi->y = (float)p.y;
    tmpi->z = (float)p.z;
    tmpi->r = (float)table->getRcovalent(elem);
    tmpi->atm = i;

    int hash = HASHX*(int)(p.x/MAXDIST) +
               HASHY*(int)(p.y/MAXDIST) + 
               HASHZ*(int)(p.z/MAXDIST);

    for (int dx = -HASHX; dx <= HASHX; dx += HASHX)
      for (int dy = -HASHY; dy <= HASHY; dy += HASHY)
        for (int dz = -HASHZ; dz <= HASHZ; dz += HASHZ) {
          int probe = hash + dx + dy + dz;
          int list = HashTable[probe & HASHMASK];
          while (list != -1) {
            ProximityEntry *tmpj = &tmp[list];
            if (tmpj->hash == probe &&
                IsBonded(tmpi,tmpj) &&
                !mol->getBondBetweenAtoms(tmpi->atm,tmpj->atm))
              mol->addBond(tmpi->atm,tmpj->atm,Bond::SINGLE);
            list = tmpj->next;
          }
        }

    int list = hash & HASHMASK;
    tmpi->next =HashTable[list];
    HashTable[list] = i;
    tmpi->hash = hash;
  }
  free(tmp);
}

void setAngleRad(Conformer &conf, unsigned int iAtomId, unsigned int jAtomId,
                 unsigned int kAtomId, double value) {
  RDGeom::POINT3D_VECT &pos = conf.getPositions();
  URANGE_CHECK(iAtomId, pos.size());
  URANGE_CHECK(jAtomId, pos.size());
  URANGE_CHECK(kAtomId, pos.size());
  ROMol &mol = conf.getOwningMol();
  Bond *bondJI = mol.getBondBetweenAtoms(jAtomId, iAtomId);
  if (!bondJI) throw ValueErrorException("atoms i and j must be bonded");
  Bond *bondJK = mol.getBondBetweenAtoms(jAtomId, kAtomId);
  if (!bondJK) throw ValueErrorException("atoms j and k must be bonded");
  if (queryIsBondInRing(bondJI) && queryIsBondInRing(bondJK))
    throw ValueErrorException(
        "bonds (i,j) and (j,k) must not both belong to a ring");

  RDGeom::Point3D rJI = pos[iAtomId] - pos[jAtomId];
  double rJISqLength = rJI.lengthSq();
  if (rJISqLength <= 1.e-16)
    throw ValueErrorException("atoms i and j have identical 3D coordinates");
  RDGeom::Point3D rJK = pos[kAtomId] - pos[jAtomId];
  double rJKSqLength = rJK.lengthSq();
  if (rJKSqLength <= 1.e-16)
    throw ValueErrorException("atoms j and k have identical 3D coordinates");

  // we only need to rotate by delta with respect to the current angle value
  value -= rJI.angleTo(rJK);
  RDGeom::Point3D &rotAxisBegin = pos[jAtomId];
  // our rotation axis is the normal to the plane of atoms i, j, k
  RDGeom::Point3D rotAxisEnd = rJI.crossProduct(rJK) + pos[jAtomId];
  RDGeom::Point3D rotAxis = rotAxisEnd - rotAxisBegin;
  rotAxis.normalize();
  // get all atoms bonded to j and loop through them
  std::list<unsigned int> alist;
  _toBeMovedIdxList(mol, jAtomId, kAtomId, alist);
  for (std::list<unsigned int>::iterator it = alist.begin(); it != alist.end();
       ++it) {
    // translate atom so that it coincides with the origin of rotation
    pos[*it] -= rotAxisBegin;
    // rotate around our rotation axis
    RDGeom::Transform3D rotByAngle;
    rotByAngle.SetRotation(value, rotAxis);
    rotByAngle.TransformPoint(pos[*it]);
    // translate atom back
    pos[*it] += rotAxisBegin;
  }
}

double getBondLength(const Conformer &conf, unsigned int iAtomId,
                     unsigned int jAtomId) {
  const RDGeom::POINT3D_VECT &pos = conf.getPositions();
  URANGE_CHECK(iAtomId, pos.size());
  URANGE_CHECK(jAtomId, pos.size());

  return (pos[iAtomId] - pos[jAtomId]).length();
}

void OptimizeLigandGeomety(Conformer& ligand, double chi_increment, std::vector<double> chi_stdev, std::string lp_atom,
                           double unit_cell_length, double angle_tolerance, std::string output_dir) {
  if (ligand.parent()->conformer_states() == NULL)
    Log->error(FLERR, "There are no conformer states for the ligand!");
  const std::vector<Linkage<std::string> >& chi_names = ligand.parent()->conformer_states()->degrees_of_freedom();

  std::vector<double> chi_values = GetChiValues(ligand);
  //adjust chi to values between [0, 360]
  for (size_t i = 0; i < chi_values.size(); ++i) {
    chi_values[i] = (chi_values[i] < 0 ? chi_values[i] += 360 : chi_values[i]);
  }

  //loop through the possible chi angles, filling the ligand geometry map
  std::map<std::vector<double>, std::vector<double> > ligand_geometry;
  for(double chi1 = floor(chi_values[0] - chi_stdev[0]); chi1 < ceil(chi_values[0] + chi_stdev[0]); chi1 += chi_increment) {
    for(double chi2 = floor(chi_values[1] - chi_stdev[1]); chi2 < ceil(chi_values[1] + chi_stdev[1]); chi2 += chi_increment) {
      ligand.AdjustLinkage(chi_names[0], chi1);
      ligand.AdjustLinkage(chi_names[1], chi2);

      std::vector<double> chis;
      chis.push_back(chi1);
      chis.push_back(chi2);
      Vector3 v = CalculateLonePairDirection(ligand, lp_atom);

      //Calculate the dot product of v and z axis unit vector
      //cos(q) = v*w / (|v|*|w|)
      double angle_z = RADIANS_TO_DEGREES * acos(v.DotProduct(Vector3(0, 0, 1)) / v.Length());
      //Calculate the dot product of v and the vector from the N to the metal (where the C4 symmetry axis is)
      Atom* N = ligand.Find(lp_atom);
      Vector3 N_metal(unit_cell_length - N->x(), - N->y(), 0);
      double angle_metal = RADIANS_TO_DEGREES * acos(v.DotProduct(N_metal) / (v.Length() * N_metal.Length()));

      if ((angle_z < (90.0 + angle_tolerance) && angle_z > (90.0 - angle_tolerance)) && (angle_metal < angle_tolerance)) {
        ligand_geometry[chis].push_back(angle_z);
        ligand_geometry[chis].push_back(angle_metal);

        //OutputStructure(ligand.parent()->parent()->parent()->parent()->parent()->parent(), chis, lp_atom, output_dir);
        Log->print("Optimized Cu binding PDB with chi1 " + Log->to_str(chi1) + ", chi2 " + Log->to_str(chi2));

        OutputLatticeStructure(ligand.parent()->parent()->parent()->parent()->parent()->parent(), chis, unit_cell_length, lp_atom, output_dir);
      }
    }
  }

  OutputStat(ligand_geometry, lp_atom, output_dir);
}

MOL_SPTR_VECT generateOneProductSet(const ChemicalReaction& rxn,
                                    const MOL_SPTR_VECT& reactants,
                                    const std::vector<MatchVectType> &reactantsMatch)
{
    PRECONDITION(reactants.size() == reactantsMatch.size(),"vector size mismatch");

    // if any of the reactants have a conformer, we'll go ahead and
    // generate conformers for the products:
    bool doConfs=false;
    BOOST_FOREACH(ROMOL_SPTR reactant,reactants) {
        if(reactant->getNumConformers()) {
            doConfs=true;
            break;
        }
    }

    MOL_SPTR_VECT res;
    res.resize(rxn.getNumProductTemplates());
    unsigned int prodId=0;
    for(MOL_SPTR_VECT::const_iterator pTemplIt = rxn.beginProductTemplates();
            pTemplIt != rxn.endProductTemplates(); ++pTemplIt) {
        // copy product template and its properties to a new product RWMol
        RWMOL_SPTR product = initProduct(*pTemplIt);
        Conformer *conf = 0;
        if(doConfs) {
            conf = new Conformer();
            conf->set3D(false);
        }

        unsigned int reactantId = 0;
        for(MOL_SPTR_VECT::const_iterator iter = rxn.beginReactantTemplates();
                iter != rxn.endReactantTemplates(); ++iter, reactantId++) {
            addReactantAtomsAndBonds(rxn, product, reactants.at(reactantId),
                                     reactantsMatch.at(reactantId),
                                     *iter, conf);
        }

        if(doConfs) {
            product->addConformer(conf,true);
        }
        res[prodId] = product;
        ++prodId;
    }
    return res;
}

void _fillAtomPositions(RDGeom::Point3DConstPtrVect &pts, const Conformer &conf,
                        const std::vector<unsigned int> *atomIds = 0) {
  unsigned int na = conf.getNumAtoms();
  pts.clear();
  if (atomIds == 0) {
    unsigned int ai;
    pts.reserve(na);
    for (ai = 0; ai < na; ++ai) {
      pts.push_back(&conf.getAtomPos(ai));
    }
  } else {
    pts.reserve(atomIds->size());
    std::vector<unsigned int>::const_iterator cai;
    for (cai = atomIds->begin(); cai != atomIds->end(); cai++) {
      pts.push_back(&conf.getAtomPos(*cai));
    }
  }
}