--- branches/development/src/brains/Thermo.cpp 2013/02/20 15:39:39 1850 +++ branches/development/src/brains/Thermo.cpp 2013/05/07 19:09:54 1870 @@ -233,7 +233,7 @@ namespace OpenMD { kinetic *= 0.5; eTemp = (2.0 * kinetic) / - (info_->getNFluctuatingCharges() * PhysicalConstants::kb ); + (info_->getNFluctuatingCharges() * PhysicalConstants::kb ); snap->setElectronicTemperature(eTemp); } @@ -619,8 +619,8 @@ namespace OpenMD { } /** - * Return intertia tensor for entire system and angular momentum - * Vector. + * \brief Return inertia tensor for entire system and angular momentum + * Vector. * * * @@ -700,7 +700,68 @@ namespace OpenMD { return; } + + + Mat3x3d Thermo::getBoundingBox(){ + + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + + if (!(snap->hasBoundingBox)) { + + SimInfo::MoleculeIterator i; + Molecule::RigidBodyIterator ri; + Molecule::AtomIterator ai; + Molecule* mol; + RigidBody* rb; + Atom* atom; + Vector3d pos, bMax, bMin; + int index = 0; + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + //change the positions of atoms which belong to the rigidbodies + for (rb = mol->beginRigidBody(ri); rb != NULL; + rb = mol->nextRigidBody(ri)) { + rb->updateAtoms(); + } + + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + pos = atom->getPos(); + + if (index == 0) { + bMax = pos; + bMin = pos; + } else { + for (int i = 0; i < 3; i++) { + bMax[i] = max(bMax[i], pos[i]); + bMin[i] = min(bMin[i], pos[i]); + } + } + index++; + } + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMax[0], 3, MPI::REALTYPE, + MPI::MAX); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMin[0], 3, MPI::REALTYPE, + MPI::MIN); +#endif + Mat3x3d bBox = Mat3x3d(0.0); + for (int i = 0; i < 3; i++) { + bBox(i,i) = bMax[i] - bMin[i]; + } + snap->setBoundingBox(bBox); + } + + return snap->getBoundingBox(); + } + + // Returns the angular momentum of the system Vector3d Thermo::getAngularMomentum(){ Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); @@ -896,10 +957,15 @@ namespace OpenMD { // Compute surface Mesh surfaceMesh_->computeHull(localSites_); snap->setHullVolume(surfaceMesh_->getVolume()); + + delete surfaceMesh_; } + return snap->getHullVolume(); #else return 0.0; #endif } + + }