C++ Conformer::getAtomPos方法代码示例

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;
}

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 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 _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));
    }
  }
}

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 EncodeShape(const Conformer &conf, RDGeom::UniformGrid3D &grid, 
                  const RDGeom::Transform3D *trans, double vdwScale, 
                  double stepSize, int maxLayers, bool ignoreHs) {
   const ROMol &mol = conf.getOwningMol();
   ROMol::ConstAtomIterator ai;
   double rad;
   unsigned int aid, anum;
   for (ai = mol.beginAtoms(); ai != mol.endAtoms(); ai++) {
     anum = (*ai)->getAtomicNum();
     if ((anum == 1) && (ignoreHs)) { //ignore hydrigens
       continue;
     }
     aid = (*ai)->getIdx();
     RDGeom::Point3D loc = conf.getAtomPos(aid);
     rad = PeriodicTable::getTable()->getRvdw(anum);
     if (trans) {
       trans->TransformPoint(loc);
     }
     grid.setSphereOccupancy(loc, vdwScale*rad, stepSize, maxLayers);
   }
 }

static void ConnectTheDots_Small(RWMol *mol)
{
  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);
    for (unsigned int j=0; j<i; j++) {
      ProximityEntry *tmpj = tmp+j;
      if (IsBonded(tmpi,tmpj) && !mol->getBondBetweenAtoms(i,j))
        mol->addBond(i,j,Bond::SINGLE);
    }
  }
  free(tmp);
}

void test1Canonicalization() {
  ROMol *mol = SmilesToMol("C", 0, 1);
  Conformer *conf = new Conformer(1);
  conf->setAtomPos(0, RDGeom::Point3D(4.0, 5.0, 6.0));
  int cid = mol->addConformer(conf, true);
  RDGeom::Point3D pt = computeCentroid(*conf);
  CHECK_INVARIANT(comparePts(pt, RDGeom::Point3D(4.0, 5.0, 6.0)), "");
  
  RDGeom::Transform3D *trans = computeCanonicalTransform(*conf);
  transformConformer(*conf, *trans);
  CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(0.0, 0.0, 0.0)), "");

  conf->setAtomPos(0, RDGeom::Point3D(4.0, 5.0, 6.0));
  canonicalizeConformer(*conf);
  CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(0.0, 0.0, 0.0)), "");

  delete mol;
  //delete conf;
  delete trans;
  // lets try two points now
  mol = SmilesToMol("CC", 0, 1);
  conf = new Conformer(2);
  conf->setAtomPos(0, RDGeom::Point3D(0.0, 0.0, 0.0));
  conf->setAtomPos(1, RDGeom::Point3D(1.5, 0.0, 0.0));
  cid = mol->addConformer(conf, true);
  trans = computeCanonicalTransform(*conf);
  canonicalizeConformer(*conf);
  CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(-0.75, 0.0, 0.0)), "");
  CHECK_INVARIANT(comparePts(conf->getAtomPos(1), RDGeom::Point3D(0.75, 0.0, 0.0)), "");

  conf->setAtomPos(0, RDGeom::Point3D(0.0, 0.0, 0.0));
  conf->setAtomPos(1, RDGeom::Point3D(0.0, 1.5, 0.0));
  trans = computeCanonicalTransform(*conf);
  canonicalizeConformer(*conf);
  
  CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(-0.75, 0.0, 0.0)), "");
  CHECK_INVARIANT(comparePts(conf->getAtomPos(1), RDGeom::Point3D(0.75, 0.0, 0.0)), "");
  delete mol;
  delete trans;
  
  mol = SmilesToMol("CC", 0, 1);
  conf = new Conformer(2);
  conf->setAtomPos(0, RDGeom::Point3D(0.0, 0.0, 0.0));
  conf->setAtomPos(1, RDGeom::Point3D(1.5, 0.0, 0.0));
  cid = mol->addConformer(conf, true);
  trans = computeCanonicalTransform(*conf);
  transformConformer(*conf, *trans);
  canonicalizeConformer(*conf);
  CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(-0.75, 0.0, 0.0)), "");
  CHECK_INVARIANT(comparePts(conf->getAtomPos(1), RDGeom::Point3D(0.75, 0.0, 0.0)), "");
  delete mol;
  delete trans;

  mol = SmilesToMol("C1CC1", 0, 1);
  conf = new Conformer(3);
  conf->setAtomPos(0, RDGeom::Point3D(0.58, -0.66, -0.08));
  conf->setAtomPos(1, RDGeom::Point3D(-0.88, -0.18, -0.04));
  conf->setAtomPos(2, RDGeom::Point3D(.26, 0.82, 0.14));
  cid = mol->addConformer(conf, true);
  //trans = computeCanonicalTransform(*conf);
  //transformConformer(*conf, *trans);
  canonicalizeConformer(*conf);
  CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(-0.6418, 0.6158, 0.0)), "");
  CHECK_INVARIANT(comparePts(conf->getAtomPos(1), RDGeom::Point3D(-0.2029, -0.8602, 0.0)), "");
  CHECK_INVARIANT(comparePts(conf->getAtomPos(2), RDGeom::Point3D(0.8447, 0.2445, 0.0)), "");
  MolToMolFile(*mol, "junk.mol", 0);
  //CHECK_INVARIANT(comparePts(conf->getAtomPos(0), RDGeom::Point3D(-0.75, 0.0, 0.0)), "");
  //CHECK_INVARIANT(comparePts(conf->getAtomPos(1), RDGeom::Point3D(0.75, 0.0, 0.0)), "");
  delete mol;
  
  std::string rdbase = getenv("RDBASE");
  std::string fname1 = rdbase + "/Code/GraphMol/MolTransforms/test_data/1oir.mol";
  mol = MolFileToMol(fname1);
  std::string fname2 = rdbase + "/Code/GraphMol/MolTransforms/test_data/1oir_canon.mol";
  ROMol *mol2 = MolFileToMol(fname2);
  
  Conformer &conf1 = mol->getConformer(0); 
  canonicalizeConformer(conf1);

  Conformer &conf2 = mol2->getConformer();
  unsigned int i, nats = mol->getNumAtoms();
  for (i = 0; i < nats; ++i) {
    CHECK_INVARIANT(comparePts(conf1.getAtomPos(i), conf2.getAtomPos(i)), "");
  }


  delete mol;
  delete mol2;
}