--- trunk/src/brains/SimCreator.cpp 2006/09/25 22:08:33 1051 +++ trunk/src/brains/SimCreator.cpp 2014/02/28 13:25:13 1971 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,15 +28,30 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** * @file SimCreator.cpp * @author tlin * @date 11/03/2004 - * @time 13:51am * @version 1.0 */ + +#ifdef IS_MPI +#include "mpi.h" +#include "math/ParallelRandNumGen.hpp" +#endif + #include #include #include @@ -55,7 +61,7 @@ #include "brains/SimCreator.hpp" #include "brains/SimSnapshotManager.hpp" #include "io/DumpReader.hpp" -#include "UseTheForce/ForceFieldFactory.hpp" +#include "brains/ForceField.hpp" #include "utils/simError.h" #include "utils/StringUtils.hpp" #include "math/SeqRandNumGen.hpp" @@ -75,13 +81,18 @@ #include "antlr/NoViableAltForCharException.hpp" #include "antlr/NoViableAltException.hpp" -#ifdef IS_MPI -#include "math/ParallelRandNumGen.hpp" -#endif +#include "types/DirectionalAdapter.hpp" +#include "types/MultipoleAdapter.hpp" +#include "types/EAMAdapter.hpp" +#include "types/SuttonChenAdapter.hpp" +#include "types/PolarizableAdapter.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/FluctuatingChargeAdapter.hpp" -namespace oopse { + +namespace OpenMD { - Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){ + Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){ Globals* simParams = NULL; try { @@ -90,34 +101,43 @@ namespace oopse { #ifdef IS_MPI int streamSize; const int masterNode = 0; - int commStatus; + if (worldRank == masterNode) { -#endif - + MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); +#endif SimplePreprocessor preprocessor; - preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); + preprocessor.preprocess(rawMetaDataStream, filename, + startOfMetaDataBlock, ppStream); #ifdef IS_MPI - //brocasting the stream size + //broadcasting the stream size streamSize = ppStream.str().size() +1; - commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); + MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); + MPI_Bcast(static_cast(const_cast(ppStream.str().c_str())), + streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); - commStatus = MPI_Bcast(static_cast(const_cast(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); - - + // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); + // MPI::COMM_WORLD.Bcast(static_cast(const_cast(ppStream.str().c_str())), + // streamSize, MPI::CHAR, masterNode); + } else { + + MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); + //get stream size - commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); - + MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); char* buf = new char[streamSize]; assert(buf); //receive file content - commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); - - ppStream.str(buf); - delete buf; + MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); + ppStream.str(buf); + delete [] buf; } #endif // Create a scanner that reads from the input stream @@ -139,13 +159,11 @@ namespace oopse { parser.initializeASTFactory(factory); parser.setASTFactory(&factory); parser.mdfile(); - // Create a tree parser that reads information into Globals MDTreeParser treeParser; treeParser.initializeASTFactory(factory); treeParser.setASTFactory(&factory); simParams = treeParser.walkTree(parser.getAST()); - } @@ -215,7 +233,7 @@ namespace oopse { painCave.isFatal = 1; simError(); } - catch (OOPSEException& e) { + catch (OpenMDException& e) { sprintf(painCave.errMsg, "%s\n", e.getMessage().c_str()); @@ -230,27 +248,48 @@ namespace oopse { simError(); } + simParams->setMDfileVersion(mdFileVersion); return simParams; } SimInfo* SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) { - + const int bufferSize = 65535; char buffer[bufferSize]; int lineNo = 0; std::string mdRawData; int metaDataBlockStart = -1; int metaDataBlockEnd = -1; - int i; - int mdOffset; + int i, j; + streamoff mdOffset; + int mdFileVersion; + // Create a string for embedding the version information in the MetaData + std::string version; + version.assign("## Last run using OpenMD Version: "); + version.append(OPENMD_VERSION_MAJOR); + version.append("."); + version.append(OPENMD_VERSION_MINOR); + + std::string svnrev; + //convert a macro from compiler to a string in c++ + STR_DEFINE(svnrev, SVN_REV ); + version.append(" Revision: "); + // If there's no SVN revision, just call this the RELEASE revision. + if (!svnrev.empty()) { + version.append(svnrev); + } else { + version.append("RELEASE"); + } + #ifdef IS_MPI const int masterNode = 0; if (worldRank == masterNode) { #endif - std::ifstream mdFile_(mdFileName.c_str()); + std::ifstream mdFile_; + mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary); if (mdFile_.fail()) { sprintf(painCave.errMsg, @@ -263,15 +302,41 @@ namespace oopse { mdFile_.getline(buffer, bufferSize); ++lineNo; std::string line = trimLeftCopy(buffer); - i = CaseInsensitiveFind(line, "(i) == string::npos) { + // try the older file strings to see if that works: + i = CaseInsensitiveFind(line, "(i) == string::npos) { + // still no luck! sprintf(painCave.errMsg, - "SimCreator: File: %s is not an OOPSE file!\n", + "SimCreator: File: %s is not a valid OpenMD file!\n", mdFileName.c_str()); painCave.isFatal = 1; simError(); } + + // found the correct opening string, now try to get the file + // format version number. + StringTokenizer tokenizer(line, "=<> \t\n\r"); + std::string fileType = tokenizer.nextToken(); + toUpper(fileType); + + mdFileVersion = 0; + + if (fileType == "OPENMD") { + while (tokenizer.hasMoreTokens()) { + std::string token(tokenizer.nextToken()); + toUpper(token); + if (token == "VERSION") { + mdFileVersion = tokenizer.nextTokenAsInt(); + break; + } + } + } + //scan through the input stream and find MetaData tag while(mdFile_.getline(buffer, bufferSize)) { ++lineNo; @@ -312,12 +377,23 @@ namespace oopse { mdRawData.clear(); + bool foundVersion = false; + for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) { mdFile_.getline(buffer, bufferSize); - mdRawData += buffer; + std::string line = trimLeftCopy(buffer); + j = CaseInsensitiveFind(line, "## Last run using OpenMD Version"); + if (static_cast(j) != string::npos) { + foundVersion = true; + mdRawData += version; + } else { + mdRawData += buffer; + } mdRawData += "\n"; } - + + if (!foundVersion) mdRawData += version + "\n"; + mdFile_.close(); #ifdef IS_MPI @@ -327,12 +403,12 @@ namespace oopse { std::stringstream rawMetaDataStream(mdRawData); //parse meta-data file - Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1); + Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion, + metaDataBlockStart + 1); //create the force field - ForceField * ff = ForceFieldFactory::getInstance() - ->createForceField(simParams->getForceField()); - + ForceField * ff = new ForceField(simParams->getForceField()); + if (ff == NULL) { sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n", @@ -342,7 +418,6 @@ namespace oopse { } if (simParams->haveForceFieldFileName()) { - std::cout<< simParams->getForceFieldFileName() << "\n"; ff->setForceFieldFileName(simParams->getForceFieldFileName()); } @@ -366,7 +441,6 @@ namespace oopse { } ff->parse(forcefieldFileName); - ff->setFortranForceOptions(); //create SimInfo SimInfo * info = new SimInfo(ff, simParams); @@ -384,10 +458,13 @@ namespace oopse { //create the molecules createMolecules(info); - + //find the storage layout + + int storageLayout = computeStorageLayout(info); + //allocate memory for DataStorage(circular reference, need to //break it) - info->setSnapshotManager(new SimSnapshotManager(info)); + info->setSnapshotManager(new SimSnapshotManager(info, storageLayout)); //set the global index of atoms, rigidbodies and cutoffgroups //(only need to be set once, the global index will never change @@ -396,21 +473,20 @@ namespace oopse { //responsibility to LocalIndexManager. setGlobalIndex(info); - //Although addExcludePairs is called inside SimInfo's addMolecule + //Although addInteractionPairs is called inside SimInfo's addMolecule //method, at that point atoms don't have the global index yet //(their global index are all initialized to -1). Therefore we - //have to call addExcludePairs explicitly here. A way to work + //have to call addInteractionPairs explicitly here. A way to work //around is that we can determine the beginning global indices of //atoms before they get created. SimInfo::MoleculeIterator mi; Molecule* mol; for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { - info->addExcludePairs(mol); + info->addInteractionPairs(mol); } if (loadInitCoords) loadCoordinates(info, mdFileName); - return info; } @@ -445,27 +521,16 @@ namespace oopse { #ifdef IS_MPI void SimCreator::divideMolecules(SimInfo *info) { - RealType numerator; - RealType denominator; - RealType precast; - RealType x; - RealType y; RealType a; - int old_atoms; - int add_atoms; - int new_atoms; - int nTarget; - int done; - int i; - int j; - int loops; - int which_proc; int nProcessors; std::vector atomsPerProc; int nGlobalMols = info->getNGlobalMolecules(); - std::vector molToProcMap(nGlobalMols, -1); // default to an error condition: + std::vector molToProcMap(nGlobalMols, -1); // default to an + // error + // condition: - MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); + MPI_Comm_size( MPI_COMM_WORLD, &nProcessors); + //nProcessors = MPI::COMM_WORLD.Get_size(); if (nProcessors > nGlobalMols) { sprintf(painCave.errMsg, @@ -474,17 +539,18 @@ namespace oopse { "\tthe number of molecules. This will not result in a \n" "\tusable division of atoms for force decomposition.\n" "\tEither try a smaller number of processors, or run the\n" - "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols); + "\tsingle-processor version of OpenMD.\n", nProcessors, + nGlobalMols); painCave.isFatal = 1; simError(); } - int seedValue; Globals * simParams = info->getSimParams(); - SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator + SeqRandNumGen* myRandom; //divide labor does not need Parallel + //random number generator if (simParams->haveSeed()) { - seedValue = simParams->getSeed(); + int seedValue = simParams->getSeed(); myRandom = new SeqRandNumGen(seedValue); }else { myRandom = new SeqRandNumGen(); @@ -497,43 +563,46 @@ namespace oopse { atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0); if (worldRank == 0) { - numerator = info->getNGlobalAtoms(); - denominator = nProcessors; - precast = numerator / denominator; - nTarget = (int)(precast + 0.5); + RealType numerator = info->getNGlobalAtoms(); + RealType denominator = nProcessors; + RealType precast = numerator / denominator; + int nTarget = (int)(precast + 0.5); - for(i = 0; i < nGlobalMols; i++) { - done = 0; - loops = 0; + for(int i = 0; i < nGlobalMols; i++) { + + int done = 0; + int loops = 0; while (!done) { loops++; // Pick a processor at random - which_proc = (int) (myRandom->rand() * nProcessors); + int which_proc = (int) (myRandom->rand() * nProcessors); //get the molecule stamp first int stampId = info->getMoleculeStampId(i); MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId); // How many atoms does this processor have so far? - old_atoms = atomsPerProc[which_proc]; - add_atoms = moleculeStamp->getNAtoms(); - new_atoms = old_atoms + add_atoms; + int old_atoms = atomsPerProc[which_proc]; + int add_atoms = moleculeStamp->getNAtoms(); + int new_atoms = old_atoms + add_atoms; // If we've been through this loop too many times, we need // to just give up and assign the molecule to this processor // and be done with it. if (loops > 100) { + sprintf(painCave.errMsg, - "I've tried 100 times to assign molecule %d to a " - " processor, but can't find a good spot.\n" - "I'm assigning it at random to processor %d.\n", - i, which_proc); - + "There have been 100 attempts to assign molecule %d to an\n" + "\tunderworked processor, but there's no good place to\n" + "\tleave it. OpenMD is assigning it at random to processor %d.\n", + i, which_proc); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; simError(); molToProcMap[i] = which_proc; @@ -562,8 +631,8 @@ namespace oopse { // Pacc(x) = exp(- a * x) // where a = penalty / (average atoms per molecule) - x = (RealType)(new_atoms - nTarget); - y = myRandom->rand(); + RealType x = (RealType)(new_atoms - nTarget); + RealType y = myRandom->rand(); if (y < exp(- a * x)) { molToProcMap[i] = which_proc; @@ -578,21 +647,22 @@ namespace oopse { } delete myRandom; - + // Spray out this nonsense to all other processors: - MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); } else { // Listen to your marching orders from processor 0: - MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); + } info->setMolToProcMap(molToProcMap); sprintf(checkPointMsg, "Successfully divided the molecules among the processors.\n"); - MPIcheckPoint(); + errorCheckPoint(); } #endif @@ -609,8 +679,10 @@ namespace oopse { #endif stampId = info->getMoleculeStampId(i); - Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId), - stampId, i, info->getLocalIndexManager()); + Molecule * mol = molCreator.createMolecule(info->getForceField(), + info->getMoleculeStamp(stampId), + stampId, i, + info->getLocalIndexManager()); info->addMolecule(mol); @@ -622,132 +694,276 @@ namespace oopse { } //end for(int i=0) } + + int SimCreator::computeStorageLayout(SimInfo* info) { + + Globals* simParams = info->getSimParams(); + int nRigidBodies = info->getNGlobalRigidBodies(); + set atomTypes = info->getSimulatedAtomTypes(); + set::iterator i; + bool hasDirectionalAtoms = false; + bool hasFixedCharge = false; + bool hasDipoles = false; + bool hasQuadrupoles = false; + bool hasPolarizable = false; + bool hasFluctuatingCharge = false; + bool hasMetallic = false; + int storageLayout = 0; + storageLayout |= DataStorage::dslPosition; + storageLayout |= DataStorage::dslVelocity; + storageLayout |= DataStorage::dslForce; + + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + + DirectionalAdapter da = DirectionalAdapter( (*i) ); + MultipoleAdapter ma = MultipoleAdapter( (*i) ); + EAMAdapter ea = EAMAdapter( (*i) ); + SuttonChenAdapter sca = SuttonChenAdapter( (*i) ); + PolarizableAdapter pa = PolarizableAdapter( (*i) ); + FixedChargeAdapter fca = FixedChargeAdapter( (*i) ); + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) ); + + if (da.isDirectional()){ + hasDirectionalAtoms = true; + } + if (ma.isDipole()){ + hasDipoles = true; + } + if (ma.isQuadrupole()){ + hasQuadrupoles = true; + } + if (ea.isEAM() || sca.isSuttonChen()){ + hasMetallic = true; + } + if ( fca.isFixedCharge() ){ + hasFixedCharge = true; + } + if ( fqa.isFluctuatingCharge() ){ + hasFluctuatingCharge = true; + } + if ( pa.isPolarizable() ){ + hasPolarizable = true; + } + } + if (nRigidBodies > 0 || hasDirectionalAtoms) { + storageLayout |= DataStorage::dslAmat; + if(storageLayout & DataStorage::dslVelocity) { + storageLayout |= DataStorage::dslAngularMomentum; + } + if (storageLayout & DataStorage::dslForce) { + storageLayout |= DataStorage::dslTorque; + } + } + if (hasDipoles) { + storageLayout |= DataStorage::dslDipole; + } + if (hasQuadrupoles) { + storageLayout |= DataStorage::dslQuadrupole; + } + if (hasFixedCharge || hasFluctuatingCharge) { + storageLayout |= DataStorage::dslSkippedCharge; + } + if (hasMetallic) { + storageLayout |= DataStorage::dslDensity; + storageLayout |= DataStorage::dslFunctional; + storageLayout |= DataStorage::dslFunctionalDerivative; + } + if (hasPolarizable) { + storageLayout |= DataStorage::dslElectricField; + } + if (hasFluctuatingCharge){ + storageLayout |= DataStorage::dslFlucQPosition; + if(storageLayout & DataStorage::dslVelocity) { + storageLayout |= DataStorage::dslFlucQVelocity; + } + if (storageLayout & DataStorage::dslForce) { + storageLayout |= DataStorage::dslFlucQForce; + } + } + + // if the user has asked for them, make sure we've got the memory for the + // objects defined. + + if (simParams->getOutputParticlePotential()) { + storageLayout |= DataStorage::dslParticlePot; + } + + if (simParams->havePrintHeatFlux()) { + if (simParams->getPrintHeatFlux()) { + storageLayout |= DataStorage::dslParticlePot; + } + } + + if (simParams->getOutputElectricField() | simParams->haveElectricField()) { + storageLayout |= DataStorage::dslElectricField; + } + + if (simParams->getOutputFluctuatingCharges()) { + storageLayout |= DataStorage::dslFlucQPosition; + storageLayout |= DataStorage::dslFlucQVelocity; + storageLayout |= DataStorage::dslFlucQForce; + } + + info->setStorageLayout(storageLayout); + + return storageLayout; + } + void SimCreator::setGlobalIndex(SimInfo *info) { SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; Molecule::RigidBodyIterator ri; Molecule::CutoffGroupIterator ci; + Molecule::BondIterator boi; + Molecule::BendIterator bei; + Molecule::TorsionIterator ti; + Molecule::InversionIterator ii; Molecule::IntegrableObjectIterator ioi; - Molecule * mol; - Atom * atom; - RigidBody * rb; - CutoffGroup * cg; + Molecule* mol; + Atom* atom; + RigidBody* rb; + CutoffGroup* cg; + Bond* bond; + Bend* bend; + Torsion* torsion; + Inversion* inversion; int beginAtomIndex; int beginRigidBodyIndex; int beginCutoffGroupIndex; + int beginBondIndex; + int beginBendIndex; + int beginTorsionIndex; + int beginInversionIndex; int nGlobalAtoms = info->getNGlobalAtoms(); - - /**@todo fixme */ -#ifndef IS_MPI + int nGlobalRigidBodies = info->getNGlobalRigidBodies(); beginAtomIndex = 0; - beginRigidBodyIndex = 0; + // The rigid body indices begin immediately after the atom indices: + beginRigidBodyIndex = info->getNGlobalAtoms(); beginCutoffGroupIndex = 0; - -#else - - int nproc; - int myNode; - - myNode = worldRank; - MPI_Comm_size(MPI_COMM_WORLD, &nproc); - - std::vector < int > tmpAtomsInProc(nproc, 0); - std::vector < int > tmpRigidBodiesInProc(nproc, 0); - std::vector < int > tmpCutoffGroupsInProc(nproc, 0); - std::vector < int > NumAtomsInProc(nproc, 0); - std::vector < int > NumRigidBodiesInProc(nproc, 0); - std::vector < int > NumCutoffGroupsInProc(nproc, 0); - - tmpAtomsInProc[myNode] = info->getNAtoms(); - tmpRigidBodiesInProc[myNode] = info->getNRigidBodies(); - tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups(); - - //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups - MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT, - MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, - MPI_INT, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, - MPI_INT, MPI_SUM, MPI_COMM_WORLD); - - beginAtomIndex = 0; - beginRigidBodyIndex = 0; - beginCutoffGroupIndex = 0; - - for(int i = 0; i < myNode; i++) { - beginAtomIndex += NumAtomsInProc[i]; - beginRigidBodyIndex += NumRigidBodiesInProc[i]; - beginCutoffGroupIndex += NumCutoffGroupsInProc[i]; - } - -#endif - - //rigidbody's index begins right after atom's - beginRigidBodyIndex += info->getNGlobalAtoms(); - - for(mol = info->beginMolecule(mi); mol != NULL; - mol = info->nextMolecule(mi)) { + beginBondIndex = 0; + beginBendIndex = 0; + beginTorsionIndex = 0; + beginInversionIndex = 0; + + for(int i = 0; i < info->getNGlobalMolecules(); i++) { - //local index(index in DataStorge) of atom is important - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - atom->setGlobalIndex(beginAtomIndex++); +#ifdef IS_MPI + if (info->getMolToProc(i) == worldRank) { +#endif + // stuff to do if I own this molecule + mol = info->getMoleculeByGlobalIndex(i); + + // The local index(index in DataStorge) of the atom is important: + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + atom->setGlobalIndex(beginAtomIndex++); + } + + for(rb = mol->beginRigidBody(ri); rb != NULL; + rb = mol->nextRigidBody(ri)) { + rb->setGlobalIndex(beginRigidBodyIndex++); + } + + // The local index of other objects only depends on the order + // of traversal: + for(cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { + cg->setGlobalIndex(beginCutoffGroupIndex++); + } + for(bond = mol->beginBond(boi); bond != NULL; + bond = mol->nextBond(boi)) { + bond->setGlobalIndex(beginBondIndex++); + } + for(bend = mol->beginBend(bei); bend != NULL; + bend = mol->nextBend(bei)) { + bend->setGlobalIndex(beginBendIndex++); + } + for(torsion = mol->beginTorsion(ti); torsion != NULL; + torsion = mol->nextTorsion(ti)) { + torsion->setGlobalIndex(beginTorsionIndex++); + } + for(inversion = mol->beginInversion(ii); inversion != NULL; + inversion = mol->nextInversion(ii)) { + inversion->setGlobalIndex(beginInversionIndex++); + } + +#ifdef IS_MPI + } else { + + // stuff to do if I don't own this molecule + + int stampId = info->getMoleculeStampId(i); + MoleculeStamp* stamp = info->getMoleculeStamp(stampId); + + beginAtomIndex += stamp->getNAtoms(); + beginRigidBodyIndex += stamp->getNRigidBodies(); + beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); + beginBondIndex += stamp->getNBonds(); + beginBendIndex += stamp->getNBends(); + beginTorsionIndex += stamp->getNTorsions(); + beginInversionIndex += stamp->getNInversions(); } - - for(rb = mol->beginRigidBody(ri); rb != NULL; - rb = mol->nextRigidBody(ri)) { - rb->setGlobalIndex(beginRigidBodyIndex++); - } - - //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)) { - cg->setGlobalIndex(beginCutoffGroupIndex++); - } - } - +#endif + + } //end for(int i=0) + //fill globalGroupMembership std::vector globalGroupMembership(info->getNGlobalAtoms(), 0); - for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { - for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { - + for(mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex(); } } } - + #ifdef IS_MPI // Since the globalGroupMembership has been zero filled and we've only // poked values into the atoms we know, we can do an Allreduce // to get the full globalGroupMembership array (We think). // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 // docs said we could. - std::vector tmpGroupMembership(nGlobalAtoms, 0); - MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, + std::vector tmpGroupMembership(info->getNGlobalAtoms(), 0); + MPI_Allreduce(&globalGroupMembership[0], + &tmpGroupMembership[0], nGlobalAtoms, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0], + // &tmpGroupMembership[0], nGlobalAtoms, + // MPI::INT, MPI::SUM); info->setGlobalGroupMembership(tmpGroupMembership); #else info->setGlobalGroupMembership(globalGroupMembership); #endif //fill molMembership - std::vector globalMolMembership(info->getNGlobalAtoms(), 0); + std::vector globalMolMembership(info->getNGlobalAtoms() + + info->getNGlobalRigidBodies(), 0); - for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { - + for(mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); } + for (rb = mol->beginRigidBody(ri); rb != NULL; + rb = mol->nextRigidBody(ri)) { + globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex(); + } } #ifdef IS_MPI - std::vector tmpMolMembership(nGlobalAtoms, 0); - - MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, + std::vector tmpMolMembership(info->getNGlobalAtoms() + + info->getNGlobalRigidBodies(), 0); + MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], + nGlobalAtoms + nGlobalRigidBodies, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0], + // nGlobalAtoms + nGlobalRigidBodies, + // MPI::INT, MPI::SUM); info->setGlobalMolMembership(tmpMolMembership); #else @@ -758,47 +974,48 @@ namespace oopse { // here the molecules are listed by their global indices. std::vector nIOPerMol(info->getNGlobalMolecules(), 0); - for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { + for (mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects(); } #ifdef IS_MPI std::vector numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0); MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], - info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); + info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], + // info->getNGlobalMolecules(), MPI::INT, MPI::SUM); #else std::vector numIntegrableObjectsPerMol = nIOPerMol; #endif - std::vector startingIOIndexForMol(info->getNGlobalMolecules()); - -int startingIndex = 0; - for (int i = 0; i < info->getNGlobalMolecules(); i++) { - startingIOIndexForMol[i] = startingIndex; - startingIndex += numIntegrableObjectsPerMol[i]; - } - - std::vector IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); - for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { + std::vector startingIOIndexForMol(info->getNGlobalMolecules()); + + int startingIndex = 0; + for (int i = 0; i < info->getNGlobalMolecules(); i++) { + startingIOIndexForMol[i] = startingIndex; + startingIndex += numIntegrableObjectsPerMol[i]; + } + + std::vector IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); + for (mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { int myGlobalIndex = mol->getGlobalIndex(); int globalIO = startingIOIndexForMol[myGlobalIndex]; - for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ioi)) { - integrableObject->setGlobalIntegrableObjectIndex(globalIO); - IOIndexToIntegrableObject[globalIO] = integrableObject; - globalIO++; + for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; + sd = mol->nextIntegrableObject(ioi)) { + sd->setGlobalIntegrableObjectIndex(globalIO); + IOIndexToIntegrableObject[globalIO] = sd; + globalIO++; } } - - info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); - + + info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); + } void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { - Globals* simParams; - simParams = info->getSimParams(); - DumpReader reader(info, mdFileName); int nframes = reader.getNFrames(); @@ -812,11 +1029,10 @@ int startingIndex = 0; painCave.isFatal = 1; simError(); } - //copy the current snapshot to previous snapshot info->getSnapshotManager()->advance(); } -} //end namespace oopse +} //end namespace OpenMD