--- branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/11/09 23:11:39 1722 +++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/11/13 05:08:12 1738 @@ -37,24 +37,82 @@ SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0 namespace oopse { -SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), - nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), sman_(NULL){ +SimInfo::SimInfo(const std::vector >& molStampPairs, + ForceField* ff, Globals* globals) : + forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0), + nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0), + nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL), + fortranInitialized_(false) { + std::vector >::iterator i; + int nCutoffAtoms; // number of atoms belong to cutoff groups + int ngroups; //total cutoff groups defined in meta-data file + MoleculeStamp* molStamp; + int nMolWithSameStamp; + CutoffGroupStamp* cgStamp; + int nAtomsIngroups; + int nCutoffGroupsInStamp; + + nGlobalAtoms_ = 0; + ngroups = 0; + + for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { + molStamp = i->first; + nMolWithSameStamp = i->second; + + addMoleculeStamp(molStamp, nMolWithSameStamp); + + nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; + + nAtomsIngroups = 0; + nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); + + for (int j=0; j < nCutoffGroupsInStamp; j++) { + cgStamp = molStamp->getCutoffGroup(j); + nAtomsIngroups += cgStamp->getNMembers(); + } + + ngroups += *nMolWithSameStamp; + nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp; + } + + //every free atom (atom does not belong to cutoff groups) is a cutoff group + //therefore the total number of cutoff groups in the system is 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 + nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + ngroups; + + //initialize globalGroupMembership_, every element of this array will be 0 + globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0); + + nGlobalMols_ = molStampIds_.size(); + +#ifdef IS_MPI + molToProcMap_.resize(nGlobalMols_); +#endif + } SimInfo::~SimInfo() { - MemoryUtils::deleteVectorOfPointer(molecules_); + //MemoryUtils::deleteVectorOfPointer(molecules_); + + MemoryUtils::deleteVectorOfPointer(moleculeStamps_); + delete sman_; + delete globals_; + delete forceField_; } bool SimInfo::addMolecule(Molecule* mol) { - std::vector::iterator i; - i = std::find(molecules_.begin(), molecules_.end(), mol); + MoleculeIterator i; + + i = molecules_.find(mol->getGlobalIndex()); if (i != molecules_.end() ) { - molecules_.push_back(mol); + molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); + nAtoms_ += mol->getNAtoms(); nBonds_ += mol->getNBonds(); nBends_ += mol->getNBends(); @@ -71,11 +129,13 @@ bool SimInfo::removeMolecule(Molecule* mol) { } bool SimInfo::removeMolecule(Molecule* mol) { - std::vector::iterator i; - i = std::find(molecules_.begin(), molecules_.end(), mol); + MoleculeIterator i; + i = molecules_.find(mol->getGlobalIndex()); if (i != molecules_.end() ) { - molecules_.push_back(mol); + + assert(mol == i->second); + nAtoms_ -= mol->getNAtoms(); nBonds_ -= mol->getNBonds(); nBends_ -= mol->getNBends(); @@ -85,6 +145,10 @@ bool SimInfo::removeMolecule(Molecule* mol) { nCutoffGroups_ -= mol->getNCutoffGroups(); nConstraints_ -= mol->getNConstraints(); + molecules_.erase(mol->getGlobalIndex()); + + delete mol; + return true; } else { return false; @@ -94,20 +158,20 @@ Molecule* SimInfo::beginMolecule(std::vector::iterator& i) { +Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { i = molecules_.begin(); - return i == molecules_.end() ? NULL : *i; + return i == molecules_.end() ? NULL : i->second; } -Molecule* SimInfo::nextMolecule(std::vector::iterator& i) { +Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { ++i; - return i == molecules_.end() ? NULL : *i; + return i == molecules_.end() ? NULL : i->second; } void SimInfo::calcNdf() { int ndf_local; - std::vector::iterator i; + MoleculeIterator i; std::vector::iterator j; Molecule* mol; StuntDouble* integrableObject; @@ -140,16 +204,16 @@ void SimInfo::calcNdf() { ndf_ = ndf_local; #endif - // nZconstraints is global, as are the 3 COM translations for the + // nZconstraints_ is global, as are the 3 COM translations for the // entire system: - ndf_ = ndf_ - 3 - nZconstraints; + ndf_ = ndf_ - 3 - nZconstraints_; } void SimInfo::calcNdfRaw() { int ndfRaw_local; - std::vector::iterator i; + MoleculeIterator i; std::vector::iterator j; Molecule* mol; StuntDouble* integrableObject; @@ -193,7 +257,7 @@ void SimInfo::calcNdfTrans() { ndfTrans_ = ndfTrans_local; #endif - ndfTrans_ = ndfTrans_ - 3 - nZconstraints; + ndfTrans_ = ndfTrans_ - 3 - nZconstraints_; } @@ -287,4 +351,415 @@ void SimInfo::removeExcludePairs(Molecule* mol) { } +void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { + int curStampId; + + //index from 0 + curStampId = molStampIds_.size(); + + moleculeStamps_.push_back(molStamp); + molStampIds_.insert(molStampIds_.end(), nmol, curStampId) +} + +void SimInfo::update() { + + setupSimType(); + +#ifdef IS_MPI + setupFortranParallel(); +#endif + + setupFortranSim(); + + setupCutoff(); + + //notify fortran whether reaction field is used or not. It is deprecated now + //int isError = 0; + //initFortranFF( &useReactionField, &isError ); + + //if(isError){ + // sprintf( painCave.errMsg, + // "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" ); + // painCave.isFatal = 1; + // simError(); + //} + + calcNdf(); + calcNdfRaw(); + calcNdfTrans(); + + fortranInitialized_ = true; +} + +std::set SimInfo::getUniqueAtomTypes() { + typename SimInfo::MoleculeIterator mi; + Molecule* mol; + typename Molecule::AtomIterator ai; + Atom* atom; + std::set atomTypes; + + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + atomTypes.insert(atom->getAtomType()); + } + + } + + return atomTypes; +} + +void SimInfo::setupSimType() { + std::set::iterator i; + std::set atomTypes; + atomTypes = getUniqueAtomTypes(); + + int useLennardJones = 0; + int useElectrostatic = 0; + int useEAM = 0; + int useCharge = 0; + int useDirectional = 0; + int useDipole = 0; + int useGayBerne = 0; + int useSticky = 0; + int useShape = 0; + int useFLARB = 0; //it is not in AtomType yet + int useDirectionalAtom = 0; + int useElectrostatics = 0; + //usePBC and useRF are from globals + bool usePBC = globals_->getPBC(); + bool useRF = globals_->getUseRF(); + + //loop over all of the atom types + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + useLennardJones |= i->isLennardJones(); + useElectrostatic |= i->isElectrostatic(); + useEAM |= i->isEAM(); + useCharge |= i->isCharge(); + useDirectional |= i->isDirectional(); + useDipole |= i->isDipole(); + useGayBerne |= i->isGayBerne(); + useSticky |= i->isSticky(); + useShape |= i->isShape(); + } + + if (useSticky || useDipole || useGayBerne || useShape) { + useDirectionalAtom = 1; + } + + if (useCharge || useDipole) { + useElectrostatics = 1; + } + +#ifdef IS_MPI + int temp; + + temp = usePBC; + MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useDirectionalAtom; + MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useLennardJones; + MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useElectrostatics; + MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useCharge; + MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useDipole; + MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useSticky; + MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useGayBerne; + MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useEAM; + MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useShape; + MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useFLARB; + MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useRF; + MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + +#endif + + fInfo_.SIM_uses_PBC = usePBC; + fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; + fInfo_.SIM_uses_LennardJones = useLennardJones; + fInfo_.SIM_uses_Electrostatics = useElectrostatics; + fInfo_.SIM_uses_Charges = useCharge; + fInfo_.SIM_uses_Dipoles = useDipole; + fInfo_.SIM_uses_Sticky = useSticky; + fInfo_.SIM_uses_GayBerne = useGayBerne; + fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_Shapes = useShape; + fInfo_.SIM_uses_FLARB = useFLARB; + fInfo_.SIM_uses_RF = useRF; + + if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { + fInfo_.dielect = dielectric; + } else { + fInfo_.dielect = 0.0; + } + +} + +void SimInfo::setupFortranSim() { + int isError; + int nExclude; + std::vector fortranGlobalGroupMembership; + + nExclude = exclude_.getSize(); + isError = 0; + + //globalGroupMembership_ is filled by SimCreator + for (int i = 0; i < nGlobalAtoms_; i++) { + fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); + } + + //calculate mass ratio of cutoff group + std::vector mfact; + typename SimInfo::MoleculeIterator mi; + Molecule* mol; + typename Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; + typename Molecule::AtomIterator ai; + Atom* atom; + double totalMass; + + //to avoid memory reallocation, reserve enough space for mfact + mfact.reserve(getNCutoffGroups()); + + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { + + totalMass = cg->getMass(); + for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { + mfact.push_back(atom->getMass()/totalMass); + } + + } + } + + //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) + std::vector identArray; + + //to avoid memory reallocation, reserve enough space identArray + identArray.reserve(getNAtoms()); + + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + identArray.push_back(atom->getIdent()); + } + } + + //fill molMembershipArray + //molMembershipArray is filled by SimCreator + std::vector molMembershipArray(nGlobalAtoms_); + for (int i = 0; i < nGlobalAtoms_; i++) { + molMembershipArray.push_back(globalMolMembership_[i] + 1); + } + + //setup fortran simulation + //gloalExcludes and molMembershipArray should go away (They are never used) + //why the hell fortran need to know molecule? + //OOPSE = Object-Obfuscated Parallel Simulation Engine + + setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(), + &nGlobalExcludes, globalExcludes, molMembershipArray, + &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); + + if( isError ){ + + sprintf( painCave.errMsg, + "There was an error setting the simulation information in fortran.\n" ); + painCave.isFatal = 1; + painCave.severity = OOPSE_ERROR; + simError(); + } + +#ifdef IS_MPI + sprintf( checkPointMsg, + "succesfully sent the simulation information to fortran.\n"); + MPIcheckPoint(); +#endif // is_mpi +} + + +#ifdef IS_MPI +void SimInfo::setupFortranParallel() { + + //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex + std::vector localToGlobalAtomIndex(getNAtoms(), 0); + std::vector localToGlobalCutoffGroupIndex; + typename SimInfo::MoleculeIterator mi; + typename Molecule::AtomIterator ai; + typename Molecule::CutoffGroupIterator ci; + Molecule* mol; + Atom* atom; + CutoffGroup* cg; + mpiSimData parallelData; + int isError; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + //local index(index in DataStorge) of atom is important + for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; + } + + //local index of cutoff group is trivial, it only depends on the order of travesing + for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { + localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); + } + + } + + //fill up mpiSimData struct + parallelData.nMolGlobal = getNGlobalMolecules(); + parallelData.nMolLocal = getNMolecules(); + parallelData.nAtomsGlobal = getNGlobalAtoms(); + parallelData.nAtomsLocal = getNAtoms(); + parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); + parallelData.nGroupsLocal = getNCutoffGroups(); + parallelData.myNode = worldRank; + MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors)); + + //pass mpiSimData struct and index arrays to fortran + setFsimParallel(parallelData, &(parallelData->nAtomsLocal), + &localToGlobalAtomIndex[0], &(parallelData->nGroupsLocal), + &localToGlobalCutoffGroupIndex[0], &isError); + + if (isError) { + sprintf(painCave.errMsg, + "mpiRefresh errror: fortran didn't like something we gave it.\n"); + painCave.isFatal = 1; + simError(); + } + + sprintf(checkPointMsg, " mpiRefresh successful.\n"); + MPIcheckPoint(); + + +} + +#endif + +double SimInfo::calcMaxCutoffRadius() { + + + std::vector atomTypes; + std::vector::iterator i; + std::vector cutoffRadius; + + //get the unique atom types + atomTypes = getUniqueAtomTypes(); + + //query the max cutoff radius among these atom types + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); + } + + double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end()); +#ifdef IS_MPI + //pick the max cutoff radius among the processors +#endif + + return maxCutoffRadius; +} + +void SimInfo::setupCutoff() { + double rcut_; //cutoff radius + double rsw_; //switching radius + + if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + + if (!globals_->haveRcut()){ + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the cutoffRadius.\n" + "\tOOPSE will use a default value of 15.0 angstroms" + "\tfor the cutoffRadius.\n"); + painCave.isFatal = 0; + simError(); + rcut_ = 15.0; + } else{ + rcut_ = globals_->getRcut(); + } + + if (!globals_->haveRsw()){ + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\n" + "\t0.95 * cutoffRadius for the switchingRadius\n"); + painCave.isFatal = 0; + simError(); + rsw_ = 0.95 * rcut_; + } else{ + rsw_ = globals_->getRsw(); + } + + } else { + // if charge, dipole or reaction field is not used and the cutofff radius is not specified in + //meta-data file, the maximum cutoff radius calculated from forcefiled will be used + + if (globals_->haveRcut()) { + rcut_ = globals_->getRcut(); + } else { + //set cutoff radius to the maximum cutoff radius based on atom types in the whole system + rcut_ = calcMaxCutoffRadius(); + } + + if (globals_->haveRsw()) { + rsw_ = globals_->getRsw() + } else { + rsw_ = rcut_; + } + + } + + double rnblist = rcut_ + 1; // skin of neighbor list + + //Pass these cutoff radius etc. to fortran. This function should be called once and only once + notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); +} + +void SimInfo::addProperty(GenericData* genData) { + properties_.addProperty(genData); +} + +void SimInfo::removeProperty(const std::string& propName) { + properties_.removeProperty(propName); +} + +void SimInfo::clearProperties() { + properties_.clearProperties(); +} + +std::vector SimInfo::getPropertyNames() { + return properties_.getPropertyNames(); +} + +std::vector SimInfo::getProperties() { + return properties_.getProperties(); +} + +GenericData* SimInfo::getPropertyByName(const std::string& propName) { + return properties_.getPropertyByName(propName); +} + + +std::ostream& operator <<(ostream& o, SimInfo& info) { + + return o; +} + }//end namespace oopse