C++ Molecule::atom方法代码示例

TEST(RingPerceiverTest, ethanol)
{
  Molecule molecule;
  molecule.addAtom(6);
  molecule.addAtom(6);
  molecule.addAtom(8);
  molecule.addBond(molecule.atom(0), molecule.atom(1), 1);
  molecule.addBond(molecule.atom(1), molecule.atom(2), 1);

  RingPerceiver perceiver(&molecule);
  std::vector<std::vector<size_t>> rings = perceiver.rings();
  EXPECT_EQ(rings.size(), static_cast<size_t>(0));
}

TEST(CjsonTest, crystal)
{
  CjsonFormat cjson;
  Molecule molecule;
  bool success = cjson.readFile(std::string(AVOGADRO_DATA) +
                                "/data/rutile.cjson", molecule);
  EXPECT_TRUE(success);
  EXPECT_EQ(cjson.error(), "");
  EXPECT_EQ(molecule.data("name").toString(), "TiO2 rutile");
  EXPECT_EQ(molecule.atomCount(), static_cast<size_t>(6));
  EXPECT_EQ(molecule.bondCount(), static_cast<size_t>(0));

  const UnitCell *unitCell = molecule.unitCell();
  ASSERT_NE(unitCell, (UnitCell*)NULL);
  EXPECT_TRUE(std::fabs((float)unitCell->a() - 2.95812f) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)unitCell->b() - 4.59373f) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)unitCell->c() - 4.59373f) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)unitCell->alpha() - (.5f * PI_F)) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)unitCell->beta()  - (.5f * PI_F)) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)unitCell->gamma() - (.5f * PI_F)) < 1e-5f);

  Atom atom = molecule.atom(5);
  EXPECT_EQ(atom.atomicNumber(), 8);
  EXPECT_TRUE(std::fabs((float)atom.position3d().x() - 1.479060f) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)atom.position3d().y() - 3.699331f) < 1e-5f);
  EXPECT_TRUE(std::fabs((float)atom.position3d().z() - 0.894399f) < 1e-5f);

  std::string cjsonStr;
  cjson.writeString(cjsonStr, molecule);
  Molecule otherMolecule;
  cjson.readString(cjsonStr, otherMolecule);

  const UnitCell *otherUnitCell = otherMolecule.unitCell();
  ASSERT_NE(otherUnitCell, (UnitCell*)NULL);
  EXPECT_FLOAT_EQ((float)otherUnitCell->a(),     (float)unitCell->a());
  EXPECT_FLOAT_EQ((float)otherUnitCell->b(),     (float)unitCell->b());
  EXPECT_FLOAT_EQ((float)otherUnitCell->c(),     (float)unitCell->c());
  EXPECT_FLOAT_EQ((float)otherUnitCell->alpha(), (float)unitCell->alpha());
  EXPECT_FLOAT_EQ((float)otherUnitCell->beta(),  (float)unitCell->beta());
  EXPECT_FLOAT_EQ((float)otherUnitCell->gamma(), (float)unitCell->gamma());

  Atom otherAtom = otherMolecule.atom(5);
  EXPECT_EQ(otherAtom.atomicNumber(), atom.atomicNumber());
  EXPECT_FLOAT_EQ((float)otherAtom.position3d().x(),
                  (float)atom.position3d().x());
  EXPECT_FLOAT_EQ((float)otherAtom.position3d().y(),
                  (float)atom.position3d().y());
  EXPECT_FLOAT_EQ((float)otherAtom.position3d().z(),
                  (float)atom.position3d().z());
}

   /*
   * Allocate arrays.
   */
   void G1MSD::setup() 
   {
 
      // Get number of molecules of this species in the System. 
      nMolecule_ = system().nMolecule(speciesId_);
      
      int nratoms;
      nratoms =  nMolecule_*nAtom_;
      // Allocate arrays of position Vectors and shifts
      truePositions_.allocate(nratoms); 
      oldPositions_.allocate(nratoms); 
      shifts_.allocate(nratoms); 

      // Initialize the AutoCorrArray object
      accumulator_.setParam(nratoms, capacity_);

      // Store initial positions, and set initial shift vectors.
      Vector     r;
      IntVector  zero(0);
      Molecule*  moleculePtr;
      int        iatom = 0;
      for (int i = 0; i < nMolecule_; ++i) {
	 moleculePtr = &system().molecule(speciesId_, i);
	 for (int j = 0 ; j < nAtom_; j++) {
            r = moleculePtr->atom(j).position();
            system().boundary().shift(r);
            oldPositions_[iatom] = r;
            shifts_[iatom] = zero;
	    iatom++;
	 }
      }

   }

   /*
   * Initialize all Angle objects for Molecules of one Species.
   *
   * This functions assigns pointers to Atoms and angle types ids within a
   * contiguous block of Angle objects, and sets a pointer in each Molecule
   * to the first Angle in the associated block.
   */
   void Simulation::initializeSpeciesAngles(int iSpecies)
   {

      if (nAngleType_ <= 0) {
         UTIL_THROW("nAngleType must be positive");
      }

      Species*  speciesPtr = 0;
      Molecule *moleculePtr = 0;
      Angle    *anglePtr = 0;
      Atom     *firstAtomPtr, *atom0Ptr, *atom1Ptr, *atom2Ptr;
      int       iMol, iAngle, atom0Id, atom1Id, atom2Id, type;
      int       capacity, nAngle;

      speciesPtr = &species(iSpecies);
      capacity   = speciesPtr->capacity();
      nAngle     = speciesPtr->nAngle();

      // Initialize pointers before loop
      moleculePtr = &molecules_[firstMoleculeIds_[iSpecies]];
      anglePtr = &angles_[firstAngleIds_[iSpecies]];

      // Loop over molecules in Species
      for (iMol = 0; iMol < capacity; ++iMol) {

         firstAtomPtr = &(moleculePtr->atom(0));
         moleculePtr->setFirstAngle(*anglePtr);
         moleculePtr->setNAngle(nAngle);

         if (nAngle > 0) {

            // Create angles for a molecule
            for (iAngle = 0; iAngle < nAngle; ++iAngle) {

               // Get pointers to atoms spanning the angle and angle type
               atom0Id  = speciesPtr->speciesAngle(iAngle).atomId(0);
               atom1Id  = speciesPtr->speciesAngle(iAngle).atomId(1);
               atom2Id  = speciesPtr->speciesAngle(iAngle).atomId(2);
               type     = speciesPtr->speciesAngle(iAngle).typeId();
               atom0Ptr = firstAtomPtr + atom0Id;
               atom1Ptr = firstAtomPtr + atom1Id;
               atom2Ptr = firstAtomPtr + atom2Id;

               // Set fields of the Angle object
               anglePtr->setAtom(0, *atom0Ptr);
               anglePtr->setAtom(1, *atom1Ptr);
               anglePtr->setAtom(2, *atom2Ptr);
               anglePtr->setTypeId(type);

               ++anglePtr;

            }

         }

         ++moleculePtr;
      }

   }

