C++ AtomList::add方法代码示例

AtomList AtomList::neighborhood(Atomo *a)
{
    AtomList al = AtomList();
    for (int i = 0; i < _size; i++) {
	Atomo *b = &contents[i];
	if (closeEnough(a, b))
	    al.add(*b);
    }
    return al;
}

AtomList* MolData::findAtomsGrid(const AtomGrid* grid) 
{
  // Scan atoms and check adjacent cells
  AtomList *al = new AtomList(count);
  const double r = getActualRadius();
  const double rsq = r*r;
  int atomi;
  for(atomi=0; atomi<count; atomi++) {
    Vector3D atom = coords[atomi].coord;
//    std::cout << "Checking atom " << atomi << "(" << atom.x << "," << atom.y << "," << atom.z <<std::endl;
    int i0, j0, k0;
    bool found = false;
    
    grid->get_cijk(atom,i0,j0,k0);
    int di;
    for(di=-1; !found && di <= 1; di++) {
      int ci = i0 + di;
      int cna = grid->get_cna();
      int wrapi = 0;
      if (ci<0) {
        ci += cna;
        wrapi = grid->get_na();
      } else if (ci >= cna) {
        ci -= cna;
        wrapi = -grid->get_na();
      }
      int dj;
      for(dj=-1; !found && dj <= 1; dj++) {
        int cj = j0 + dj;
        int cnb = grid->get_cnb();
        int wrapj = 0;
        if (cj<0) {
          cj += cnb;
          wrapj = grid->get_nb();
        } else if (cj >= cnb) {
          cj -= cnb;
          wrapj = grid->get_nb();
        }
        int dk;
        for(dk=-1; !found && dk <= 1; dk++) {
          int ck = k0 + dk;
          int cnc = grid->get_cnc();
          int wrapk = 0;
          if (ck<0) {
            ck += cnc;
            wrapk = grid->get_nc();
          } else if (ck >= cnc) {
            ck -= cnc;
            wrapk = -grid->get_nc();
          }
          const GridList* list = grid->get_cell(ci,cj,ck);
//          std::cout << "Found " << list->getNum() << " atoms in cell " 
//            << ci << "," << cj << "," << ck << std::endl;
          int gli;
          for(gli=0; !found && gli < list->getNum(); gli++) {
            Vector3D gridpt = list->get(gli);
            gridpt.x += wrapi;
            gridpt.y += wrapj;
            gridpt.z += wrapk;
            const Vector3D v = grid->get_xyz(gridpt);
            const Vector3D dv = v - atom;
            if (dv.length2() <= rsq) {
              al->add(coords[atomi].index);
              found = true;
//              std::cout << "Found vac " << gli << std::endl;
            }
          }
        }
      }
    }
  }
  return al;
}

