--- trunk/src/selection/DistanceFinder.cpp 2009/11/25 20:02:06 1390 +++ trunk/src/selection/DistanceFinder.cpp 2015/01/09 19:06:35 2052 @@ -1,5 +1,5 @@ /* - * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * Copyright (c) 2005, 2010 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a * non-exclusive, royalty free, license to use, modify and @@ -35,61 +35,314 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ +#ifdef IS_MPI +#include +#endif + #include "selection/DistanceFinder.hpp" #include "primitives/Molecule.hpp" -namespace OpenMD { +namespace OpenMD { + DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) { + nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies()); + nObjects_.push_back(info_->getNGlobalBonds()); + nObjects_.push_back(info_->getNGlobalBends()); + nObjects_.push_back(info_->getNGlobalTorsions()); + nObjects_.push_back(info_->getNGlobalInversions()); + nObjects_.push_back(info_->getNGlobalMolecules()); - nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies(); - stuntdoubles_.resize(nStuntDoubles_); + stuntdoubles_.resize(nObjects_[STUNTDOUBLE]); + bonds_.resize(nObjects_[BOND]); + bends_.resize(nObjects_[BEND]); + torsions_.resize(nObjects_[TORSION]); + inversions_.resize(nObjects_[INVERSION]); + molecules_.resize(nObjects_[MOLECULE]); SimInfo::MoleculeIterator mi; - Molecule* mol; Molecule::AtomIterator ai; - Atom* atom; Molecule::RigidBodyIterator rbIter; - RigidBody* rb; + Molecule::BondIterator bondIter; + Molecule::BendIterator bendIter; + Molecule::TorsionIterator torsionIter; + Molecule::InversionIterator inversionIter; + Molecule* mol; + Atom* atom; + RigidBody* rb; + Bond* bond; + Bend* bend; + Torsion* torsion; + Inversion* inversion; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + + molecules_[mol->getGlobalIndex()] = mol; + + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { stuntdoubles_[atom->getGlobalIndex()] = atom; } - - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { stuntdoubles_[rb->getGlobalIndex()] = rb; } - - } + for (bond = mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { + bonds_[bond->getGlobalIndex()] = bond; + } + for (bend = mol->beginBend(bendIter); bend != NULL; + bend = mol->nextBend(bendIter)) { + bends_[bend->getGlobalIndex()] = bend; + } + for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; + torsion = mol->nextTorsion(torsionIter)) { + torsions_[torsion->getGlobalIndex()] = torsion; + } + for (inversion = mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { + inversions_[inversion->getGlobalIndex()] = inversion; + } + } } - OpenMDBitSet DistanceFinder::find(const OpenMDBitSet& bs, RealType distance) { + SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance) { StuntDouble * center; Vector3d centerPos; Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); - OpenMDBitSet bsResult(nStuntDoubles_); + SelectionSet bsResult(nObjects_); assert(bsResult.size() == bs.size()); - - for (int i = bs.firstOnBit(); i != -1; i = bs.nextOnBit(i)) { + +#ifdef IS_MPI + int mol; + int proc; + RealType data[3]; + int worldRank; + MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); +#endif + + for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) { + if (stuntdoubles_[j]->isRigidBody()) { + RigidBody* rb = static_cast(stuntdoubles_[j]); + rb->updateAtoms(); + } + } + + SelectionSet bsTemp(nObjects_); + bsTemp = bs; + bsTemp.parallelReduce(); + + for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1; + i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) { + + // Now, if we own stuntdouble i, we can use the position, but in + // parallel, we'll need to let everyone else know what that + // position is! + +#ifdef IS_MPI + mol = info_->getGlobalMolMembership(i); + proc = info_->getMolToProc(mol); + + if (proc == worldRank) { + center = stuntdoubles_[i]; + centerPos = center->getPos(); + data[0] = centerPos.x(); + data[1] = centerPos.y(); + data[2] = centerPos.z(); + MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD); + } else { + MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD); + centerPos = Vector3d(data); + } +#else center = stuntdoubles_[i]; centerPos = center->getPos(); - for (int j = 0; j < stuntdoubles_.size(); ++j) { +#endif + + for (unsigned int j = 0; j < molecules_.size(); ++j) { + Vector3d r =centerPos - molecules_[j]->getCom(); + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[MOLECULE].setBitOn(j); + } + } + for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) { Vector3d r =centerPos - stuntdoubles_[j]->getPos(); currSnapshot->wrapVector(r); if (r.length() <= distance) { - bsResult.setBitOn(j); + bsResult.bitsets_[STUNTDOUBLE].setBitOn(j); } } + for (unsigned int j = 0; j < bonds_.size(); ++j) { + Vector3d loc = bonds_[j]->getAtomA()->getPos(); + loc += bonds_[j]->getAtomB()->getPos(); + loc = loc / 2.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[BOND].setBitOn(j); + } + } + for (unsigned int j = 0; j < bends_.size(); ++j) { + Vector3d loc = bends_[j]->getAtomA()->getPos(); + loc += bends_[j]->getAtomB()->getPos(); + loc += bends_[j]->getAtomC()->getPos(); + loc = loc / 3.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[BEND].setBitOn(j); + } + } + for (unsigned int j = 0; j < torsions_.size(); ++j) { + Vector3d loc = torsions_[j]->getAtomA()->getPos(); + loc += torsions_[j]->getAtomB()->getPos(); + loc += torsions_[j]->getAtomC()->getPos(); + loc += torsions_[j]->getAtomD()->getPos(); + loc = loc / 4.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[TORSION].setBitOn(j); + } + } + for (unsigned int j = 0; j < inversions_.size(); ++j) { + Vector3d loc = inversions_[j]->getAtomA()->getPos(); + loc += inversions_[j]->getAtomB()->getPos(); + loc += inversions_[j]->getAtomC()->getPos(); + loc += inversions_[j]->getAtomD()->getPos(); + loc = loc / 4.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[INVERSION].setBitOn(j); + } + } } - return bsResult; } + + + SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance, int frame ) { + StuntDouble * center; + Vector3d centerPos; + Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame); + SelectionSet bsResult(nObjects_); + assert(bsResult.size() == bs.size()); +#ifdef IS_MPI + int mol; + int proc; + RealType data[3]; + int worldRank; + MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); +#endif + + for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) { + if (stuntdoubles_[j]->isRigidBody()) { + RigidBody* rb = static_cast(stuntdoubles_[j]); + rb->updateAtoms(frame); + } + } + + SelectionSet bsTemp(nObjects_); + bsTemp = bs; + bsTemp.parallelReduce(); + + for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1; + i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) { + + // Now, if we own stuntdouble i, we can use the position, but in + // parallel, we'll need to let everyone else know what that + // position is! + +#ifdef IS_MPI + mol = info_->getGlobalMolMembership(i); + proc = info_->getMolToProc(mol); + + if (proc == worldRank) { + center = stuntdoubles_[i]; + centerPos = center->getPos(frame); + data[0] = centerPos.x(); + data[1] = centerPos.y(); + data[2] = centerPos.z(); + MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD); + } else { + MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD); + centerPos = Vector3d(data); + } +#else + center = stuntdoubles_[i]; + centerPos = center->getPos(frame); +#endif + for (unsigned int j = 0; j < molecules_.size(); ++j) { + Vector3d r =centerPos - molecules_[j]->getCom(frame); + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[MOLECULE].setBitOn(j); + } + } + + for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) { + Vector3d r =centerPos - stuntdoubles_[j]->getPos(frame); + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[STUNTDOUBLE].setBitOn(j); + } + } + for (unsigned int j = 0; j < bonds_.size(); ++j) { + Vector3d loc = bonds_[j]->getAtomA()->getPos(frame); + loc += bonds_[j]->getAtomB()->getPos(frame); + loc = loc / 2.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[BOND].setBitOn(j); + } + } + for (unsigned int j = 0; j < bends_.size(); ++j) { + Vector3d loc = bends_[j]->getAtomA()->getPos(frame); + loc += bends_[j]->getAtomB()->getPos(frame); + loc += bends_[j]->getAtomC()->getPos(frame); + loc = loc / 3.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[BEND].setBitOn(j); + } + } + for (unsigned int j = 0; j < torsions_.size(); ++j) { + Vector3d loc = torsions_[j]->getAtomA()->getPos(frame); + loc += torsions_[j]->getAtomB()->getPos(frame); + loc += torsions_[j]->getAtomC()->getPos(frame); + loc += torsions_[j]->getAtomD()->getPos(frame); + loc = loc / 4.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[TORSION].setBitOn(j); + } + } + for (unsigned int j = 0; j < inversions_.size(); ++j) { + Vector3d loc = inversions_[j]->getAtomA()->getPos(frame); + loc += inversions_[j]->getAtomB()->getPos(frame); + loc += inversions_[j]->getAtomC()->getPos(frame); + loc += inversions_[j]->getAtomD()->getPos(frame); + loc = loc / 4.0; + Vector3d r = centerPos - loc; + currSnapshot->wrapVector(r); + if (r.length() <= distance) { + bsResult.bitsets_[INVERSION].setBitOn(j); + } + } + } + return bsResult; + } }