| 6 |  | * redistribute this software in source and binary code form, provided | 
| 7 |  | * that the following conditions are met: | 
| 8 |  | * | 
| 9 | < | * 1. Acknowledgement of the program authors must be made in any | 
| 10 | < | *    publication of scientific results based in part on use of the | 
| 11 | < | *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | < | *    the article in which the program was described (Matthew | 
| 13 | < | *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | < | *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | < | *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | < | *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | < | * | 
| 18 | < | * 2. Redistributions of source code must retain the above copyright | 
| 9 | > | * 1. Redistributions of source code must retain the above copyright | 
| 10 |  | *    notice, this list of conditions and the following disclaimer. | 
| 11 |  | * | 
| 12 | < | * 3. Redistributions in binary form must reproduce the above copyright | 
| 12 | > | * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 |  | *    notice, this list of conditions and the following disclaimer in the | 
| 14 |  | *    documentation and/or other materials provided with the | 
| 15 |  | *    distribution. | 
| 28 |  | * arising out of the use of or inability to use software, even if the | 
| 29 |  | * University of Notre Dame has been advised of the possibility of | 
| 30 |  | * such damages. | 
| 31 | + | * | 
| 32 | + | * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your | 
| 33 | + | * research, please cite the appropriate papers when you publish your | 
| 34 | + | * work.  Good starting points are: | 
| 35 | + | * | 
| 36 | + | * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | 
| 37 | + | * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | 
| 38 | + | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | 
| 39 | + | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | + | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 |  | */ | 
| 42 |  |  | 
| 43 |  | /** | 
| 56 |  | #include "brains/SimCreator.hpp" | 
| 57 |  | #include "brains/SimSnapshotManager.hpp" | 
| 58 |  | #include "io/DumpReader.hpp" | 
| 59 | < | #include "UseTheForce/ForceFieldFactory.hpp" | 
| 59 | > | #include "brains/ForceField.hpp" | 
| 60 |  | #include "utils/simError.h" | 
| 61 |  | #include "utils/StringUtils.hpp" | 
| 62 |  | #include "math/SeqRandNumGen.hpp" | 
| 76 |  | #include "antlr/NoViableAltForCharException.hpp" | 
| 77 |  | #include "antlr/NoViableAltException.hpp" | 
| 78 |  |  | 
| 79 | + | #include "types/DirectionalAdapter.hpp" | 
| 80 | + | #include "types/MultipoleAdapter.hpp" | 
| 81 | + | #include "types/EAMAdapter.hpp" | 
| 82 | + | #include "types/SuttonChenAdapter.hpp" | 
| 83 | + | #include "types/PolarizableAdapter.hpp" | 
| 84 | + | #include "types/FixedChargeAdapter.hpp" | 
| 85 | + | #include "types/FluctuatingChargeAdapter.hpp" | 
| 86 | + |  | 
| 87 |  | #ifdef IS_MPI | 
| 88 | + | #include "mpi.h" | 
| 89 |  | #include "math/ParallelRandNumGen.hpp" | 
| 90 |  | #endif | 
| 91 |  |  | 
| 92 | < | namespace oopse { | 
| 92 | > | namespace OpenMD { | 
| 93 |  |  | 
| 94 | < | Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){ | 
| 94 | > | Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){ | 
| 95 |  | Globals* simParams = NULL; | 
| 96 |  | try { | 
| 97 |  |  | 
| 102 |  | const int masterNode = 0; | 
| 103 |  | int commStatus; | 
| 104 |  | if (worldRank == masterNode) { | 
| 105 | < | #endif | 
| 106 | < |  | 
| 105 | > | commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 106 | > | #endif | 
| 107 |  | SimplePreprocessor preprocessor; | 
| 108 |  | preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); | 
| 109 |  |  | 
| 116 |  |  | 
| 117 |  |  | 
| 118 |  | } else { | 
| 119 | + |  | 
| 120 | + | commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 121 | + |  | 
| 122 |  | //get stream size | 
| 123 |  | commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); | 
| 124 | < |  | 
| 124 | > |  | 
| 125 |  | char* buf = new char[streamSize]; | 
| 126 |  | assert(buf); | 
| 127 |  |  | 
| 129 |  | commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 130 |  |  | 
| 131 |  | ppStream.str(buf); | 
| 132 | < | delete buf; | 
| 132 | > | delete [] buf; | 
| 133 |  |  | 
| 134 |  | } | 
| 135 |  | #endif | 
| 158 |  | treeParser.initializeASTFactory(factory); | 
| 159 |  | treeParser.setASTFactory(&factory); | 
| 160 |  | simParams = treeParser.walkTree(parser.getAST()); | 
| 148 | – |  | 
| 161 |  | } | 
| 162 |  |  | 
| 163 |  |  | 
| 227 |  | painCave.isFatal = 1; | 
| 228 |  | simError(); | 
| 229 |  | } | 
| 230 | < | catch (OOPSEException& e) { | 
| 230 | > | catch (OpenMDException& e) { | 
| 231 |  | sprintf(painCave.errMsg, | 
| 232 |  | "%s\n", | 
| 233 |  | e.getMessage().c_str()); | 
| 242 |  | simError(); | 
| 243 |  | } | 
| 244 |  |  | 
| 245 | + | simParams->setMDfileVersion(mdFileVersion); | 
| 246 |  | return simParams; | 
| 247 |  | } | 
| 248 |  |  | 
| 249 |  | SimInfo*  SimCreator::createSim(const std::string & mdFileName, | 
| 250 |  | bool loadInitCoords) { | 
| 251 | < |  | 
| 251 | > |  | 
| 252 |  | const int bufferSize = 65535; | 
| 253 |  | char buffer[bufferSize]; | 
| 254 |  | int lineNo = 0; | 
| 256 |  | int metaDataBlockStart = -1; | 
| 257 |  | int metaDataBlockEnd = -1; | 
| 258 |  | int i; | 
| 259 | < | int mdOffset; | 
| 259 | > | streamoff mdOffset(0); | 
| 260 | > | int mdFileVersion; | 
| 261 |  |  | 
| 262 | + |  | 
| 263 |  | #ifdef IS_MPI | 
| 264 |  | const int masterNode = 0; | 
| 265 |  | if (worldRank == masterNode) { | 
| 266 |  | #endif | 
| 267 |  |  | 
| 268 | < | std::ifstream mdFile_(mdFileName.c_str()); | 
| 268 | > | std::ifstream mdFile_; | 
| 269 | > | mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary); | 
| 270 |  |  | 
| 271 |  | if (mdFile_.fail()) { | 
| 272 |  | sprintf(painCave.errMsg, | 
| 279 |  | mdFile_.getline(buffer, bufferSize); | 
| 280 |  | ++lineNo; | 
| 281 |  | std::string line = trimLeftCopy(buffer); | 
| 282 | < | i = CaseInsensitiveFind(line, "<OOPSE"); | 
| 283 | < | if (i == string::npos) { | 
| 282 | > | i = CaseInsensitiveFind(line, "<OpenMD"); | 
| 283 | > | if (static_cast<size_t>(i) == string::npos) { | 
| 284 | > | // try the older file strings to see if that works: | 
| 285 | > | i = CaseInsensitiveFind(line, "<OOPSE"); | 
| 286 | > | } | 
| 287 | > |  | 
| 288 | > | if (static_cast<size_t>(i) == string::npos) { | 
| 289 | > | // still no luck! | 
| 290 |  | sprintf(painCave.errMsg, | 
| 291 | < | "SimCreator: File: %s is not an OOPSE file!\n", | 
| 291 | > | "SimCreator: File: %s is not a valid OpenMD file!\n", | 
| 292 |  | mdFileName.c_str()); | 
| 293 |  | painCave.isFatal = 1; | 
| 294 |  | simError(); | 
| 295 |  | } | 
| 296 | + |  | 
| 297 | + | // found the correct opening string, now try to get the file | 
| 298 | + | // format version number. | 
| 299 |  |  | 
| 300 | + | StringTokenizer tokenizer(line, "=<> \t\n\r"); | 
| 301 | + | std::string fileType = tokenizer.nextToken(); | 
| 302 | + | toUpper(fileType); | 
| 303 | + |  | 
| 304 | + | mdFileVersion = 0; | 
| 305 | + |  | 
| 306 | + | if (fileType == "OPENMD") { | 
| 307 | + | while (tokenizer.hasMoreTokens()) { | 
| 308 | + | std::string token(tokenizer.nextToken()); | 
| 309 | + | toUpper(token); | 
| 310 | + | if (token == "VERSION") { | 
| 311 | + | mdFileVersion = tokenizer.nextTokenAsInt(); | 
| 312 | + | break; | 
| 313 | + | } | 
| 314 | + | } | 
| 315 | + | } | 
| 316 | + |  | 
| 317 |  | //scan through the input stream and find MetaData tag | 
| 318 |  | while(mdFile_.getline(buffer, bufferSize)) { | 
| 319 |  | ++lineNo; | 
| 369 |  | std::stringstream rawMetaDataStream(mdRawData); | 
| 370 |  |  | 
| 371 |  | //parse meta-data file | 
| 372 | < | Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1); | 
| 372 | > | Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion, | 
| 373 | > | metaDataBlockStart + 1); | 
| 374 |  |  | 
| 375 |  | //create the force field | 
| 376 | < | ForceField * ff = ForceFieldFactory::getInstance() | 
| 377 | < | ->createForceField(simParams->getForceField()); | 
| 335 | < |  | 
| 376 | > | ForceField * ff = new ForceField(simParams->getForceField()); | 
| 377 | > |  | 
| 378 |  | if (ff == NULL) { | 
| 379 |  | sprintf(painCave.errMsg, | 
| 380 |  | "ForceField Factory can not create %s force field\n", | 
| 407 |  | } | 
| 408 |  |  | 
| 409 |  | ff->parse(forcefieldFileName); | 
| 368 | – | ff->setFortranForceOptions(); | 
| 410 |  | //create SimInfo | 
| 411 |  | SimInfo * info = new SimInfo(ff, simParams); | 
| 412 |  |  | 
| 424 |  | //create the molecules | 
| 425 |  | createMolecules(info); | 
| 426 |  |  | 
| 427 | < |  | 
| 427 | > | //find the storage layout | 
| 428 | > |  | 
| 429 | > | int storageLayout = computeStorageLayout(info); | 
| 430 | > |  | 
| 431 |  | //allocate memory for DataStorage(circular reference, need to | 
| 432 |  | //break it) | 
| 433 | < | info->setSnapshotManager(new SimSnapshotManager(info)); | 
| 433 | > | info->setSnapshotManager(new SimSnapshotManager(info, storageLayout)); | 
| 434 |  |  | 
| 435 |  | //set the global index of atoms, rigidbodies and cutoffgroups | 
| 436 |  | //(only need to be set once, the global index will never change | 
| 439 |  | //responsibility to LocalIndexManager. | 
| 440 |  | setGlobalIndex(info); | 
| 441 |  |  | 
| 442 | < | //Although addExcludePairs is called inside SimInfo's addMolecule | 
| 442 | > | //Although addInteractionPairs is called inside SimInfo's addMolecule | 
| 443 |  | //method, at that point atoms don't have the global index yet | 
| 444 |  | //(their global index are all initialized to -1).  Therefore we | 
| 445 | < | //have to call addExcludePairs explicitly here. A way to work | 
| 445 | > | //have to call addInteractionPairs explicitly here. A way to work | 
| 446 |  | //around is that we can determine the beginning global indices of | 
| 447 |  | //atoms before they get created. | 
| 448 |  | SimInfo::MoleculeIterator mi; | 
| 449 |  | Molecule* mol; | 
| 450 |  | for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 451 | < | info->addExcludePairs(mol); | 
| 451 | > | info->addInteractionPairs(mol); | 
| 452 |  | } | 
| 453 |  |  | 
| 454 |  | if (loadInitCoords) | 
| 455 |  | loadCoordinates(info, mdFileName); | 
| 412 | – |  | 
| 456 |  | return info; | 
| 457 |  | } | 
| 458 |  |  | 
| 516 |  | "\tthe number of molecules.  This will not result in a \n" | 
| 517 |  | "\tusable division of atoms for force decomposition.\n" | 
| 518 |  | "\tEither try a smaller number of processors, or run the\n" | 
| 519 | < | "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols); | 
| 519 | > | "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols); | 
| 520 |  |  | 
| 521 |  | painCave.isFatal = 1; | 
| 522 |  | simError(); | 
| 651 |  | #endif | 
| 652 |  |  | 
| 653 |  | stampId = info->getMoleculeStampId(i); | 
| 654 | < | Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId), | 
| 655 | < | stampId, i, info->getLocalIndexManager()); | 
| 654 | > | Molecule * mol = molCreator.createMolecule(info->getForceField(), | 
| 655 | > | info->getMoleculeStamp(stampId), | 
| 656 | > | stampId, i, | 
| 657 | > | info->getLocalIndexManager()); | 
| 658 |  |  | 
| 659 |  | info->addMolecule(mol); | 
| 660 |  |  | 
| 667 |  | } //end for(int i=0) | 
| 668 |  | } | 
| 669 |  |  | 
| 670 | + | int SimCreator::computeStorageLayout(SimInfo* info) { | 
| 671 | + |  | 
| 672 | + | Globals* simParams = info->getSimParams(); | 
| 673 | + | int nRigidBodies = info->getNGlobalRigidBodies(); | 
| 674 | + | set<AtomType*> atomTypes = info->getSimulatedAtomTypes(); | 
| 675 | + | set<AtomType*>::iterator i; | 
| 676 | + | bool hasDirectionalAtoms = false; | 
| 677 | + | bool hasFixedCharge = false; | 
| 678 | + | bool hasMultipoles = false; | 
| 679 | + | bool hasPolarizable = false; | 
| 680 | + | bool hasFluctuatingCharge = false; | 
| 681 | + | bool hasMetallic = false; | 
| 682 | + | int storageLayout = 0; | 
| 683 | + | storageLayout |= DataStorage::dslPosition; | 
| 684 | + | storageLayout |= DataStorage::dslVelocity; | 
| 685 | + | storageLayout |= DataStorage::dslForce; | 
| 686 | + |  | 
| 687 | + | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 688 | + |  | 
| 689 | + | DirectionalAdapter da = DirectionalAdapter( (*i) ); | 
| 690 | + | MultipoleAdapter ma = MultipoleAdapter( (*i) ); | 
| 691 | + | EAMAdapter ea = EAMAdapter( (*i) ); | 
| 692 | + | SuttonChenAdapter sca = SuttonChenAdapter( (*i) ); | 
| 693 | + | PolarizableAdapter pa = PolarizableAdapter( (*i) ); | 
| 694 | + | FixedChargeAdapter fca = FixedChargeAdapter( (*i) ); | 
| 695 | + | FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) ); | 
| 696 | + |  | 
| 697 | + | if (da.isDirectional()){ | 
| 698 | + | hasDirectionalAtoms = true; | 
| 699 | + | } | 
| 700 | + | if (ma.isMultipole()){ | 
| 701 | + | hasMultipoles = true; | 
| 702 | + | } | 
| 703 | + | if (ea.isEAM() || sca.isSuttonChen()){ | 
| 704 | + | hasMetallic = true; | 
| 705 | + | } | 
| 706 | + | if ( fca.isFixedCharge() ){ | 
| 707 | + | hasFixedCharge = true; | 
| 708 | + | } | 
| 709 | + | if ( fqa.isFluctuatingCharge() ){ | 
| 710 | + | hasFluctuatingCharge = true; | 
| 711 | + | } | 
| 712 | + | if ( pa.isPolarizable() ){ | 
| 713 | + | hasPolarizable = true; | 
| 714 | + | } | 
| 715 | + | } | 
| 716 | + |  | 
| 717 | + | if (nRigidBodies > 0 || hasDirectionalAtoms) { | 
| 718 | + | storageLayout |= DataStorage::dslAmat; | 
| 719 | + | if(storageLayout & DataStorage::dslVelocity) { | 
| 720 | + | storageLayout |= DataStorage::dslAngularMomentum; | 
| 721 | + | } | 
| 722 | + | if (storageLayout & DataStorage::dslForce) { | 
| 723 | + | storageLayout |= DataStorage::dslTorque; | 
| 724 | + | } | 
| 725 | + | } | 
| 726 | + | if (hasMultipoles) { | 
| 727 | + | storageLayout |= DataStorage::dslElectroFrame; | 
| 728 | + | } | 
| 729 | + | if (hasFixedCharge || hasFluctuatingCharge) { | 
| 730 | + | storageLayout |= DataStorage::dslSkippedCharge; | 
| 731 | + | } | 
| 732 | + | if (hasMetallic) { | 
| 733 | + | storageLayout |= DataStorage::dslDensity; | 
| 734 | + | storageLayout |= DataStorage::dslFunctional; | 
| 735 | + | storageLayout |= DataStorage::dslFunctionalDerivative; | 
| 736 | + | } | 
| 737 | + | if (hasPolarizable) { | 
| 738 | + | storageLayout |= DataStorage::dslElectricField; | 
| 739 | + | } | 
| 740 | + | if (hasFluctuatingCharge){ | 
| 741 | + | storageLayout |= DataStorage::dslFlucQPosition; | 
| 742 | + | if(storageLayout & DataStorage::dslVelocity) { | 
| 743 | + | storageLayout |= DataStorage::dslFlucQVelocity; | 
| 744 | + | } | 
| 745 | + | if (storageLayout & DataStorage::dslForce) { | 
| 746 | + | storageLayout |= DataStorage::dslFlucQForce; | 
| 747 | + | } | 
| 748 | + | } | 
| 749 | + |  | 
| 750 | + | // if the user has asked for them, make sure we've got the memory for the | 
| 751 | + | // objects defined. | 
| 752 | + |  | 
| 753 | + | if (simParams->getOutputParticlePotential()) { | 
| 754 | + | storageLayout |= DataStorage::dslParticlePot; | 
| 755 | + | } | 
| 756 | + |  | 
| 757 | + | if (simParams->havePrintHeatFlux()) { | 
| 758 | + | if (simParams->getPrintHeatFlux()) { | 
| 759 | + | storageLayout |= DataStorage::dslParticlePot; | 
| 760 | + | } | 
| 761 | + | } | 
| 762 | + |  | 
| 763 | + | if (simParams->getOutputElectricField()) { | 
| 764 | + | storageLayout |= DataStorage::dslElectricField; | 
| 765 | + | } | 
| 766 | + | if (simParams->getOutputFluctuatingCharges()) { | 
| 767 | + | storageLayout |= DataStorage::dslFlucQPosition; | 
| 768 | + | storageLayout |= DataStorage::dslFlucQVelocity; | 
| 769 | + | storageLayout |= DataStorage::dslFlucQForce; | 
| 770 | + | } | 
| 771 | + |  | 
| 772 | + | return storageLayout; | 
| 773 | + | } | 
| 774 | + |  | 
| 775 |  | void SimCreator::setGlobalIndex(SimInfo *info) { | 
| 776 |  | SimInfo::MoleculeIterator mi; | 
| 777 |  | Molecule::AtomIterator ai; | 
| 786 |  | int beginRigidBodyIndex; | 
| 787 |  | int beginCutoffGroupIndex; | 
| 788 |  | int nGlobalAtoms = info->getNGlobalAtoms(); | 
| 639 | – |  | 
| 640 | – | /**@todo fixme */ | 
| 641 | – | #ifndef IS_MPI | 
| 789 |  |  | 
| 790 |  | beginAtomIndex = 0; | 
| 644 | – | beginRigidBodyIndex = 0; | 
| 645 | – | beginCutoffGroupIndex = 0; | 
| 646 | – |  | 
| 647 | – | #else | 
| 648 | – |  | 
| 649 | – | int nproc; | 
| 650 | – | int myNode; | 
| 651 | – |  | 
| 652 | – | myNode = worldRank; | 
| 653 | – | MPI_Comm_size(MPI_COMM_WORLD, &nproc); | 
| 654 | – |  | 
| 655 | – | std::vector < int > tmpAtomsInProc(nproc, 0); | 
| 656 | – | std::vector < int > tmpRigidBodiesInProc(nproc, 0); | 
| 657 | – | std::vector < int > tmpCutoffGroupsInProc(nproc, 0); | 
| 658 | – | std::vector < int > NumAtomsInProc(nproc, 0); | 
| 659 | – | std::vector < int > NumRigidBodiesInProc(nproc, 0); | 
| 660 | – | std::vector < int > NumCutoffGroupsInProc(nproc, 0); | 
| 661 | – |  | 
| 662 | – | tmpAtomsInProc[myNode] = info->getNAtoms(); | 
| 663 | – | tmpRigidBodiesInProc[myNode] = info->getNRigidBodies(); | 
| 664 | – | tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups(); | 
| 665 | – |  | 
| 666 | – | //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups | 
| 667 | – | MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT, | 
| 668 | – | MPI_SUM, MPI_COMM_WORLD); | 
| 669 | – | MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, | 
| 670 | – | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 671 | – | MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, | 
| 672 | – | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 673 | – |  | 
| 674 | – | beginAtomIndex = 0; | 
| 675 | – | beginRigidBodyIndex = 0; | 
| 676 | – | beginCutoffGroupIndex = 0; | 
| 677 | – |  | 
| 678 | – | for(int i = 0; i < myNode; i++) { | 
| 679 | – | beginAtomIndex += NumAtomsInProc[i]; | 
| 680 | – | beginRigidBodyIndex += NumRigidBodiesInProc[i]; | 
| 681 | – | beginCutoffGroupIndex += NumCutoffGroupsInProc[i]; | 
| 682 | – | } | 
| 683 | – |  | 
| 684 | – | #endif | 
| 685 | – |  | 
| 791 |  | //rigidbody's index begins right after atom's | 
| 792 | < | beginRigidBodyIndex += info->getNGlobalAtoms(); | 
| 793 | < |  | 
| 794 | < | for(mol = info->beginMolecule(mi); mol != NULL; | 
| 795 | < | mol = info->nextMolecule(mi)) { | 
| 792 | > | beginRigidBodyIndex = info->getNGlobalAtoms(); | 
| 793 | > | beginCutoffGroupIndex = 0; | 
| 794 | > |  | 
| 795 | > | for(int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 796 |  |  | 
| 797 | < | //local index(index in DataStorge) of atom is important | 
| 798 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 799 | < | atom->setGlobalIndex(beginAtomIndex++); | 
| 797 | > | #ifdef IS_MPI | 
| 798 | > | if (info->getMolToProc(i) == worldRank) { | 
| 799 | > | #endif | 
| 800 | > | // stuff to do if I own this molecule | 
| 801 | > | mol = info->getMoleculeByGlobalIndex(i); | 
| 802 | > |  | 
| 803 | > | //local index(index in DataStorge) of atom is important | 
| 804 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 805 | > | atom->setGlobalIndex(beginAtomIndex++); | 
| 806 | > | } | 
| 807 | > |  | 
| 808 | > | for(rb = mol->beginRigidBody(ri); rb != NULL; | 
| 809 | > | rb = mol->nextRigidBody(ri)) { | 
| 810 | > | rb->setGlobalIndex(beginRigidBodyIndex++); | 
| 811 | > | } | 
| 812 | > |  | 
| 813 | > | //local index of cutoff group is trivial, it only depends on | 
| 814 | > | //the order of travesing | 
| 815 | > | for(cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 816 | > | cg = mol->nextCutoffGroup(ci)) { | 
| 817 | > | cg->setGlobalIndex(beginCutoffGroupIndex++); | 
| 818 | > | } | 
| 819 | > |  | 
| 820 | > | #ifdef IS_MPI | 
| 821 | > | }  else { | 
| 822 | > |  | 
| 823 | > | // stuff to do if I don't own this molecule | 
| 824 | > |  | 
| 825 | > | int stampId = info->getMoleculeStampId(i); | 
| 826 | > | MoleculeStamp* stamp = info->getMoleculeStamp(stampId); | 
| 827 | > |  | 
| 828 | > | beginAtomIndex += stamp->getNAtoms(); | 
| 829 | > | beginRigidBodyIndex += stamp->getNRigidBodies(); | 
| 830 | > | beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); | 
| 831 |  | } | 
| 832 | < |  | 
| 833 | < | for(rb = mol->beginRigidBody(ri); rb != NULL; | 
| 834 | < | rb = mol->nextRigidBody(ri)) { | 
| 835 | < | rb->setGlobalIndex(beginRigidBodyIndex++); | 
| 700 | < | } | 
| 701 | < |  | 
| 702 | < | //local index of cutoff group is trivial, it only depends on the order of travesing | 
| 703 | < | for(cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 704 | < | cg = mol->nextCutoffGroup(ci)) { | 
| 705 | < | cg->setGlobalIndex(beginCutoffGroupIndex++); | 
| 706 | < | } | 
| 707 | < | } | 
| 708 | < |  | 
| 832 | > | #endif | 
| 833 | > |  | 
| 834 | > | } //end for(int i=0) | 
| 835 | > |  | 
| 836 |  | //fill globalGroupMembership | 
| 837 |  | std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); | 
| 838 |  | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 844 |  |  | 
| 845 |  | } | 
| 846 |  | } | 
| 847 | < |  | 
| 847 | > |  | 
| 848 |  | #ifdef IS_MPI | 
| 849 |  | // Since the globalGroupMembership has been zero filled and we've only | 
| 850 |  | // poked values into the atoms we know, we can do an Allreduce | 
| 851 |  | // to get the full globalGroupMembership array (We think). | 
| 852 |  | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | 
| 853 |  | // docs said we could. | 
| 854 | < | std::vector<int> tmpGroupMembership(nGlobalAtoms, 0); | 
| 854 | > | std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); | 
| 855 |  | MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, | 
| 856 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 857 |  | info->setGlobalGroupMembership(tmpGroupMembership); | 
| 863 |  | std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); | 
| 864 |  |  | 
| 865 |  | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 739 | – |  | 
| 866 |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 867 |  | globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 868 |  | } | 
| 869 |  | } | 
| 870 |  |  | 
| 871 |  | #ifdef IS_MPI | 
| 872 | < | std::vector<int> tmpMolMembership(nGlobalAtoms, 0); | 
| 872 | > | std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0); | 
| 873 |  |  | 
| 874 |  | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, | 
| 875 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 895 |  | std::vector<int> numIntegrableObjectsPerMol = nIOPerMol; | 
| 896 |  | #endif | 
| 897 |  |  | 
| 898 | < | std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules()); | 
| 899 | < |  | 
| 900 | < | int startingIndex = 0; | 
| 901 | < | for (int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 902 | < | startingIOIndexForMol[i] = startingIndex; | 
| 903 | < | startingIndex += numIntegrableObjectsPerMol[i]; | 
| 904 | < | } | 
| 905 | < |  | 
| 906 | < | std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | 
| 907 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 898 | > | std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules()); | 
| 899 | > |  | 
| 900 | > | int startingIndex = 0; | 
| 901 | > | for (int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 902 | > | startingIOIndexForMol[i] = startingIndex; | 
| 903 | > | startingIndex += numIntegrableObjectsPerMol[i]; | 
| 904 | > | } | 
| 905 | > |  | 
| 906 | > | std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | 
| 907 | > | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 908 |  | int myGlobalIndex = mol->getGlobalIndex(); | 
| 909 |  | int globalIO = startingIOIndexForMol[myGlobalIndex]; | 
| 910 | < | for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL; | 
| 911 | < | integrableObject = mol->nextIntegrableObject(ioi)) { | 
| 912 | < | integrableObject->setGlobalIntegrableObjectIndex(globalIO); | 
| 913 | < | IOIndexToIntegrableObject[globalIO] = integrableObject; | 
| 914 | < | globalIO++; | 
| 910 | > | for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; | 
| 911 | > | sd = mol->nextIntegrableObject(ioi)) { | 
| 912 | > | sd->setGlobalIntegrableObjectIndex(globalIO); | 
| 913 | > | IOIndexToIntegrableObject[globalIO] = sd; | 
| 914 | > | globalIO++; | 
| 915 |  | } | 
| 916 |  | } | 
| 917 | < |  | 
| 918 | < | info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); | 
| 919 | < |  | 
| 917 | > |  | 
| 918 | > | info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); | 
| 919 | > |  | 
| 920 |  | } | 
| 921 |  |  | 
| 922 |  | void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { | 
| 923 |  | Globals* simParams; | 
| 924 | + |  | 
| 925 |  | simParams = info->getSimParams(); | 
| 926 |  |  | 
| 800 | – |  | 
| 927 |  | DumpReader reader(info, mdFileName); | 
| 928 |  | int nframes = reader.getNFrames(); | 
| 929 | < |  | 
| 929 | > |  | 
| 930 |  | if (nframes > 0) { | 
| 931 |  | reader.readFrame(nframes - 1); | 
| 932 |  | } else { | 
| 937 |  | painCave.isFatal = 1; | 
| 938 |  | simError(); | 
| 939 |  | } | 
| 814 | – |  | 
| 940 |  | //copy the current snapshot to previous snapshot | 
| 941 |  | info->getSnapshotManager()->advance(); | 
| 942 |  | } | 
| 943 |  |  | 
| 944 | < | } //end namespace oopse | 
| 944 | > | } //end namespace OpenMD | 
| 945 |  |  | 
| 946 |  |  |