| 1 |  | /* | 
| 2 | < | * copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 2 | > | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 |  | * | 
| 4 |  | * The University of Notre Dame grants you ("Licensee") a | 
| 5 |  | * non-exclusive, royalty free, license to use, modify and | 
| 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). | 
| 38 | > | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). | 
| 39 |  | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 |  | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 |  | */ | 
| 44 |  | * @file SimCreator.cpp | 
| 45 |  | * @author tlin | 
| 46 |  | * @date 11/03/2004 | 
| 47 | – | * @time 13:51am | 
| 47 |  | * @version 1.0 | 
| 48 |  | */ | 
| 49 | + |  | 
| 50 | + | #ifdef IS_MPI | 
| 51 | + | #include "mpi.h" | 
| 52 | + | #include "math/ParallelRandNumGen.hpp" | 
| 53 | + | #endif | 
| 54 | + |  | 
| 55 |  | #include <exception> | 
| 56 |  | #include <iostream> | 
| 57 |  | #include <sstream> | 
| 64 |  | #include "brains/ForceField.hpp" | 
| 65 |  | #include "utils/simError.h" | 
| 66 |  | #include "utils/StringUtils.hpp" | 
| 67 | + | #include "utils/Revision.hpp" | 
| 68 |  | #include "math/SeqRandNumGen.hpp" | 
| 69 |  | #include "mdParser/MDLexer.hpp" | 
| 70 |  | #include "mdParser/MDParser.hpp" | 
| 90 |  | #include "types/FixedChargeAdapter.hpp" | 
| 91 |  | #include "types/FluctuatingChargeAdapter.hpp" | 
| 92 |  |  | 
| 87 | – | #ifdef IS_MPI | 
| 88 | – | #include "mpi.h" | 
| 89 | – | #include "math/ParallelRandNumGen.hpp" | 
| 90 | – | #endif | 
| 93 |  |  | 
| 94 |  | namespace OpenMD { | 
| 95 |  |  | 
| 104 |  | const int masterNode = 0; | 
| 105 |  |  | 
| 106 |  | if (worldRank == masterNode) { | 
| 107 | < | MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 107 | > | MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 108 |  | #endif | 
| 109 |  | SimplePreprocessor preprocessor; | 
| 110 | < | preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, | 
| 111 | < | ppStream); | 
| 110 | > | preprocessor.preprocess(rawMetaDataStream, filename, | 
| 111 | > | startOfMetaDataBlock, ppStream); | 
| 112 |  |  | 
| 113 |  | #ifdef IS_MPI | 
| 114 | < | //brocasting the stream size | 
| 114 | > | //broadcasting the stream size | 
| 115 |  | streamSize = ppStream.str().size() +1; | 
| 116 | < | MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 117 | < | MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), | 
| 118 | < | streamSize, MPI::CHAR, masterNode); | 
| 117 | < |  | 
| 116 | > | MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 117 | > | MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), | 
| 118 | > | streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 119 |  | } else { | 
| 120 |  |  | 
| 121 | < | MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 121 | > | MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 122 |  |  | 
| 123 |  | //get stream size | 
| 124 | < | MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 124 | < |  | 
| 124 | > | MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 125 |  | char* buf = new char[streamSize]; | 
| 126 |  | assert(buf); | 
| 127 |  |  | 
| 128 |  | //receive file content | 
| 129 | < | MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); | 
| 130 | < |  | 
| 129 | > | MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 130 | > |  | 
| 131 |  | ppStream.str(buf); | 
| 132 |  | delete [] buf; | 
| 133 | – |  | 
| 133 |  | } | 
| 134 |  | #endif | 
| 135 |  | // Create a scanner that reads from the input stream | 
| 151 |  | parser.initializeASTFactory(factory); | 
| 152 |  | parser.setASTFactory(&factory); | 
| 153 |  | parser.mdfile(); | 
| 155 | – |  | 
| 154 |  | // Create a tree parser that reads information into Globals | 
| 155 |  | MDTreeParser treeParser; | 
| 156 |  | treeParser.initializeASTFactory(factory); | 
| 253 |  | std::string mdRawData; | 
| 254 |  | int metaDataBlockStart = -1; | 
| 255 |  | int metaDataBlockEnd = -1; | 
| 256 | < | int i; | 
| 257 | < | streamoff mdOffset(0); | 
| 256 | > | int i, j; | 
| 257 | > | streamoff mdOffset; | 
| 258 |  | int mdFileVersion; | 
| 259 |  |  | 
| 260 | + | // Create a string for embedding the version information in the MetaData | 
| 261 | + | std::string version; | 
| 262 | + | version.assign("## Last run using OpenMD Version: "); | 
| 263 | + | version.append(OPENMD_VERSION_MAJOR); | 
| 264 | + | version.append("."); | 
| 265 | + | version.append(OPENMD_VERSION_MINOR); | 
| 266 |  |  | 
| 267 | + | std::string svnrev(g_REVISION, strnlen(g_REVISION, 20)); | 
| 268 | + | //convert a macro from compiler to a string in c++ | 
| 269 | + | // STR_DEFINE(svnrev, SVN_REV ); | 
| 270 | + | version.append(" Revision: "); | 
| 271 | + | // If there's no SVN revision, just call this the RELEASE revision. | 
| 272 | + | if (!svnrev.empty()) { | 
| 273 | + | version.append(svnrev); | 
| 274 | + | } else { | 
| 275 | + | version.append("RELEASE"); | 
| 276 | + | } | 
| 277 | + |  | 
| 278 |  | #ifdef IS_MPI | 
| 279 |  | const int masterNode = 0; | 
| 280 |  | if (worldRank == masterNode) { | 
| 369 |  |  | 
| 370 |  | mdRawData.clear(); | 
| 371 |  |  | 
| 372 | + | bool foundVersion = false; | 
| 373 | + |  | 
| 374 |  | for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) { | 
| 375 |  | mdFile_.getline(buffer, bufferSize); | 
| 376 | < | mdRawData += buffer; | 
| 376 | > | std::string line = trimLeftCopy(buffer); | 
| 377 | > | j = CaseInsensitiveFind(line, "## Last run using OpenMD Version"); | 
| 378 | > | if (static_cast<size_t>(j) != string::npos) { | 
| 379 | > | foundVersion = true; | 
| 380 | > | mdRawData += version; | 
| 381 | > | } else { | 
| 382 | > | mdRawData += buffer; | 
| 383 | > | } | 
| 384 |  | mdRawData += "\n"; | 
| 385 |  | } | 
| 386 | < |  | 
| 386 | > |  | 
| 387 | > | if (!foundVersion) mdRawData += version + "\n"; | 
| 388 | > |  | 
| 389 |  | mdFile_.close(); | 
| 390 |  |  | 
| 391 |  | #ifdef IS_MPI | 
| 513 |  |  | 
| 514 |  | #ifdef IS_MPI | 
| 515 |  | void SimCreator::divideMolecules(SimInfo *info) { | 
| 490 | – | RealType numerator; | 
| 491 | – | RealType denominator; | 
| 492 | – | RealType precast; | 
| 493 | – | RealType x; | 
| 494 | – | RealType y; | 
| 516 |  | RealType a; | 
| 496 | – | int old_atoms; | 
| 497 | – | int add_atoms; | 
| 498 | – | int new_atoms; | 
| 499 | – | int nTarget; | 
| 500 | – | int done; | 
| 501 | – | int i; | 
| 502 | – | int loops; | 
| 503 | – | int which_proc; | 
| 517 |  | int nProcessors; | 
| 518 |  | std::vector<int> atomsPerProc; | 
| 519 |  | int nGlobalMols = info->getNGlobalMolecules(); | 
| 520 | < | std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition: | 
| 520 | > | std::vector<int> molToProcMap(nGlobalMols, -1); // default to an | 
| 521 | > | // error | 
| 522 | > | // condition: | 
| 523 |  |  | 
| 524 | < | MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); | 
| 524 | > | MPI_Comm_size( MPI_COMM_WORLD, &nProcessors); | 
| 525 |  |  | 
| 526 |  | if (nProcessors > nGlobalMols) { | 
| 527 |  | sprintf(painCave.errMsg, | 
| 530 |  | "\tthe number of molecules.  This will not result in a \n" | 
| 531 |  | "\tusable division of atoms for force decomposition.\n" | 
| 532 |  | "\tEither try a smaller number of processors, or run the\n" | 
| 533 | < | "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols); | 
| 533 | > | "\tsingle-processor version of OpenMD.\n", nProcessors, | 
| 534 | > | nGlobalMols); | 
| 535 |  |  | 
| 536 |  | painCave.isFatal = 1; | 
| 537 |  | simError(); | 
| 538 |  | } | 
| 539 |  |  | 
| 524 | – | int seedValue; | 
| 540 |  | Globals * simParams = info->getSimParams(); | 
| 541 | < | SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator | 
| 541 | > | SeqRandNumGen* myRandom; //divide labor does not need Parallel | 
| 542 | > | //random number generator | 
| 543 |  | if (simParams->haveSeed()) { | 
| 544 | < | seedValue = simParams->getSeed(); | 
| 544 | > | int seedValue = simParams->getSeed(); | 
| 545 |  | myRandom = new SeqRandNumGen(seedValue); | 
| 546 |  | }else { | 
| 547 |  | myRandom = new SeqRandNumGen(); | 
| 554 |  | atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0); | 
| 555 |  |  | 
| 556 |  | if (worldRank == 0) { | 
| 557 | < | numerator = info->getNGlobalAtoms(); | 
| 558 | < | denominator = nProcessors; | 
| 559 | < | precast = numerator / denominator; | 
| 560 | < | nTarget = (int)(precast + 0.5); | 
| 557 | > | RealType numerator = info->getNGlobalAtoms(); | 
| 558 | > | RealType denominator = nProcessors; | 
| 559 | > | RealType precast = numerator / denominator; | 
| 560 | > | int nTarget = (int)(precast + 0.5); | 
| 561 |  |  | 
| 562 | < | for(i = 0; i < nGlobalMols; i++) { | 
| 563 | < | done = 0; | 
| 564 | < | loops = 0; | 
| 562 | > | for(int i = 0; i < nGlobalMols; i++) { | 
| 563 | > |  | 
| 564 | > | int done = 0; | 
| 565 | > | int loops = 0; | 
| 566 |  |  | 
| 567 |  | while (!done) { | 
| 568 |  | loops++; | 
| 569 |  |  | 
| 570 |  | // Pick a processor at random | 
| 571 |  |  | 
| 572 | < | which_proc = (int) (myRandom->rand() * nProcessors); | 
| 572 | > | int which_proc = (int) (myRandom->rand() * nProcessors); | 
| 573 |  |  | 
| 574 |  | //get the molecule stamp first | 
| 575 |  | int stampId = info->getMoleculeStampId(i); | 
| 576 |  | MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId); | 
| 577 |  |  | 
| 578 |  | // How many atoms does this processor have so far? | 
| 579 | < | old_atoms = atomsPerProc[which_proc]; | 
| 580 | < | add_atoms = moleculeStamp->getNAtoms(); | 
| 581 | < | new_atoms = old_atoms + add_atoms; | 
| 579 | > | int old_atoms = atomsPerProc[which_proc]; | 
| 580 | > | int add_atoms = moleculeStamp->getNAtoms(); | 
| 581 | > | int new_atoms = old_atoms + add_atoms; | 
| 582 |  |  | 
| 583 |  | // If we've been through this loop too many times, we need | 
| 584 |  | // to just give up and assign the molecule to this processor | 
| 585 |  | // and be done with it. | 
| 586 |  |  | 
| 587 |  | if (loops > 100) { | 
| 588 | + |  | 
| 589 |  | sprintf(painCave.errMsg, | 
| 590 | < | "I've tried 100 times to assign molecule %d to a " | 
| 591 | < | " processor, but can't find a good spot.\n" | 
| 592 | < | "I'm assigning it at random to processor %d.\n", | 
| 590 | > | "There have been 100 attempts to assign molecule %d to an\n" | 
| 591 | > | "\tunderworked processor, but there's no good place to\n" | 
| 592 | > | "\tleave it.  OpenMD is assigning it at random to processor %d.\n", | 
| 593 |  | i, which_proc); | 
| 594 | < |  | 
| 594 | > |  | 
| 595 |  | painCave.isFatal = 0; | 
| 596 | + | painCave.severity = OPENMD_INFO; | 
| 597 |  | simError(); | 
| 598 |  |  | 
| 599 |  | molToProcMap[i] = which_proc; | 
| 622 |  | //           Pacc(x) = exp(- a * x) | 
| 623 |  | // where a = penalty / (average atoms per molecule) | 
| 624 |  |  | 
| 625 | < | x = (RealType)(new_atoms - nTarget); | 
| 626 | < | y = myRandom->rand(); | 
| 625 | > | RealType x = (RealType)(new_atoms - nTarget); | 
| 626 | > | RealType y = myRandom->rand(); | 
| 627 |  |  | 
| 628 |  | if (y < exp(- a * x)) { | 
| 629 |  | molToProcMap[i] = which_proc; | 
| 638 |  | } | 
| 639 |  |  | 
| 640 |  | delete myRandom; | 
| 641 | < |  | 
| 641 | > |  | 
| 642 |  | // Spray out this nonsense to all other processors: | 
| 624 | – |  | 
| 643 |  | MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | 
| 644 | + |  | 
| 645 |  | } else { | 
| 646 |  |  | 
| 647 |  | // Listen to your marching orders from processor 0: | 
| 629 | – |  | 
| 648 |  | MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | 
| 649 | + |  | 
| 650 |  | } | 
| 651 |  |  | 
| 652 |  | info->setMolToProcMap(molToProcMap); | 
| 693 |  | set<AtomType*>::iterator i; | 
| 694 |  | bool hasDirectionalAtoms = false; | 
| 695 |  | bool hasFixedCharge = false; | 
| 696 | < | bool hasMultipoles = false; | 
| 696 | > | bool hasDipoles = false; | 
| 697 | > | bool hasQuadrupoles = false; | 
| 698 |  | bool hasPolarizable = false; | 
| 699 |  | bool hasFluctuatingCharge = false; | 
| 700 |  | bool hasMetallic = false; | 
| 716 |  | if (da.isDirectional()){ | 
| 717 |  | hasDirectionalAtoms = true; | 
| 718 |  | } | 
| 719 | < | if (ma.isMultipole()){ | 
| 720 | < | hasMultipoles = true; | 
| 719 | > | if (ma.isDipole()){ | 
| 720 | > | hasDipoles = true; | 
| 721 |  | } | 
| 722 | + | if (ma.isQuadrupole()){ | 
| 723 | + | hasQuadrupoles = true; | 
| 724 | + | } | 
| 725 |  | if (ea.isEAM() || sca.isSuttonChen()){ | 
| 726 |  | hasMetallic = true; | 
| 727 |  | } | 
| 745 |  | storageLayout |= DataStorage::dslTorque; | 
| 746 |  | } | 
| 747 |  | } | 
| 748 | < | if (hasMultipoles) { | 
| 749 | < | storageLayout |= DataStorage::dslElectroFrame; | 
| 748 | > | if (hasDipoles) { | 
| 749 | > | storageLayout |= DataStorage::dslDipole; | 
| 750 |  | } | 
| 751 | + | if (hasQuadrupoles) { | 
| 752 | + | storageLayout |= DataStorage::dslQuadrupole; | 
| 753 | + | } | 
| 754 |  | if (hasFixedCharge || hasFluctuatingCharge) { | 
| 755 |  | storageLayout |= DataStorage::dslSkippedCharge; | 
| 756 |  | } | 
| 785 |  | } | 
| 786 |  | } | 
| 787 |  |  | 
| 788 | < | if (simParams->getOutputElectricField()) { | 
| 788 | > | if (simParams->getOutputElectricField() | simParams->haveElectricField()) { | 
| 789 |  | storageLayout |= DataStorage::dslElectricField; | 
| 790 |  | } | 
| 791 | + |  | 
| 792 | + | if (simParams->getOutputSitePotential() ) { | 
| 793 | + | storageLayout |= DataStorage::dslSitePotential; | 
| 794 | + | } | 
| 795 | + |  | 
| 796 |  | if (simParams->getOutputFluctuatingCharges()) { | 
| 797 |  | storageLayout |= DataStorage::dslFlucQPosition; | 
| 798 |  | storageLayout |= DataStorage::dslFlucQVelocity; | 
| 799 |  | storageLayout |= DataStorage::dslFlucQForce; | 
| 800 |  | } | 
| 801 |  |  | 
| 802 | + | info->setStorageLayout(storageLayout); | 
| 803 | + |  | 
| 804 |  | return storageLayout; | 
| 805 |  | } | 
| 806 |  |  | 
| 809 |  | Molecule::AtomIterator ai; | 
| 810 |  | Molecule::RigidBodyIterator ri; | 
| 811 |  | Molecule::CutoffGroupIterator ci; | 
| 812 | + | Molecule::BondIterator boi; | 
| 813 | + | Molecule::BendIterator bei; | 
| 814 | + | Molecule::TorsionIterator ti; | 
| 815 | + | Molecule::InversionIterator ii; | 
| 816 |  | Molecule::IntegrableObjectIterator  ioi; | 
| 817 | < | Molecule * mol; | 
| 818 | < | Atom * atom; | 
| 819 | < | RigidBody * rb; | 
| 820 | < | CutoffGroup * cg; | 
| 817 | > | Molecule* mol; | 
| 818 | > | Atom* atom; | 
| 819 | > | RigidBody* rb; | 
| 820 | > | CutoffGroup* cg; | 
| 821 | > | Bond* bond; | 
| 822 | > | Bend* bend; | 
| 823 | > | Torsion* torsion; | 
| 824 | > | Inversion* inversion; | 
| 825 |  | int beginAtomIndex; | 
| 826 |  | int beginRigidBodyIndex; | 
| 827 |  | int beginCutoffGroupIndex; | 
| 828 | + | int beginBondIndex; | 
| 829 | + | int beginBendIndex; | 
| 830 | + | int beginTorsionIndex; | 
| 831 | + | int beginInversionIndex; | 
| 832 |  | int nGlobalAtoms = info->getNGlobalAtoms(); | 
| 833 | + | int nGlobalRigidBodies = info->getNGlobalRigidBodies(); | 
| 834 |  |  | 
| 835 |  | beginAtomIndex = 0; | 
| 836 | < | //rigidbody's index begins right after atom's | 
| 836 | > | // The rigid body indices begin immediately after the atom indices: | 
| 837 |  | beginRigidBodyIndex = info->getNGlobalAtoms(); | 
| 838 |  | beginCutoffGroupIndex = 0; | 
| 839 | < |  | 
| 839 | > | beginBondIndex = 0; | 
| 840 | > | beginBendIndex = 0; | 
| 841 | > | beginTorsionIndex = 0; | 
| 842 | > | beginInversionIndex = 0; | 
| 843 | > |  | 
| 844 |  | for(int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 845 |  |  | 
| 846 |  | #ifdef IS_MPI | 
| 849 |  | // stuff to do if I own this molecule | 
| 850 |  | mol = info->getMoleculeByGlobalIndex(i); | 
| 851 |  |  | 
| 852 | < | //local index(index in DataStorge) of atom is important | 
| 853 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 852 | > | // The local index(index in DataStorge) of the atom is important: | 
| 853 | > | for(atom = mol->beginAtom(ai); atom != NULL; | 
| 854 | > | atom = mol->nextAtom(ai)) { | 
| 855 |  | atom->setGlobalIndex(beginAtomIndex++); | 
| 856 |  | } | 
| 857 |  |  | 
| 860 |  | rb->setGlobalIndex(beginRigidBodyIndex++); | 
| 861 |  | } | 
| 862 |  |  | 
| 863 | < | //local index of cutoff group is trivial, it only depends on | 
| 864 | < | //the order of travesing | 
| 863 | > | // The local index of other objects only depends on the order | 
| 864 | > | // of traversal: | 
| 865 |  | for(cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 866 |  | cg = mol->nextCutoffGroup(ci)) { | 
| 867 |  | cg->setGlobalIndex(beginCutoffGroupIndex++); | 
| 868 |  | } | 
| 869 | + | for(bond = mol->beginBond(boi); bond != NULL; | 
| 870 | + | bond = mol->nextBond(boi)) { | 
| 871 | + | bond->setGlobalIndex(beginBondIndex++); | 
| 872 | + | } | 
| 873 | + | for(bend = mol->beginBend(bei); bend != NULL; | 
| 874 | + | bend = mol->nextBend(bei)) { | 
| 875 | + | bend->setGlobalIndex(beginBendIndex++); | 
| 876 | + | } | 
| 877 | + | for(torsion = mol->beginTorsion(ti); torsion != NULL; | 
| 878 | + | torsion = mol->nextTorsion(ti)) { | 
| 879 | + | torsion->setGlobalIndex(beginTorsionIndex++); | 
| 880 | + | } | 
| 881 | + | for(inversion = mol->beginInversion(ii); inversion != NULL; | 
| 882 | + | inversion = mol->nextInversion(ii)) { | 
| 883 | + | inversion->setGlobalIndex(beginInversionIndex++); | 
| 884 | + | } | 
| 885 |  |  | 
| 886 |  | #ifdef IS_MPI | 
| 887 |  | }  else { | 
| 894 |  | beginAtomIndex += stamp->getNAtoms(); | 
| 895 |  | beginRigidBodyIndex += stamp->getNRigidBodies(); | 
| 896 |  | beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); | 
| 897 | + | beginBondIndex += stamp->getNBonds(); | 
| 898 | + | beginBendIndex += stamp->getNBends(); | 
| 899 | + | beginTorsionIndex += stamp->getNTorsions(); | 
| 900 | + | beginInversionIndex += stamp->getNInversions(); | 
| 901 |  | } | 
| 902 |  | #endif | 
| 903 |  |  | 
| 905 |  |  | 
| 906 |  | //fill globalGroupMembership | 
| 907 |  | std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); | 
| 908 | < | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 909 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 910 | < |  | 
| 908 | > | for(mol = info->beginMolecule(mi); mol != NULL; | 
| 909 | > | mol = info->nextMolecule(mi)) { | 
| 910 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 911 | > | cg = mol->nextCutoffGroup(ci)) { | 
| 912 |  | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 913 |  | globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex(); | 
| 914 |  | } | 
| 923 |  | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | 
| 924 |  | // docs said we could. | 
| 925 |  | std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); | 
| 926 | < | MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, | 
| 926 | > | MPI_Allreduce(&globalGroupMembership[0], | 
| 927 | > | &tmpGroupMembership[0], nGlobalAtoms, | 
| 928 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 929 | + |  | 
| 930 |  | info->setGlobalGroupMembership(tmpGroupMembership); | 
| 931 |  | #else | 
| 932 |  | info->setGlobalGroupMembership(globalGroupMembership); | 
| 933 |  | #endif | 
| 934 |  |  | 
| 935 |  | //fill molMembership | 
| 936 | < | std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); | 
| 936 | > | std::vector<int> globalMolMembership(info->getNGlobalAtoms() + | 
| 937 | > | info->getNGlobalRigidBodies(), 0); | 
| 938 |  |  | 
| 939 | < | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 939 | > | for(mol = info->beginMolecule(mi); mol != NULL; | 
| 940 | > | mol = info->nextMolecule(mi)) { | 
| 941 |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 942 |  | globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 943 |  | } | 
| 944 | + | for (rb = mol->beginRigidBody(ri); rb != NULL; | 
| 945 | + | rb = mol->nextRigidBody(ri)) { | 
| 946 | + | globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 947 | + | } | 
| 948 |  | } | 
| 949 |  |  | 
| 950 |  | #ifdef IS_MPI | 
| 951 | < | std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0); | 
| 952 | < |  | 
| 953 | < | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, | 
| 951 | > | std::vector<int> tmpMolMembership(info->getNGlobalAtoms() + | 
| 952 | > | info->getNGlobalRigidBodies(), 0); | 
| 953 | > | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], | 
| 954 | > | nGlobalAtoms + nGlobalRigidBodies, | 
| 955 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 956 |  |  | 
| 957 |  | info->setGlobalMolMembership(tmpMolMembership); | 
| 963 |  | // here the molecules are listed by their global indices. | 
| 964 |  |  | 
| 965 |  | std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0); | 
| 966 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 966 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 967 | > | mol = info->nextMolecule(mi)) { | 
| 968 |  | nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects(); | 
| 969 |  | } | 
| 970 |  |  | 
| 971 |  | #ifdef IS_MPI | 
| 972 |  | std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0); | 
| 973 |  | MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], | 
| 974 | < | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 974 | > | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 975 |  | #else | 
| 976 |  | std::vector<int> numIntegrableObjectsPerMol = nIOPerMol; | 
| 977 |  | #endif | 
| 985 |  | } | 
| 986 |  |  | 
| 987 |  | std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | 
| 988 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 988 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 989 | > | mol = info->nextMolecule(mi)) { | 
| 990 |  | int myGlobalIndex = mol->getGlobalIndex(); | 
| 991 |  | int globalIO = startingIOIndexForMol[myGlobalIndex]; | 
| 992 |  | for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; | 
| 1002 |  | } | 
| 1003 |  |  | 
| 1004 |  | void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { | 
| 1005 | < |  | 
| 1005 | > |  | 
| 1006 |  | DumpReader reader(info, mdFileName); | 
| 1007 |  | int nframes = reader.getNFrames(); | 
| 1008 | < |  | 
| 1008 | > |  | 
| 1009 |  | if (nframes > 0) { | 
| 1010 |  | reader.readFrame(nframes - 1); | 
| 1011 |  | } else { |