--- branches/development/src/brains/SimInfo.cpp 2011/05/26 13:55:04 1569 +++ branches/development/src/brains/SimInfo.cpp 2011/08/04 20:04:35 1601 @@ -125,13 +125,8 @@ namespace OpenMD { //equal to the total number of atoms minus number of atoms belong to //cutoff group defined in meta-data file plus the number of cutoff //groups defined in meta-data file - std::cerr << "nGA = " << nGlobalAtoms_ << "\n"; - std::cerr << "nCA = " << nCutoffAtoms << "\n"; - std::cerr << "nG = " << nGroups << "\n"; nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; - - std::cerr << "nGCG = " << nGlobalCutoffGroups_ << "\n"; //every free atom (atom does not belong to rigid bodies) is an //integrable object therefore the total number of integrable objects @@ -274,7 +269,26 @@ namespace OpenMD { #endif return fdf_; } + + unsigned int SimInfo::getNLocalCutoffGroups(){ + int nLocalCutoffAtoms = 0; + Molecule* mol; + MoleculeIterator mi; + CutoffGroup* cg; + Molecule::CutoffGroupIterator ci; + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { + nLocalCutoffAtoms += cg->getNumAtom(); + + } + } + + return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; + } + void SimInfo::calcNdfRaw() { int ndfRaw_local; @@ -680,17 +694,18 @@ namespace OpenMD { Atom* atom; set atomTypes; - for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { atomTypes.insert(atom->getAtomType()); } } - + #ifdef IS_MPI // loop over the found atom types on this processor, and add their // numerical idents to a vector: - + vector foundTypes; set::iterator i; for (i = atomTypes.begin(); i != atomTypes.end(); ++i) @@ -699,41 +714,50 @@ namespace OpenMD { // count_local holds the number of found types on this processor int count_local = foundTypes.size(); - // count holds the total number of found types on all processors - // (some will be redundant with the ones found locally): - int count; - MPI::COMM_WORLD.Allreduce(&count_local, &count, 1, MPI::INT, MPI::SUM); + int nproc = MPI::COMM_WORLD.Get_size(); - // create a vector to hold the globally found types, and resize it: - vector ftGlobal; - ftGlobal.resize(count); - vector counts; + // we need arrays to hold the counts and displacement vectors for + // all processors + vector counts(nproc, 0); + vector disps(nproc, 0); - int nproc = MPI::COMM_WORLD.Get_size(); - counts.resize(nproc); - vector disps; - disps.resize(nproc); + // fill the counts array + MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], + 1, MPI::INT); + + // use the processor counts to compute the displacement array + disps[0] = 0; + int totalCount = counts[0]; + for (int iproc = 1; iproc < nproc; iproc++) { + disps[iproc] = disps[iproc-1] + counts[iproc-1]; + totalCount += counts[iproc]; + } - // now spray out the foundTypes to all the other processors: + // we need a (possibly redundant) set of all found types: + vector ftGlobal(totalCount); + // now spray out the foundTypes to all the other processors: MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, - &ftGlobal[0], &counts[0], &disps[0], MPI::INT); + &ftGlobal[0], &counts[0], &disps[0], + MPI::INT); + vector::iterator j; + // foundIdents is a stl set, so inserting an already found ident // will have no effect. set foundIdents; - vector::iterator j; + for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) foundIdents.insert((*j)); // now iterate over the foundIdents and get the actual atom types // that correspond to these: set::iterator it; - for (it = foundIdents.begin(); it != foundIdents.end(); ++it) + for (it = foundIdents.begin(); it != foundIdents.end(); ++it) atomTypes.insert( forceField_->getAtomType((*it)) ); #endif - + return atomTypes; } @@ -745,7 +769,7 @@ namespace OpenMD { if ( simParams_->getAccumulateBoxDipole() ) { calcBoxDipole_ = true; } - + set::iterator i; set atomTypes; atomTypes = getSimulatedAtomTypes(); @@ -758,18 +782,28 @@ namespace OpenMD { usesMetallic |= (*i)->isMetal(); usesDirectional |= (*i)->isDirectional(); } - + #ifdef IS_MPI int temp; temp = usesDirectional; MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + temp = usesMetallic; MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + temp = usesElectrostatic; MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); +#else + + usesDirectionalAtoms_ = usesDirectional; + usesMetallicAtoms_ = usesMetallic; + usesElectrostaticAtoms_ = usesElectrostatic; + #endif + + requiresPrepair_ = usesMetallicAtoms_ ? true : false; + requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; + requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; } @@ -824,9 +858,16 @@ namespace OpenMD { Atom* atom; RealType totalMass; - //to avoid memory reallocation, reserve enough space for massFactors_ + /** + * The mass factor is the relative mass of an atom to the total + * mass of the cutoff group it belongs to. By default, all atoms + * are their own cutoff groups, and therefore have mass factors of + * 1. We need some special handling for massless atoms, which + * will be treated as carrying the entire mass of the cutoff + * group. + */ massFactors_.clear(); - massFactors_.reserve(getNCutoffGroups()); + massFactors_.resize(getNAtoms(), 1.0); for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { for (cg = mol->beginCutoffGroup(ci); cg != NULL; @@ -835,10 +876,10 @@ namespace OpenMD { totalMass = cg->getMass(); for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { // Check for massless groups - set mfact to 1 if true - if (totalMass != 0) - massFactors_.push_back(atom->getMass()/totalMass); + if (totalMass != 0) + massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; else - massFactors_.push_back( 1.0 ); + massFactors_[atom->getLocalIndex()] = 1.0; } } } @@ -865,14 +906,6 @@ namespace OpenMD { int* oneThreeList = oneThreeInteractions_.getPairList(); int* oneFourList = oneFourInteractions_.getPairList(); - //setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray_[0], - // &nExclude, excludeList, - // &nOneTwo, oneTwoList, - // &nOneThree, oneThreeList, - // &nOneFour, oneFourList, - // &molMembershipArray[0], &mfact[0], &nCutoffGroups_, - // &fortranGlobalGroupMembership[0], &isError); - topologyDone_ = true; }