ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/SimCreator.cpp
(Generate patch)

Comparing trunk/src/brains/SimCreator.cpp (file contents):
Revision 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1976 by gezelter, Wed Mar 12 20:01:15 2014 UTC

# Line 35 | Line 35
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]  Vardeman & Gezelter, in progress (2009).                        
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   */
42  
43   /**
44   * @file SimCreator.cpp
45   * @author tlin
46   * @date 11/03/2004
46 – * @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>
# Line 55 | Line 61
61   #include "brains/SimCreator.hpp"
62   #include "brains/SimSnapshotManager.hpp"
63   #include "io/DumpReader.hpp"
64 < #include "UseTheForce/ForceFieldFactory.hpp"
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"
# Line 75 | Line 82
82   #include "antlr/NoViableAltForCharException.hpp"
83   #include "antlr/NoViableAltException.hpp"
84  
85 < #ifdef IS_MPI
86 < #include "math/ParallelRandNumGen.hpp"
87 < #endif
85 > #include "types/DirectionalAdapter.hpp"
86 > #include "types/MultipoleAdapter.hpp"
87 > #include "types/EAMAdapter.hpp"
88 > #include "types/SuttonChenAdapter.hpp"
89 > #include "types/PolarizableAdapter.hpp"
90 > #include "types/FixedChargeAdapter.hpp"
91 > #include "types/FluctuatingChargeAdapter.hpp"
92  
93 +
94   namespace OpenMD {
95    
96 <  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
96 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
97      Globals* simParams = NULL;
98      try {
99  
# Line 90 | Line 102 | namespace OpenMD {
102   #ifdef IS_MPI            
103        int streamSize;
104        const int masterNode = 0;
105 <      int commStatus;
105 >
106        if (worldRank == masterNode) {
107 < #endif
108 <                
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, 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 <        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
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 <        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
122 <            
123 <                
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 {
109 –        //get stream size
110 –        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
126  
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_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 <        commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
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;
120 –
142        }
143   #endif            
144        // Create a scanner that reads from the input stream
# Line 139 | Line 160 | namespace OpenMD {
160        parser.initializeASTFactory(factory);
161        parser.setASTFactory(&factory);
162        parser.mdfile();
142 –
163        // Create a tree parser that reads information into Globals
164        MDTreeParser treeParser;
165        treeParser.initializeASTFactory(factory);
# Line 229 | Line 249 | namespace OpenMD {
249        simError();
250      }
251  
252 +    simParams->setMDfileVersion(mdFileVersion);
253      return simParams;
254    }
255    
256    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
257                                    bool loadInitCoords) {
258 <
258 >    
259      const int bufferSize = 65535;
260      char buffer[bufferSize];
261      int lineNo = 0;
262      std::string mdRawData;
263      int metaDataBlockStart = -1;
264      int metaDataBlockEnd = -1;
265 <    int i;
266 <    int mdOffset;
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) {
290   #endif
291  
292 <      std::ifstream mdFile_(mdFileName.c_str());
292 >      std::ifstream mdFile_;
293 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
294        
295        if (mdFile_.fail()) {
296          sprintf(painCave.errMsg,
# Line 276 | Line 317 | namespace OpenMD {
317          painCave.isFatal = 1;
318          simError();
319        }
320 +      
321 +      // found the correct opening string, now try to get the file
322 +      // format version number.
323  
324 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
325 +      std::string fileType = tokenizer.nextToken();
326 +      toUpper(fileType);
327 +
328 +      mdFileVersion = 0;
329 +
330 +      if (fileType == "OPENMD") {
331 +        while (tokenizer.hasMoreTokens()) {
332 +          std::string token(tokenizer.nextToken());
333 +          toUpper(token);
334 +          if (token == "VERSION") {
335 +            mdFileVersion = tokenizer.nextTokenAsInt();
336 +            break;
337 +          }
338 +        }
339 +      }
340 +            
341        //scan through the input stream and find MetaData tag        
342        while(mdFile_.getline(buffer, bufferSize)) {
343          ++lineNo;
# Line 317 | Line 378 | namespace OpenMD {
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
# Line 332 | Line 404 | namespace OpenMD {
404      std::stringstream rawMetaDataStream(mdRawData);
405  
406      //parse meta-data file
407 <    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
407 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
408 >                                   metaDataBlockStart + 1);
409      
410      //create the force field
411 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
411 >    ForceField * ff = new ForceField(simParams->getForceField());
412  
413      if (ff == NULL) {
414        sprintf(painCave.errMsg,
# Line 369 | Line 442 | namespace OpenMD {
442      }
443      
444      ff->parse(forcefieldFileName);
372 –    ff->setFortranForceOptions();
445      //create SimInfo
446      SimInfo * info = new SimInfo(ff, simParams);
447  
# Line 387 | Line 459 | namespace OpenMD {
459      //create the molecules
460      createMolecules(info);
461      
462 <    
462 >    //find the storage layout
463 >
464 >    int storageLayout = computeStorageLayout(info);
465 >
466      //allocate memory for DataStorage(circular reference, need to
467      //break it)
468 <    info->setSnapshotManager(new SimSnapshotManager(info));
468 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
469      
470      //set the global index of atoms, rigidbodies and cutoffgroups
471      //(only need to be set once, the global index will never change
# Line 413 | Line 488 | namespace OpenMD {
488      
489      if (loadInitCoords)
490        loadCoordinates(info, mdFileName);    
416 –    
491      return info;
492    }
493    
# Line 448 | Line 522 | namespace OpenMD {
522    
523   #ifdef IS_MPI
524    void SimCreator::divideMolecules(SimInfo *info) {
451 –    RealType numerator;
452 –    RealType denominator;
453 –    RealType precast;
454 –    RealType x;
455 –    RealType y;
525      RealType a;
457 –    int old_atoms;
458 –    int add_atoms;
459 –    int new_atoms;
460 –    int nTarget;
461 –    int done;
462 –    int i;
463 –    int j;
464 –    int loops;
465 –    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,
# Line 477 | Line 540 | namespace OpenMD {
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      
486 –    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();
# Line 500 | Line 564 | namespace OpenMD {
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;
# Line 565 | Line 632 | namespace OpenMD {
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;
# Line 581 | Line 648 | namespace OpenMD {
648        }
649        
650        delete myRandom;
651 <      
651 >
652        // Spray out this nonsense to all other processors:
586 –      
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:
591 –      
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);
# Line 612 | Line 680 | namespace OpenMD {
680   #endif
681          
682          stampId = info->getMoleculeStampId(i);
683 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
684 <                                                   stampId, i, info->getLocalIndexManager());
683 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
684 >                                                   info->getMoleculeStamp(stampId),
685 >                                                   stampId, i,
686 >                                                   info->getLocalIndexManager());
687          
688          info->addMolecule(mol);
689          
# Line 626 | Line 696 | namespace OpenMD {
696      } //end for(int i=0)  
697    }
698      
699 +  int SimCreator::computeStorageLayout(SimInfo* info) {
700 +
701 +    Globals* simParams = info->getSimParams();
702 +    int nRigidBodies = info->getNGlobalRigidBodies();
703 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
704 +    set<AtomType*>::iterator i;
705 +    bool hasDirectionalAtoms = false;
706 +    bool hasFixedCharge = false;
707 +    bool hasDipoles = false;    
708 +    bool hasQuadrupoles = false;    
709 +    bool hasPolarizable = false;    
710 +    bool hasFluctuatingCharge = false;    
711 +    bool hasMetallic = false;
712 +    int storageLayout = 0;
713 +    storageLayout |= DataStorage::dslPosition;
714 +    storageLayout |= DataStorage::dslVelocity;
715 +    storageLayout |= DataStorage::dslForce;
716 +
717 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
718 +
719 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
720 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
721 +      EAMAdapter ea = EAMAdapter( (*i) );
722 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
723 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
724 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
725 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
726 +
727 +      if (da.isDirectional()){
728 +        hasDirectionalAtoms = true;
729 +      }
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 +      }
739 +      if ( fca.isFixedCharge() ){
740 +        hasFixedCharge = true;
741 +      }
742 +      if ( fqa.isFluctuatingCharge() ){
743 +        hasFluctuatingCharge = true;
744 +      }
745 +      if ( pa.isPolarizable() ){
746 +        hasPolarizable = true;
747 +      }
748 +    }
749 +    
750 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
751 +      storageLayout |= DataStorage::dslAmat;
752 +      if(storageLayout & DataStorage::dslVelocity) {
753 +        storageLayout |= DataStorage::dslAngularMomentum;
754 +      }
755 +      if (storageLayout & DataStorage::dslForce) {
756 +        storageLayout |= DataStorage::dslTorque;
757 +      }
758 +    }
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 +    }
768 +    if (hasMetallic) {
769 +      storageLayout |= DataStorage::dslDensity;
770 +      storageLayout |= DataStorage::dslFunctional;
771 +      storageLayout |= DataStorage::dslFunctionalDerivative;
772 +    }
773 +    if (hasPolarizable) {
774 +      storageLayout |= DataStorage::dslElectricField;
775 +    }
776 +    if (hasFluctuatingCharge){
777 +      storageLayout |= DataStorage::dslFlucQPosition;
778 +      if(storageLayout & DataStorage::dslVelocity) {
779 +        storageLayout |= DataStorage::dslFlucQVelocity;
780 +      }
781 +      if (storageLayout & DataStorage::dslForce) {
782 +        storageLayout |= DataStorage::dslFlucQForce;
783 +      }
784 +    }
785 +    
786 +    // if the user has asked for them, make sure we've got the memory for the
787 +    // objects defined.
788 +
789 +    if (simParams->getOutputParticlePotential()) {
790 +      storageLayout |= DataStorage::dslParticlePot;
791 +    }
792 +
793 +    if (simParams->havePrintHeatFlux()) {
794 +      if (simParams->getPrintHeatFlux()) {
795 +        storageLayout |= DataStorage::dslParticlePot;
796 +      }
797 +    }
798 +
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 +
814    void SimCreator::setGlobalIndex(SimInfo *info) {
815      SimInfo::MoleculeIterator mi;
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 <
644 <    /**@todo fixme */
645 < #ifndef IS_MPI
840 >    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
841      
842      beginAtomIndex = 0;
843 <    beginRigidBodyIndex = 0;
843 >    // The rigid body indices begin immediately after the atom indices:
844 >    beginRigidBodyIndex = info->getNGlobalAtoms();
845      beginCutoffGroupIndex = 0;
846 <    
847 < #else
848 <    
849 <    int nproc;
850 <    int myNode;
851 <    
656 <    myNode = worldRank;
657 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
658 <    
659 <    std::vector < int > tmpAtomsInProc(nproc, 0);
660 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
661 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
662 <    std::vector < int > NumAtomsInProc(nproc, 0);
663 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
664 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
665 <    
666 <    tmpAtomsInProc[myNode] = info->getNAtoms();
667 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
668 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
669 <    
670 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
671 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
672 <                  MPI_SUM, MPI_COMM_WORLD);
673 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
674 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
675 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
676 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
677 <    
678 <    beginAtomIndex = 0;
679 <    beginRigidBodyIndex = 0;
680 <    beginCutoffGroupIndex = 0;
681 <    
682 <    for(int i = 0; i < myNode; i++) {
683 <      beginAtomIndex += NumAtomsInProc[i];
684 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
685 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
686 <    }
687 <    
688 < #endif
689 <    
690 <    //rigidbody's index begins right after atom's
691 <    beginRigidBodyIndex += info->getNGlobalAtoms();
692 <    
693 <    for(mol = info->beginMolecule(mi); mol != NULL;
694 <        mol = info->nextMolecule(mi)) {
846 >    beginBondIndex = 0;
847 >    beginBendIndex = 0;
848 >    beginTorsionIndex = 0;
849 >    beginInversionIndex = 0;
850 >  
851 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
852        
853 <      //local index(index in DataStorge) of atom is important
854 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
855 <        atom->setGlobalIndex(beginAtomIndex++);
853 > #ifdef IS_MPI      
854 >      if (info->getMolToProc(i) == worldRank) {
855 > #endif        
856 >        // stuff to do if I own this molecule
857 >        mol = info->getMoleculeByGlobalIndex(i);
858 >
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 >        
865 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
866 >            rb = mol->nextRigidBody(ri)) {
867 >          rb->setGlobalIndex(beginRigidBodyIndex++);
868 >        }
869 >        
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 {
895 >
896 >        // stuff to do if I don't own this molecule
897 >        
898 >        int stampId = info->getMoleculeStampId(i);
899 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
900 >
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 <      
910 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
911 <          rb = mol->nextRigidBody(ri)) {
912 <        rb->setGlobalIndex(beginRigidBodyIndex++);
704 <      }
705 <      
706 <      //local index of cutoff group is trivial, it only depends on the order of travesing
707 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
708 <          cg = mol->nextCutoffGroup(ci)) {
709 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
710 <      }
711 <    }
712 <    
909 > #endif          
910 >
911 >    } //end for(int i=0)  
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          }
922          
923        }      
924      }
925 <    
925 >  
926   #ifdef IS_MPI    
927      // Since the globalGroupMembership has been zero filled and we've only
928      // poked values into the atoms we know, we can do an Allreduce
# Line 729 | Line 930 | namespace OpenMD {
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
# Line 760 | Line 975 | namespace OpenMD {
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    
# Line 781 | Line 999 | namespace OpenMD {
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* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
1007 <           integrableObject = mol->nextIntegrableObject(ioi)) {
1008 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
1009 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
1006 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
1007 >           sd = mol->nextIntegrableObject(ioi)) {
1008 >        sd->setGlobalIntegrableObjectIndex(globalIO);
1009 >        IOIndexToIntegrableObject[globalIO] = sd;
1010          globalIO++;
1011        }
1012      }
1013 <    
1013 >      
1014      info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
1015      
1016    }
1017    
1018    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
800 –    Globals* simParams;
801 –    simParams = info->getSimParams();
1019      
803 –    
1020      DumpReader reader(info, mdFileName);
1021      int nframes = reader.getNFrames();
1022      
# Line 814 | Line 1030 | namespace OpenMD {
1030        painCave.isFatal = 1;
1031        simError();
1032      }
817 –    
1033      //copy the current snapshot to previous snapshot
1034      info->getSnapshotManager()->advance();
1035    }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines