--- trunk/src/brains/SimCreator.cpp 2012/09/10 18:38:44 1796 +++ trunk/src/brains/SimCreator.cpp 2014/03/12 20:01:15 1976 @@ -1,5 +1,5 @@ /* - * copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * Copyright (c) 2005 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,7 +35,7 @@ * * [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). + * [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). */ @@ -44,9 +44,14 @@ * @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 @@ -59,6 +64,7 @@ #include "brains/ForceField.hpp" #include "utils/simError.h" #include "utils/StringUtils.hpp" +#include "utils/Revision.hpp" #include "math/SeqRandNumGen.hpp" #include "mdParser/MDLexer.hpp" #include "mdParser/MDParser.hpp" @@ -84,10 +90,6 @@ #include "types/FixedChargeAdapter.hpp" #include "types/FluctuatingChargeAdapter.hpp" -#ifdef IS_MPI -#include "mpi.h" -#include "math/ParallelRandNumGen.hpp" -#endif namespace OpenMD { @@ -102,35 +104,41 @@ namespace OpenMD { const int masterNode = 0; if (worldRank == masterNode) { - MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); + 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; - 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); - + 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); + + // 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::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); + MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); + // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); //get stream size - MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); - + 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 - MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); - + 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 @@ -152,7 +160,6 @@ namespace OpenMD { parser.initializeASTFactory(factory); parser.setASTFactory(&factory); parser.mdfile(); - // Create a tree parser that reads information into Globals MDTreeParser treeParser; treeParser.initializeASTFactory(factory); @@ -255,11 +262,28 @@ namespace OpenMD { std::string mdRawData; int metaDataBlockStart = -1; int metaDataBlockEnd = -1; - int i; - streamoff mdOffset(0); + 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(g_REVISION, strnlen(g_REVISION, 20)); + //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) { @@ -354,12 +378,23 @@ namespace OpenMD { 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 @@ -487,26 +522,16 @@ namespace OpenMD { #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 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: - nProcessors = MPI::COMM_WORLD.Get_size(); + MPI_Comm_size( MPI_COMM_WORLD, &nProcessors); + //nProcessors = MPI::COMM_WORLD.Get_size(); if (nProcessors > nGlobalMols) { sprintf(painCave.errMsg, @@ -515,17 +540,18 @@ namespace OpenMD { "\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 OpenMD.\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(); @@ -538,43 +564,46 @@ namespace OpenMD { 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", + "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; @@ -603,8 +632,8 @@ namespace OpenMD { // 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; @@ -619,13 +648,16 @@ namespace OpenMD { } delete myRandom; - + // Spray out this nonsense to all other processors: - MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); + 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::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 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); @@ -672,7 +704,8 @@ namespace OpenMD { set::iterator i; bool hasDirectionalAtoms = false; bool hasFixedCharge = false; - bool hasMultipoles = false; + bool hasDipoles = false; + bool hasQuadrupoles = false; bool hasPolarizable = false; bool hasFluctuatingCharge = false; bool hasMetallic = false; @@ -694,9 +727,12 @@ namespace OpenMD { if (da.isDirectional()){ hasDirectionalAtoms = true; } - if (ma.isMultipole()){ - hasMultipoles = true; + if (ma.isDipole()){ + hasDipoles = true; } + if (ma.isQuadrupole()){ + hasQuadrupoles = true; + } if (ea.isEAM() || sca.isSuttonChen()){ hasMetallic = true; } @@ -720,9 +756,12 @@ namespace OpenMD { storageLayout |= DataStorage::dslTorque; } } - if (hasMultipoles) { - storageLayout |= DataStorage::dslElectroFrame; + if (hasDipoles) { + storageLayout |= DataStorage::dslDipole; } + if (hasQuadrupoles) { + storageLayout |= DataStorage::dslQuadrupole; + } if (hasFixedCharge || hasFluctuatingCharge) { storageLayout |= DataStorage::dslSkippedCharge; } @@ -757,15 +796,18 @@ namespace OpenMD { } } - if (simParams->getOutputElectricField()) { + 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; } @@ -774,21 +816,38 @@ namespace OpenMD { 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(); + int nGlobalRigidBodies = info->getNGlobalRigidBodies(); beginAtomIndex = 0; - //rigidbody's index begins right after atom's + // The rigid body indices begin immediately after the atom indices: beginRigidBodyIndex = info->getNGlobalAtoms(); beginCutoffGroupIndex = 0; - + beginBondIndex = 0; + beginBendIndex = 0; + beginTorsionIndex = 0; + beginInversionIndex = 0; + for(int i = 0; i < info->getNGlobalMolecules(); i++) { #ifdef IS_MPI @@ -797,8 +856,9 @@ namespace OpenMD { // stuff to do if I own this molecule mol = info->getMoleculeByGlobalIndex(i); - //local index(index in DataStorge) of atom is important - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + // 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++); } @@ -807,12 +867,28 @@ namespace OpenMD { rb->setGlobalIndex(beginRigidBodyIndex++); } - //local index of cutoff group is trivial, it only depends on - //the order of travesing + // 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 { @@ -825,6 +901,10 @@ namespace OpenMD { beginAtomIndex += stamp->getNAtoms(); beginRigidBodyIndex += stamp->getNRigidBodies(); beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); + beginBondIndex += stamp->getNBonds(); + beginBendIndex += stamp->getNBends(); + beginTorsionIndex += stamp->getNTorsions(); + beginInversionIndex += stamp->getNInversions(); } #endif @@ -832,9 +912,10 @@ namespace OpenMD { //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(); } @@ -849,28 +930,41 @@ namespace OpenMD { // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 // docs said we could. std::vector tmpGroupMembership(info->getNGlobalAtoms(), 0); - MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0], - &tmpGroupMembership[0], nGlobalAtoms, - MPI::INT, MPI::SUM); + 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(info->getNGlobalAtoms(), 0); - MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0], - nGlobalAtoms, - MPI::INT, MPI::SUM); + 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 @@ -881,14 +975,17 @@ namespace OpenMD { // 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::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], - info->getNGlobalMolecules(), MPI::INT, MPI::SUM); + MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], + 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 @@ -902,7 +999,8 @@ namespace OpenMD { } std::vector IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); - for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { + for (mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { int myGlobalIndex = mol->getGlobalIndex(); int globalIO = startingIOIndexForMol[myGlobalIndex]; for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; @@ -918,10 +1016,10 @@ namespace OpenMD { } void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { - + DumpReader reader(info, mdFileName); int nframes = reader.getNFrames(); - + if (nframes > 0) { reader.readFrame(nframes - 1); } else {