| 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 | > | // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 109 |  | #endif | 
| 110 |  | SimplePreprocessor preprocessor; | 
| 111 | < | preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, | 
| 112 | < | ppStream); | 
| 111 | > | preprocessor.preprocess(rawMetaDataStream, filename, | 
| 112 | > | startOfMetaDataBlock, ppStream); | 
| 113 |  |  | 
| 114 |  | #ifdef IS_MPI | 
| 115 | < | //brocasting the stream size | 
| 115 | > | //broadcasting the stream size | 
| 116 |  | streamSize = ppStream.str().size() +1; | 
| 117 | < | MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 118 | < | MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), | 
| 119 | < | streamSize, MPI::CHAR, masterNode); | 
| 120 | < |  | 
| 117 | > | MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 118 | > | MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), | 
| 119 | > | streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 120 | > |  | 
| 121 | > | // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 122 | > | // MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), | 
| 123 | > | //                       streamSize, MPI::CHAR, masterNode); | 
| 124 | > |  | 
| 125 |  | } else { | 
| 126 |  |  | 
| 127 | < | MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 127 | > | MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 128 | > | // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); | 
| 129 |  |  | 
| 130 |  | //get stream size | 
| 131 | < | MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 132 | < |  | 
| 131 | > | MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); | 
| 132 | > | // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); | 
| 133 |  | char* buf = new char[streamSize]; | 
| 134 |  | assert(buf); | 
| 135 |  |  | 
| 136 |  | //receive file content | 
| 137 | < | MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); | 
| 138 | < |  | 
| 137 | > | MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 138 | > | // MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); | 
| 139 | > |  | 
| 140 |  | ppStream.str(buf); | 
| 141 |  | delete [] buf; | 
| 133 | – |  | 
| 142 |  | } | 
| 143 |  | #endif | 
| 144 |  | // Create a scanner that reads from the input stream | 
| 160 |  | parser.initializeASTFactory(factory); | 
| 161 |  | parser.setASTFactory(&factory); | 
| 162 |  | parser.mdfile(); | 
| 155 | – |  | 
| 163 |  | // Create a tree parser that reads information into Globals | 
| 164 |  | MDTreeParser treeParser; | 
| 165 |  | treeParser.initializeASTFactory(factory); | 
| 262 |  | std::string mdRawData; | 
| 263 |  | int metaDataBlockStart = -1; | 
| 264 |  | int metaDataBlockEnd = -1; | 
| 265 | < | int i; | 
| 266 | < | streamoff mdOffset(0); | 
| 265 | > | int i, j; | 
| 266 | > | streamoff mdOffset; | 
| 267 |  | int mdFileVersion; | 
| 268 |  |  | 
| 269 | + | // Create a string for embedding the version information in the MetaData | 
| 270 | + | std::string version; | 
| 271 | + | version.assign("## Last run using OpenMD Version: "); | 
| 272 | + | version.append(OPENMD_VERSION_MAJOR); | 
| 273 | + | version.append("."); | 
| 274 | + | version.append(OPENMD_VERSION_MINOR); | 
| 275 |  |  | 
| 276 | + | std::string svnrev(g_REVISION, strnlen(g_REVISION, 20)); | 
| 277 | + | //convert a macro from compiler to a string in c++ | 
| 278 | + | // STR_DEFINE(svnrev, SVN_REV ); | 
| 279 | + | version.append(" Revision: "); | 
| 280 | + | // If there's no SVN revision, just call this the RELEASE revision. | 
| 281 | + | if (!svnrev.empty()) { | 
| 282 | + | version.append(svnrev); | 
| 283 | + | } else { | 
| 284 | + | version.append("RELEASE"); | 
| 285 | + | } | 
| 286 | + |  | 
| 287 |  | #ifdef IS_MPI | 
| 288 |  | const int masterNode = 0; | 
| 289 |  | if (worldRank == masterNode) { | 
| 378 |  |  | 
| 379 |  | mdRawData.clear(); | 
| 380 |  |  | 
| 381 | + | bool foundVersion = false; | 
| 382 | + |  | 
| 383 |  | for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) { | 
| 384 |  | mdFile_.getline(buffer, bufferSize); | 
| 385 | < | mdRawData += buffer; | 
| 385 | > | std::string line = trimLeftCopy(buffer); | 
| 386 | > | j = CaseInsensitiveFind(line, "## Last run using OpenMD Version"); | 
| 387 | > | if (static_cast<size_t>(j) != string::npos) { | 
| 388 | > | foundVersion = true; | 
| 389 | > | mdRawData += version; | 
| 390 | > | } else { | 
| 391 | > | mdRawData += buffer; | 
| 392 | > | } | 
| 393 |  | mdRawData += "\n"; | 
| 394 |  | } | 
| 395 | < |  | 
| 395 | > |  | 
| 396 | > | if (!foundVersion) mdRawData += version + "\n"; | 
| 397 | > |  | 
| 398 |  | mdFile_.close(); | 
| 399 |  |  | 
| 400 |  | #ifdef IS_MPI | 
| 522 |  |  | 
| 523 |  | #ifdef IS_MPI | 
| 524 |  | void SimCreator::divideMolecules(SimInfo *info) { | 
| 490 | – | RealType numerator; | 
| 491 | – | RealType denominator; | 
| 492 | – | RealType precast; | 
| 493 | – | RealType x; | 
| 494 | – | RealType y; | 
| 525 |  | 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; | 
| 526 |  | int nProcessors; | 
| 527 |  | std::vector<int> atomsPerProc; | 
| 528 |  | int nGlobalMols = info->getNGlobalMolecules(); | 
| 529 | < | std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition: | 
| 529 | > | std::vector<int> molToProcMap(nGlobalMols, -1); // default to an | 
| 530 | > | // error | 
| 531 | > | // condition: | 
| 532 |  |  | 
| 533 | < | MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); | 
| 533 | > | MPI_Comm_size( MPI_COMM_WORLD, &nProcessors); | 
| 534 | > | //nProcessors = MPI::COMM_WORLD.Get_size(); | 
| 535 |  |  | 
| 536 |  | if (nProcessors > nGlobalMols) { | 
| 537 |  | sprintf(painCave.errMsg, | 
| 540 |  | "\tthe number of molecules.  This will not result in a \n" | 
| 541 |  | "\tusable division of atoms for force decomposition.\n" | 
| 542 |  | "\tEither try a smaller number of processors, or run the\n" | 
| 543 | < | "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols); | 
| 543 | > | "\tsingle-processor version of OpenMD.\n", nProcessors, | 
| 544 | > | nGlobalMols); | 
| 545 |  |  | 
| 546 |  | painCave.isFatal = 1; | 
| 547 |  | simError(); | 
| 548 |  | } | 
| 549 |  |  | 
| 524 | – | int seedValue; | 
| 550 |  | Globals * simParams = info->getSimParams(); | 
| 551 | < | SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator | 
| 551 | > | SeqRandNumGen* myRandom; //divide labor does not need Parallel | 
| 552 | > | //random number generator | 
| 553 |  | if (simParams->haveSeed()) { | 
| 554 | < | seedValue = simParams->getSeed(); | 
| 554 | > | int seedValue = simParams->getSeed(); | 
| 555 |  | myRandom = new SeqRandNumGen(seedValue); | 
| 556 |  | }else { | 
| 557 |  | myRandom = new SeqRandNumGen(); | 
| 564 |  | atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0); | 
| 565 |  |  | 
| 566 |  | if (worldRank == 0) { | 
| 567 | < | numerator = info->getNGlobalAtoms(); | 
| 568 | < | denominator = nProcessors; | 
| 569 | < | precast = numerator / denominator; | 
| 570 | < | nTarget = (int)(precast + 0.5); | 
| 567 | > | RealType numerator = info->getNGlobalAtoms(); | 
| 568 | > | RealType denominator = nProcessors; | 
| 569 | > | RealType precast = numerator / denominator; | 
| 570 | > | int nTarget = (int)(precast + 0.5); | 
| 571 |  |  | 
| 572 | < | for(i = 0; i < nGlobalMols; i++) { | 
| 573 | < | done = 0; | 
| 574 | < | loops = 0; | 
| 572 | > | for(int i = 0; i < nGlobalMols; i++) { | 
| 573 | > |  | 
| 574 | > | int done = 0; | 
| 575 | > | int loops = 0; | 
| 576 |  |  | 
| 577 |  | while (!done) { | 
| 578 |  | loops++; | 
| 579 |  |  | 
| 580 |  | // Pick a processor at random | 
| 581 |  |  | 
| 582 | < | which_proc = (int) (myRandom->rand() * nProcessors); | 
| 582 | > | int which_proc = (int) (myRandom->rand() * nProcessors); | 
| 583 |  |  | 
| 584 |  | //get the molecule stamp first | 
| 585 |  | int stampId = info->getMoleculeStampId(i); | 
| 586 |  | MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId); | 
| 587 |  |  | 
| 588 |  | // How many atoms does this processor have so far? | 
| 589 | < | old_atoms = atomsPerProc[which_proc]; | 
| 590 | < | add_atoms = moleculeStamp->getNAtoms(); | 
| 591 | < | new_atoms = old_atoms + add_atoms; | 
| 589 | > | int old_atoms = atomsPerProc[which_proc]; | 
| 590 | > | int add_atoms = moleculeStamp->getNAtoms(); | 
| 591 | > | int new_atoms = old_atoms + add_atoms; | 
| 592 |  |  | 
| 593 |  | // If we've been through this loop too many times, we need | 
| 594 |  | // to just give up and assign the molecule to this processor | 
| 595 |  | // and be done with it. | 
| 596 |  |  | 
| 597 |  | if (loops > 100) { | 
| 598 | + |  | 
| 599 |  | sprintf(painCave.errMsg, | 
| 600 | < | "I've tried 100 times to assign molecule %d to a " | 
| 601 | < | " processor, but can't find a good spot.\n" | 
| 602 | < | "I'm assigning it at random to processor %d.\n", | 
| 600 | > | "There have been 100 attempts to assign molecule %d to an\n" | 
| 601 | > | "\tunderworked processor, but there's no good place to\n" | 
| 602 | > | "\tleave it.  OpenMD is assigning it at random to processor %d.\n", | 
| 603 |  | i, which_proc); | 
| 604 | < |  | 
| 604 | > |  | 
| 605 |  | painCave.isFatal = 0; | 
| 606 | + | painCave.severity = OPENMD_INFO; | 
| 607 |  | simError(); | 
| 608 |  |  | 
| 609 |  | molToProcMap[i] = which_proc; | 
| 632 |  | //           Pacc(x) = exp(- a * x) | 
| 633 |  | // where a = penalty / (average atoms per molecule) | 
| 634 |  |  | 
| 635 | < | x = (RealType)(new_atoms - nTarget); | 
| 636 | < | y = myRandom->rand(); | 
| 635 | > | RealType x = (RealType)(new_atoms - nTarget); | 
| 636 | > | RealType y = myRandom->rand(); | 
| 637 |  |  | 
| 638 |  | if (y < exp(- a * x)) { | 
| 639 |  | molToProcMap[i] = which_proc; | 
| 648 |  | } | 
| 649 |  |  | 
| 650 |  | delete myRandom; | 
| 651 | < |  | 
| 651 | > |  | 
| 652 |  | // Spray out this nonsense to all other processors: | 
| 624 | – |  | 
| 653 |  | MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | 
| 654 | + | // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); | 
| 655 |  | } else { | 
| 656 |  |  | 
| 657 |  | // Listen to your marching orders from processor 0: | 
| 629 | – |  | 
| 658 |  | MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | 
| 659 | + | // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); | 
| 660 | + |  | 
| 661 |  | } | 
| 662 |  |  | 
| 663 |  | info->setMolToProcMap(molToProcMap); | 
| 704 |  | set<AtomType*>::iterator i; | 
| 705 |  | bool hasDirectionalAtoms = false; | 
| 706 |  | bool hasFixedCharge = false; | 
| 707 | < | bool hasMultipoles = false; | 
| 707 | > | bool hasDipoles = false; | 
| 708 | > | bool hasQuadrupoles = false; | 
| 709 |  | bool hasPolarizable = false; | 
| 710 |  | bool hasFluctuatingCharge = false; | 
| 711 |  | bool hasMetallic = false; | 
| 727 |  | if (da.isDirectional()){ | 
| 728 |  | hasDirectionalAtoms = true; | 
| 729 |  | } | 
| 730 | < | if (ma.isMultipole()){ | 
| 731 | < | hasMultipoles = true; | 
| 730 | > | if (ma.isDipole()){ | 
| 731 | > | hasDipoles = true; | 
| 732 |  | } | 
| 733 | + | if (ma.isQuadrupole()){ | 
| 734 | + | hasQuadrupoles = true; | 
| 735 | + | } | 
| 736 |  | if (ea.isEAM() || sca.isSuttonChen()){ | 
| 737 |  | hasMetallic = true; | 
| 738 |  | } | 
| 756 |  | storageLayout |= DataStorage::dslTorque; | 
| 757 |  | } | 
| 758 |  | } | 
| 759 | < | if (hasMultipoles) { | 
| 760 | < | storageLayout |= DataStorage::dslElectroFrame; | 
| 759 | > | if (hasDipoles) { | 
| 760 | > | storageLayout |= DataStorage::dslDipole; | 
| 761 |  | } | 
| 762 | + | if (hasQuadrupoles) { | 
| 763 | + | storageLayout |= DataStorage::dslQuadrupole; | 
| 764 | + | } | 
| 765 |  | if (hasFixedCharge || hasFluctuatingCharge) { | 
| 766 |  | storageLayout |= DataStorage::dslSkippedCharge; | 
| 767 |  | } | 
| 796 |  | } | 
| 797 |  | } | 
| 798 |  |  | 
| 799 | < | if (simParams->getOutputElectricField()) { | 
| 799 | > | if (simParams->getOutputElectricField() | simParams->haveElectricField()) { | 
| 800 |  | storageLayout |= DataStorage::dslElectricField; | 
| 801 |  | } | 
| 802 | + |  | 
| 803 |  | if (simParams->getOutputFluctuatingCharges()) { | 
| 804 |  | storageLayout |= DataStorage::dslFlucQPosition; | 
| 805 |  | storageLayout |= DataStorage::dslFlucQVelocity; | 
| 806 |  | storageLayout |= DataStorage::dslFlucQForce; | 
| 807 |  | } | 
| 808 |  |  | 
| 809 | + | info->setStorageLayout(storageLayout); | 
| 810 | + |  | 
| 811 |  | return storageLayout; | 
| 812 |  | } | 
| 813 |  |  | 
| 816 |  | Molecule::AtomIterator ai; | 
| 817 |  | Molecule::RigidBodyIterator ri; | 
| 818 |  | Molecule::CutoffGroupIterator ci; | 
| 819 | + | Molecule::BondIterator boi; | 
| 820 | + | Molecule::BendIterator bei; | 
| 821 | + | Molecule::TorsionIterator ti; | 
| 822 | + | Molecule::InversionIterator ii; | 
| 823 |  | Molecule::IntegrableObjectIterator  ioi; | 
| 824 | < | Molecule * mol; | 
| 825 | < | Atom * atom; | 
| 826 | < | RigidBody * rb; | 
| 827 | < | CutoffGroup * cg; | 
| 824 | > | Molecule* mol; | 
| 825 | > | Atom* atom; | 
| 826 | > | RigidBody* rb; | 
| 827 | > | CutoffGroup* cg; | 
| 828 | > | Bond* bond; | 
| 829 | > | Bend* bend; | 
| 830 | > | Torsion* torsion; | 
| 831 | > | Inversion* inversion; | 
| 832 |  | int beginAtomIndex; | 
| 833 |  | int beginRigidBodyIndex; | 
| 834 |  | int beginCutoffGroupIndex; | 
| 835 | + | int beginBondIndex; | 
| 836 | + | int beginBendIndex; | 
| 837 | + | int beginTorsionIndex; | 
| 838 | + | int beginInversionIndex; | 
| 839 |  | int nGlobalAtoms = info->getNGlobalAtoms(); | 
| 840 | + | int nGlobalRigidBodies = info->getNGlobalRigidBodies(); | 
| 841 |  |  | 
| 842 |  | beginAtomIndex = 0; | 
| 843 | < | //rigidbody's index begins right after atom's | 
| 843 | > | // The rigid body indices begin immediately after the atom indices: | 
| 844 |  | beginRigidBodyIndex = info->getNGlobalAtoms(); | 
| 845 |  | beginCutoffGroupIndex = 0; | 
| 846 | < |  | 
| 846 | > | beginBondIndex = 0; | 
| 847 | > | beginBendIndex = 0; | 
| 848 | > | beginTorsionIndex = 0; | 
| 849 | > | beginInversionIndex = 0; | 
| 850 | > |  | 
| 851 |  | for(int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 852 |  |  | 
| 853 |  | #ifdef IS_MPI | 
| 856 |  | // stuff to do if I own this molecule | 
| 857 |  | mol = info->getMoleculeByGlobalIndex(i); | 
| 858 |  |  | 
| 859 | < | //local index(index in DataStorge) of atom is important | 
| 860 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 859 | > | // The local index(index in DataStorge) of the atom is important: | 
| 860 | > | for(atom = mol->beginAtom(ai); atom != NULL; | 
| 861 | > | atom = mol->nextAtom(ai)) { | 
| 862 |  | atom->setGlobalIndex(beginAtomIndex++); | 
| 863 |  | } | 
| 864 |  |  | 
| 867 |  | rb->setGlobalIndex(beginRigidBodyIndex++); | 
| 868 |  | } | 
| 869 |  |  | 
| 870 | < | //local index of cutoff group is trivial, it only depends on | 
| 871 | < | //the order of travesing | 
| 870 | > | // The local index of other objects only depends on the order | 
| 871 | > | // of traversal: | 
| 872 |  | for(cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 873 |  | cg = mol->nextCutoffGroup(ci)) { | 
| 874 |  | cg->setGlobalIndex(beginCutoffGroupIndex++); | 
| 875 |  | } | 
| 876 | + | for(bond = mol->beginBond(boi); bond != NULL; | 
| 877 | + | bond = mol->nextBond(boi)) { | 
| 878 | + | bond->setGlobalIndex(beginBondIndex++); | 
| 879 | + | } | 
| 880 | + | for(bend = mol->beginBend(bei); bend != NULL; | 
| 881 | + | bend = mol->nextBend(bei)) { | 
| 882 | + | bend->setGlobalIndex(beginBendIndex++); | 
| 883 | + | } | 
| 884 | + | for(torsion = mol->beginTorsion(ti); torsion != NULL; | 
| 885 | + | torsion = mol->nextTorsion(ti)) { | 
| 886 | + | torsion->setGlobalIndex(beginTorsionIndex++); | 
| 887 | + | } | 
| 888 | + | for(inversion = mol->beginInversion(ii); inversion != NULL; | 
| 889 | + | inversion = mol->nextInversion(ii)) { | 
| 890 | + | inversion->setGlobalIndex(beginInversionIndex++); | 
| 891 | + | } | 
| 892 |  |  | 
| 893 |  | #ifdef IS_MPI | 
| 894 |  | }  else { | 
| 901 |  | beginAtomIndex += stamp->getNAtoms(); | 
| 902 |  | beginRigidBodyIndex += stamp->getNRigidBodies(); | 
| 903 |  | beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); | 
| 904 | + | beginBondIndex += stamp->getNBonds(); | 
| 905 | + | beginBendIndex += stamp->getNBends(); | 
| 906 | + | beginTorsionIndex += stamp->getNTorsions(); | 
| 907 | + | beginInversionIndex += stamp->getNInversions(); | 
| 908 |  | } | 
| 909 |  | #endif | 
| 910 |  |  | 
| 912 |  |  | 
| 913 |  | //fill globalGroupMembership | 
| 914 |  | std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); | 
| 915 | < | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 916 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 917 | < |  | 
| 915 | > | for(mol = info->beginMolecule(mi); mol != NULL; | 
| 916 | > | mol = info->nextMolecule(mi)) { | 
| 917 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; | 
| 918 | > | cg = mol->nextCutoffGroup(ci)) { | 
| 919 |  | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 920 |  | globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex(); | 
| 921 |  | } | 
| 930 |  | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | 
| 931 |  | // docs said we could. | 
| 932 |  | std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); | 
| 933 | < | MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, | 
| 933 | > | MPI_Allreduce(&globalGroupMembership[0], | 
| 934 | > | &tmpGroupMembership[0], nGlobalAtoms, | 
| 935 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 936 | + | // MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0], | 
| 937 | + | //                           &tmpGroupMembership[0], nGlobalAtoms, | 
| 938 | + | //                           MPI::INT, MPI::SUM); | 
| 939 |  | info->setGlobalGroupMembership(tmpGroupMembership); | 
| 940 |  | #else | 
| 941 |  | info->setGlobalGroupMembership(globalGroupMembership); | 
| 942 |  | #endif | 
| 943 |  |  | 
| 944 |  | //fill molMembership | 
| 945 | < | std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); | 
| 945 | > | std::vector<int> globalMolMembership(info->getNGlobalAtoms() + | 
| 946 | > | info->getNGlobalRigidBodies(), 0); | 
| 947 |  |  | 
| 948 | < | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 948 | > | for(mol = info->beginMolecule(mi); mol != NULL; | 
| 949 | > | mol = info->nextMolecule(mi)) { | 
| 950 |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 951 |  | globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 952 |  | } | 
| 953 | + | for (rb = mol->beginRigidBody(ri); rb != NULL; | 
| 954 | + | rb = mol->nextRigidBody(ri)) { | 
| 955 | + | globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 956 | + | } | 
| 957 |  | } | 
| 958 |  |  | 
| 959 |  | #ifdef IS_MPI | 
| 960 | < | std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0); | 
| 961 | < |  | 
| 962 | < | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, | 
| 960 | > | std::vector<int> tmpMolMembership(info->getNGlobalAtoms() + | 
| 961 | > | info->getNGlobalRigidBodies(), 0); | 
| 962 | > | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], | 
| 963 | > | nGlobalAtoms + nGlobalRigidBodies, | 
| 964 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 965 | + | // MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0], | 
| 966 | + | //                           nGlobalAtoms + nGlobalRigidBodies, | 
| 967 | + | //                           MPI::INT, MPI::SUM); | 
| 968 |  |  | 
| 969 |  | info->setGlobalMolMembership(tmpMolMembership); | 
| 970 |  | #else | 
| 975 |  | // here the molecules are listed by their global indices. | 
| 976 |  |  | 
| 977 |  | std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0); | 
| 978 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 978 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 979 | > | mol = info->nextMolecule(mi)) { | 
| 980 |  | nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects(); | 
| 981 |  | } | 
| 982 |  |  | 
| 983 |  | #ifdef IS_MPI | 
| 984 |  | std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0); | 
| 985 |  | MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], | 
| 986 | < | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 986 | > | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 987 | > | // MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], | 
| 988 | > | //                           info->getNGlobalMolecules(), MPI::INT, MPI::SUM); | 
| 989 |  | #else | 
| 990 |  | std::vector<int> numIntegrableObjectsPerMol = nIOPerMol; | 
| 991 |  | #endif | 
| 999 |  | } | 
| 1000 |  |  | 
| 1001 |  | std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | 
| 1002 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 1002 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 1003 | > | mol = info->nextMolecule(mi)) { | 
| 1004 |  | int myGlobalIndex = mol->getGlobalIndex(); | 
| 1005 |  | int globalIO = startingIOIndexForMol[myGlobalIndex]; | 
| 1006 |  | for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; | 
| 1016 |  | } | 
| 1017 |  |  | 
| 1018 |  | void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { | 
| 1019 | < |  | 
| 1019 | > |  | 
| 1020 |  | DumpReader reader(info, mdFileName); | 
| 1021 |  | int nframes = reader.getNFrames(); | 
| 1022 | < |  | 
| 1022 | > |  | 
| 1023 |  | if (nframes > 0) { | 
| 1024 |  | reader.readFrame(nframes - 1); | 
| 1025 |  | } else { |