AtomList* MolData::findAtoms(const AtomGrid* grid) 
{
  // Find the bounding box for the selection distance
  Vector3D cornerx[8];
  Vector3D dx, dy, dz;
  
  dx.y = dx.z = dy.x = dy.z = dz.x = dz.y = 0;
  dx.x = dy.y = dz.z = getActualRadius();
  cornerx[0] = origin;
  cornerx[1] = origin + dx;
  cornerx[2] = origin + dy;
  cornerx[3] = origin + dz;
  cornerx[4] = cornerx[1] + dy;
  cornerx[5] = cornerx[1] + dz;
  cornerx[6] = cornerx[2] + dz;
  cornerx[7] = cornerx[4] + dz;  
  
  // Find the max i,j,k for those corners
  Vector3D ijk = grid->get_ijk(cornerx[0]);
  double dimax = ijk.x;
  double djmax = ijk.y;
  double dkmax = ijk.z;

  int i;  
  for(i=1; i < 8; i++) {
    ijk = grid->get_ijk(cornerx[i]);
    if (ijk.x > dimax)
      dimax = ijk.x;
    if (ijk.y > djmax)
      djmax = ijk.y;
    if (ijk.z > dkmax)
      dkmax = ijk.z;
  }
  
  const int imax = (int)ceil(dimax);
  const int jmax = (int)ceil(djmax);
  const int kmax = (int)ceil(dkmax);

  // Find the bounding box for the atom radius
  Vector3D bb = grid->findAtomBox(getActualRadius());
  
  int* ai_s = AtomGrid::scramble_array(2*imax+1);
  int* aj_s = AtomGrid::scramble_array(2*jmax+1);
  int* ak_s = AtomGrid::scramble_array(2*kmax+1);
    
//  std::cout << "Max Box=" << imax << " " << jmax << " " << kmax << std::endl;
  
  // Scan atoms and check adjacent cells
  AtomList *al = new AtomList(count);
  const double r = getActualRadius();
  const double rsq = r*r;
  int atomi;
  for(atomi=0; atomi<count; atomi++) {
    bool found = false;
    Vector3D atom = coords[atomi].coord;
    Vector3D atomijk = grid->get_ijk(atom);
    const int ai = (int)atomijk.x;
    const int aj = (int)atomijk.y;
    const int ak = (int)atomijk.z;

    const int na = grid->get_na();
    const int nb = grid->get_nb();
    const int nc = grid->get_nc();

//    for(int is=-imax; !found && is <= imax; is++) {
//      int i = ai+is;
    int is;
    for(is=0; !found && is < 2*imax+1; is++) {
      int i = ai+(ai_s[is] - imax);
      int ii = i;
      if (i < 0)
        ii += na;
      else if (i >= na)
        ii -= na;
//      for(int js=-jmax; !found && js <= jmax; js++) {
//        int j = aj+js;
      int js;
      for(js=0; !found && js < 2*jmax+1; js++) {
        int j = aj + (aj_s[js] - jmax);
        int jj = j;
        if (j < 0)
          jj += nb;
        else if (j >= nb)
          jj -= nb;
        int ks;
        for(ks=0; !found && ks <= kmax; ks++) {
          int k = ak + ks;
          int kk = k;
          if (k < 0)
            kk += nc;
          else if (kk >= nc)
            kk -= nc;

          // See if its a vacuum
          if (!grid->get(ii,jj,kk)) {
            const Vector3D v = grid->get_xyz(i,j,k);
            const Vector3D dv = v - atom;
            if (dv.length2() <= rsq) {
              al->add(coords[atomi].index);
              found = true;
//.........这里部分代码省略.........

AtomList* MolData::findAtomsGrid(const AtomGrid* grid) 
{
   // Scan atoms and check adjacent cells
   AtomList *al = new AtomList(count);
   const double r = getActualRadius();
   const double rsq = r*r;
   const Vector3D atom_box = grid->findAtomBox(radius);
   int atomi;
   for(atomi=0; atomi<count; atomi++) {
      Vector3D atom = coords[atomi].coord;
      //    std::cout << "Checking atom " << atomi << "(" << atom.x << "," << atom.y << "," << atom.z <<std::endl;
      int i0, j0, k0;
      bool found = false;
      
      grid->get_cijk(atom,i0,j0,k0);
      int di;
      bool cell_outside = false;
      for(di=-1; !found && di <= 1; di++) {
         int ci = i0 + di;
         int cna = grid->get_cna();
         int wrapi = 0;
         if (ci<0) {
           ci += cna;
           wrapi = grid->get_na();
           if (!periodic)
             cell_outside = true;
         } else if (ci >= cna) {
           ci -= cna;
           wrapi = -grid->get_na();
           if (!periodic)
             cell_outside = true;
         }
         int dj;
         for(dj=-1; !found && dj <= 1; dj++) {
            int cj = j0 + dj;
            int cnb = grid->get_cnb();
            int wrapj = 0;
            if (cj<0) {
              cj += cnb;
              wrapj = grid->get_nb();
              if (!periodic)
                cell_outside = true;
            } else if (cj >= cnb) {
              cj -= cnb;
              wrapj = grid->get_nb();
              if (!periodic)
                cell_outside = true;
            }
            int dk;
            for(dk=-1; !found && dk <= 1; dk++) {
               int ck = k0 + dk;
               int cnc = grid->get_cnc();
               int wrapk = 0;
               if (ck<0) {
                 ck += cnc;
                 wrapk = grid->get_nc();
                 if (!periodic)
                   cell_outside = true;
               } else if (ck >= cnc) {
                 ck -= cnc;
                 wrapk = -grid->get_nc();
                 if (!periodic)
                   cell_outside = true;
               }
               // If the cell is outside the box, then there must be a
               // vacuum there
               if (!periodic && cell_outside) {
                 const Vector3D atom_grid = grid->get_ijk(atom);
                 const Vector3D corner0 = atom_grid - atom_box;
                 bool addit = false;
                 if (corner0.x < 0 || corner0.y < 0 || corner0.z < 0) {
                   addit = true;
                 } else {
                   const Vector3D corner1 = atom_grid + atom_box;
                   if (corner1.x >= grid->get_na() 
                       || corner1.y >= grid->get_nb() 
                       || corner1.z >= grid->get_nc() )
                     addit = true;
                 }
                 if (addit) {
                   al->add(coords[atomi].index);
                   found = true;
                 }
                 cell_outside = false;
               }
               const GridList* list = grid->get_cell(ci,cj,ck);
//             std::cout << "Found " << list->getNum() 
//                       << " atoms in cell " 
//                       << ci << "," << cj << "," << ck << std::endl;
               int gli;
               for(gli=0; !found && gli < list->getNum(); gli++) {
//                XXX: this does not handle get() receiving an illegal argument.
                  Vector3D gridpt = list->get(gli);
                  gridpt.x += wrapi;
                  gridpt.y += wrapj;
                  gridpt.z += wrapk;
                  const Vector3D v = grid->get_xyz(gridpt);
                  const Vector3D dv = v - atom;
                  if (dv.length2() <= rsq) {
                     al->add(coords[atomi].index);
//.........这里部分代码省略.........