void BallAndStick::process(const Molecule &molecule,
                           Rendering::GroupNode &node)
{
  // Add a sphere node to contain all of the spheres.
  GeometryNode *geometry = new GeometryNode;
  node.addChild(geometry);
  SphereGeometry *spheres = new SphereGeometry;
  spheres->identifier().molecule = &molecule;
  spheres->identifier().type = Rendering::AtomType;
  geometry->addDrawable(spheres);

  for (Index i = 0; i < molecule.atomCount(); ++i) {
    Core::Atom atom = molecule.atom(i);
    unsigned char atomicNumber = atom.atomicNumber();
    const unsigned char *c = Elements::color(atomicNumber);
    Vector3ub color(c[0], c[1], c[2]);
    spheres->addSphere(atom.position3d().cast<float>(), color,
                       static_cast<float>(Elements::radiusVDW(atomicNumber))
                       * 0.3f);
  }

  float bondRadius = 0.1f;
  CylinderGeometry *cylinders = new CylinderGeometry;
  cylinders->identifier().molecule = &molecule;
  cylinders->identifier().type = Rendering::BondType;
  geometry->addDrawable(cylinders);
  for (Index i = 0; i < molecule.bondCount(); ++i) {
    Core::Bond bond = molecule.bond(i);
    Vector3f pos1 = bond.atom1().position3d().cast<float>();
    Vector3f pos2 = bond.atom2().position3d().cast<float>();
    Vector3ub color1(Elements::color(bond.atom1().atomicNumber()));
    Vector3ub color2(Elements::color(bond.atom2().atomicNumber()));
    Vector3f bondVector = pos2 - pos1;
    float bondLength = bondVector.norm();
    bondVector /= bondLength;
    switch (bond.order()) {
    case 3: {
      Vector3f delta = bondVector.unitOrthogonal() * (2.0f * bondRadius);
      cylinders->addCylinder(pos1 + delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
      cylinders->addCylinder(pos1 - delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
    }
    default:
    case 1:
      cylinders->addCylinder(pos1, bondVector, bondLength, bondRadius,
                             color1, color2, i);
      break;
    case 2: {
      Vector3f delta = bondVector.unitOrthogonal() * bondRadius;
      cylinders->addCylinder(pos1 + delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
      cylinders->addCylinder(pos1 - delta, bondVector, bondLength, bondRadius,
                             color1, color2, i);
    }
    }
  }
}

static PyObject *setbonds(PyObject *self, PyObject *args) {
  int molid;
  PyObject *atomlist, *bondlist; 

  if (!PyArg_ParseTuple(args, (char *)"iO!O!:setbonds", &molid, 
                        &PyTuple_Type, &atomlist, &PyList_Type, &bondlist)) 
    return NULL;  // bad args
 
  Molecule *mol = get_vmdapp()->moleculeList->mol_from_id(molid);
  if (!mol) {
    PyErr_SetString(PyExc_ValueError, "molecule no longer exists");
    return NULL;
  }
  int num_atoms = mol->nAtoms;
  int num_selected = PyTuple_Size(atomlist);
  if (PyList_Size(bondlist) != num_selected) {
    PyErr_SetString(PyExc_ValueError, 
      (char *)"setbonds: atomlist and bondlist must have the same size");
    return NULL;
  }
  mol->force_recalc(DrawMolItem::MOL_REGEN); // many reps ignore bonds
  for (int i=0; i<num_selected; i++) {
    int id = PyInt_AsLong(PyTuple_GET_ITEM(atomlist, i));
    if (PyErr_Occurred()) {
      return NULL;
    }
    if (id < 0 || id >= num_atoms) {
      PyErr_SetString(PyExc_ValueError, (char *)"invalid atom id found");
      return NULL;
    }
    MolAtom *atom = mol->atom(id);
   
    PyObject *atomids = PyList_GET_ITEM(bondlist, i);
    if (!PyList_Check(atomids)) {
      PyErr_SetString(PyExc_TypeError, 
        (char *)"bondlist must contain lists");
      return NULL;
    }
    int numbonds = PyList_Size(atomids);
    int k=0;
    for (int j=0; j<numbonds; j++) {
      int bond = PyInt_AsLong(PyList_GET_ITEM(atomids, j));
      if (PyErr_Occurred())
        return NULL;
      if (bond >= 0 && bond < mol->nAtoms) {
        atom->bondTo[k++] = bond;
      } else {
        char buf[40];
        sprintf(buf, "Invalid atom id in bondlist: %d", bond);
        PyErr_SetString(PyExc_ValueError, buf);
        return NULL;
      }
    }
    atom->bonds = k;
  }
  Py_INCREF(Py_None);
  return Py_None;
}

   /*
   * Initialize all Dihedral objects for Molecules of one Species.
   *
   * This functions assigns pointers to Atoms and Dihedral types ids within
   * a contiguous block of Dihedral objects, and sets a pointer in each 
   * Molecule to the first Dihedral in the associated block.
   */
   void Simulation::initializeSpeciesDihedrals(int iSpecies)
   {

      Species*  speciesPtr = 0;
      Molecule *moleculePtr = 0;
      Dihedral  *dihedralPtr = 0;
      Atom     *firstAtomPtr, *atom0Ptr, *atom1Ptr, *atom2Ptr, *atom3Ptr;
      int       iMol, iDihedral, atom0Id, atom1Id, atom2Id, atom3Id, type;
      int       capacity, nDihedral;

      speciesPtr = &species(iSpecies);
      capacity   = speciesPtr->capacity();
      nDihedral  = speciesPtr->nDihedral();

      // Initialize pointers before loop
      moleculePtr = &molecules_[firstMoleculeIds_[iSpecies]];
      dihedralPtr = &dihedrals_[firstDihedralIds_[iSpecies]];

      // Loop over molecules in Species
      for (iMol = 0; iMol < capacity; ++iMol) {

         firstAtomPtr = &(moleculePtr->atom(0));
         moleculePtr->setFirstDihedral(*dihedralPtr);
         moleculePtr->setNDihedral(nDihedral);

         if (nDihedral > 0) {

            // Create dihedrals for a molecule
            for (iDihedral = 0; iDihedral < nDihedral; ++iDihedral) {

               // Get local indices for atoms and dihedral type
               atom0Id  = speciesPtr->speciesDihedral(iDihedral).atomId(0);
               atom1Id  = speciesPtr->speciesDihedral(iDihedral).atomId(1);
               atom2Id  = speciesPtr->speciesDihedral(iDihedral).atomId(2);
               atom3Id  = speciesPtr->speciesDihedral(iDihedral).atomId(3);
               type     = speciesPtr->speciesDihedral(iDihedral).typeId();

               // Calculate atom pointers
               atom0Ptr = firstAtomPtr + atom0Id;
               atom1Ptr = firstAtomPtr + atom1Id;
               atom2Ptr = firstAtomPtr + atom2Id;
               atom3Ptr = firstAtomPtr + atom3Id;

               // Set fields of the Dihedral object
               dihedralPtr->setAtom(0, *atom0Ptr);
               dihedralPtr->setAtom(1, *atom1Ptr);
               dihedralPtr->setAtom(2, *atom2Ptr);
               dihedralPtr->setAtom(3, *atom3Ptr);
               dihedralPtr->setTypeId(type);

               ++dihedralPtr;
            }
         }
         ++moleculePtr;
      }

   }

   /// Evaluate end-to-end vectors of all chains.
   void EndtoEndXYZ::sample(long iStep) 
   { 
      if (isAtInterval(iStep))  {

         Molecule* moleculePtr;
         Vector    r1, r2, dR;
         double    dxSq, dySq, dzSq;
         int       i, j, nMolecule;

         dxSq = 0.0;
         dySq = 0.0; 
         dzSq = 0.0; 
         nMolecule = system().nMolecule(speciesId_);
         for (i = 0; i < system().nMolecule(speciesId_); i++) {
            moleculePtr = &system().molecule(speciesId_, i);

            // Construct map of molecule with no periodic boundary conditions
            positions_[0] = moleculePtr->atom(0).position();
            for (j = 1 ; j < nAtom_; j++) {
               r1 = moleculePtr->atom(j-1).position();
               r2 = moleculePtr->atom(j).position();
               system().boundary().distanceSq(r1, r2, dR);
               positions_[j]  = positions_[j-1];
               positions_[j] += dR;
            }

            dR.subtract(positions_[0], positions_[nAtom_-1]);
            dxSq += dR[0]*dR[0];
            dySq += dR[1]*dR[1];
            dzSq += dR[2]*dR[2];    
            
     
         }
         dxSq /= double(nMolecule);
         dySq /= double(nMolecule);
         dzSq /= double(nMolecule); 
      
         accumulatorX_.sample(dxSq, outputFileX_);
         accumulatorY_.sample(dySq, outputFileY_); 
         accumulatorZ_.sample(dzSq, outputFileZ_); 

      } // if isAtInterval

   }

static PyObject *contacts(PyObject *self, PyObject *args) {
  
  int mol1, frame1, mol2, frame2;
  PyObject *selected1, *selected2;
  float cutoff;
  if (!PyArg_ParseTuple(args, (char *)"iiO!iiO!f:atomselection.contacts",
        &mol1, &frame1, &PyTuple_Type, &selected1,
        &mol2, &frame2, &PyTuple_Type, &selected2,
        &cutoff))
    return NULL;
  VMDApp *app = get_vmdapp();
  AtomSel *sel1 = sel_from_py(mol1, frame1, selected1, app);
  AtomSel *sel2 = sel_from_py(mol2, frame2, selected2, app);
  if (!sel1 || !sel2) {
    delete sel1;
    delete sel2;
    return NULL;
  }
  const float *ts1 = sel1->coordinates(app->moleculeList);
  const float *ts2 = sel2->coordinates(app->moleculeList);
  if (!ts1 || !ts2) {
    PyErr_SetString(PyExc_ValueError, "No coordinates in selection");
    delete sel1;
    delete sel2;
    return NULL;
  }
  Molecule *mol = app->moleculeList->mol_from_id(mol1);

  GridSearchPair *pairlist = vmd_gridsearch3(
      ts1, sel1->num_atoms, sel1->on,
      ts2, sel2->num_atoms, sel2->on,
      cutoff, -1, (sel1->num_atoms + sel2->num_atoms) * 27);

  delete sel1;
  delete sel2;
  GridSearchPair *p, *tmp;
  PyObject *list1 = PyList_New(0);
  PyObject *list2 = PyList_New(0);
  for (p=pairlist; p != NULL; p=tmp) {
    // throw out pairs that are already bonded
    MolAtom *a1 = mol->atom(p->ind1);
    if (mol1 != mol2 || !a1->bonded(p->ind2)) {
      PyList_Append(list1, PyInt_FromLong(p->ind1));
      PyList_Append(list2, PyInt_FromLong(p->ind2));
    }
    tmp = p->next;
    free(p);
  }
  PyObject *result = PyList_New(2);
  PyList_SET_ITEM(result, 0, list1);
  PyList_SET_ITEM(result, 1, list2);
  return result;
}

   /*
   * Evaluate Rosenbluth weight, and add to accumulator.
   */
   void McChemicalPotential::sample(long iStep)
   {
      if (isAtInterval(iStep))  {

         Species* speciesPtr;
         Molecule* molPtr;
         Molecule::BondIterator bondIter;
         Atom* endPtr;
         double w;
         double rosenbluth = 1;
         double de;
         double e = 0;

         speciesPtr = &(simulation().species(speciesId_));

         // Pop a new molecule off the species reservoir
         molPtr = &(speciesPtr->reservoir().pop());
         system().addMolecule(*molPtr);

         // Loop over molecule growth trials
         for (int i = 0; i < nMoleculeTrial_; i++) {

            // Pick a random position for the first atom
            endPtr = &molPtr->atom(0);
            boundary().randomPosition(random(), endPtr->position());

            e = system().pairPotential().atomEnergy(*endPtr);
            rosenbluth = boltzmann(e);
            system().pairPotential().addAtom(*endPtr);

            for (molPtr->begin(bondIter); bondIter.notEnd(); ++bondIter) {
                addEndAtom(&(bondIter->atom(1)), &(bondIter->atom(0)), bondIter->typeId(), w, de);
                e += de;
                rosenbluth *= w;
                system().pairPotential().addAtom(bondIter->atom(1));
            }

            rosenbluth = rosenbluth / pow(nTrial_,molPtr->nAtom()-1);
            accumulator_.sample(rosenbluth, outputFile_);

            system().pairPotential().deleteAtom(*endPtr);
            for (molPtr->begin(bondIter); bondIter.notEnd(); ++bondIter) {
                system().pairPotential().deleteAtom(bondIter->atom(1));
            }
         }

         // Return additional molecule to reservoir
         system().removeMolecule(*molPtr);
         speciesPtr->reservoir().push(*molPtr);
      }

   }

void HydrogenTools::adjustHydrogens(Molecule &molecule, Adjustment adjustment)
{
  // This vector stores indices of hydrogens that need to be removed. Additions
  // are made first, followed by removals to keep indexing sane.
  std::vector<size_t> badHIndices;

  // Temporary container for calls to generateNewHydrogenPositions.
  std::vector<Vector3> newHPos;

  // Convert the adjustment option to a couple of booleans
  bool doAdd(adjustment == Add || adjustment == AddAndRemove);
  bool doRemove(adjustment == Remove || adjustment == AddAndRemove);

  // Limit to only the original atoms:
  const size_t numAtoms = molecule.atomCount();

  // Iterate through all atoms in the molecule, adding hydrogens as needed
  // and building up a list of hydrogens that should be removed.
  for (size_t atomIndex = 0; atomIndex < numAtoms; ++atomIndex) {
    const Atom atom(molecule.atom(atomIndex));
    int hDiff = valencyAdjustment(atom);
    // Add hydrogens:
    if (doAdd && hDiff > 0) {
      newHPos.clear();
      generateNewHydrogenPositions(atom, hDiff, newHPos);
      for (std::vector<Vector3>::const_iterator it = newHPos.begin(),
           itEnd = newHPos.end(); it != itEnd; ++it) {
        Atom newH(molecule.addAtom(1));
        newH.setPosition3d(*it);
        molecule.addBond(atom, newH, 1);
      }
    }
    // Add bad hydrogens to our list of hydrogens to remove:
    else if (doRemove && hDiff < 0) {
      extraHydrogenIndices(atom, -hDiff, badHIndices);
    }
  }

  // Remove dead hydrogens now. Remove them in reverse-index order to keep
  // indexing sane.
  if (doRemove && !badHIndices.empty()) {
    std::sort(badHIndices.begin(), badHIndices.end());
    std::vector<size_t>::iterator newEnd(std::unique(badHIndices.begin(),
                                                     badHIndices.end()));
    badHIndices.resize(std::distance(badHIndices.begin(), newEnd));
    for (std::vector<size_t>::const_reverse_iterator it = badHIndices.rbegin(),
         itEnd = badHIndices.rend(); it != itEnd; ++it) {
      molecule.removeAtom(*it);
    }
  }
}

   /*
   * Evaluate end-to-end vectors of all chains, add to ensemble.
   */
   void RingRouseAutoCorr::sample(long iStep) 
   { 
      if (isAtInterval(iStep))  {

         Molecule* moleculePtr;
         Vector    r1, r2, dR, atomPos;
         int       i, j;

         // Confirm that nMolecule has remained constant
         if (nMolecule_ != system().nMolecule(speciesId_)) {
            UTIL_THROW("Number of molecules has changed.");
         }

         // Loop over molecules
         for (i=0; i < nMolecule_; i++) {
            moleculePtr = &(system().molecule(speciesId_, i));

            // Retrace non-periodic shape and calculate coefficients
            atomPos = moleculePtr->atom(0).position();
            dR.multiply(atomPos, projector_[0]);
            data_[i] = dR;
            for (j = 1; j < nAtom_; j++) {
               r1 = moleculePtr->atom(j-1).position();
               r2 = moleculePtr->atom(j).position();
               system().boundary().distanceSq(r2, r1, dR);
               atomPos += dR;
               dR.multiply(atomPos, projector_[j]);
               data_[i] += dR;
            }
         }
      
         accumulator_.sample(data_);

      } // if isAtInterval

   }

TEST(RWMoleculeTest, MoleculeToRWMolecule)
{
  Molecule mol;
  typedef Molecule::AtomType Atom;
  typedef Molecule::BondType Bond;
  Atom a0 = mol.addAtom(1);
  Atom a1 = mol.addAtom(6);
  Atom a2 = mol.addAtom(9);
  Bond b0 = mol.addBond(a0, a2);
  a1.setPosition3d(Vector3(0, 6, 9));
  b0.setOrder(3);

  RWMolecule rwmol(mol, 0);
  EXPECT_EQ(rwmol.atomCount(), mol.atomCount());
  EXPECT_EQ(rwmol.bondCount(), mol.bondCount());
  EXPECT_EQ(rwmol.atom(2).atomicNumber(), mol.atom(2).atomicNumber());
  EXPECT_EQ(rwmol.bond(0).order(), mol.bond(0).order());
}

   /*
   * Evaluate end-to-end vectors of all chains, add to ensemble.
   */
   void G1MSD::sample(long iStep) 
   { 
      if (!isAtInterval(iStep)) return;

      // Confirm that nMolecule has remained constant, and nMolecule > 0.
      if (nMolecule_ <= 0) {
         UTIL_THROW("nMolecule <= 0");
      }
      if (nMolecule_ != system().nMolecule(speciesId_)) {
         UTIL_THROW("Number of molecules has changed.");
      }

      Vector     r;
      IntVector  shift;
      Molecule*  moleculePtr;
      Vector     lengths = system().boundary().lengths();
      int        i, j, k;
      int        iatom = 0;
      for (i = 0; i < nMolecule_; ++i) {
	 moleculePtr = &system().molecule(speciesId_, i);
	 for (j = 0 ; j < nAtom_; j++) {	 
            r = moleculePtr->atom(j).position();
            system().boundary().shift(r);

            // Compare current r to previous position, oldPositions_[iatom]
            system().boundary().distanceSq(r, oldPositions_[iatom], shift);

            // If this atom crossed a boundary, increment its shift vector
            shifts_[iatom] += shift;

            // Reconstruct true position
            for (k = 0; k < Dimension; ++k) {
               truePositions_[iatom][k] = r[k] + shifts_[iatom][k]*lengths[k];
            }

            // Store current position in box for comparison to next one
            oldPositions_[iatom] = r;
	    
	    iatom++;
	 }
      }
      accumulator_.sample(truePositions_);

   }

static PyObject *getbonds(PyObject *self, PyObject *args) {
  int molid;
  PyObject *atomlist;
   
  if (!PyArg_ParseTuple(args, (char *)"iO!:getbonds", &molid, 
                        &PyTuple_Type, &atomlist)) 
    return NULL;  // bad args
 
  Molecule *mol = get_vmdapp()->moleculeList->mol_from_id(molid);
  if (!mol) {
    PyErr_SetString(PyExc_ValueError, "molecule no longer exists");
    return NULL;
  }
  int num_atoms = mol->nAtoms;
  int num_selected = PyTuple_Size(atomlist);
  PyObject *newlist = PyList_New(num_selected);
  for (int i=0; i< num_selected; i++) {
    int id = PyInt_AsLong(PyTuple_GET_ITEM(atomlist, i));
    if (PyErr_Occurred()) {
      Py_DECREF(newlist);
      return NULL;
    }
    if (id < 0 || id >= num_atoms) {
      PyErr_SetString(PyExc_ValueError, (char *)"invalid atom id found");
      Py_DECREF(newlist);
      return NULL;
    }
    const MolAtom *atom = mol->atom(id);
    PyObject *bondlist = PyList_New(atom->bonds);
    for (int j=0; j<atom->bonds; j++) {
      PyList_SET_ITEM(bondlist, j, PyInt_FromLong(atom->bondTo[j]));
    }
    PyList_SET_ITEM(newlist, i, bondlist);
  }
  return newlist;
}