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 1993 by gezelter, Tue Apr 29 17:32:31 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 > #endif                
109          SimplePreprocessor preprocessor;
110 <        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, 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 <        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
117 <
118 <        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
106 <            
107 <                
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 {
109 –        //get stream size
110 –        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
120  
121 +        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
122 +
123 +        //get stream size
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 <        commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
130 <                
129 >        MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
130 >
131          ppStream.str(buf);
132          delete [] buf;
120 –
133        }
134   #endif            
135        // Create a scanner that reads from the input stream
# Line 139 | Line 151 | namespace OpenMD {
151        parser.initializeASTFactory(factory);
152        parser.setASTFactory(&factory);
153        parser.mdfile();
142 –
154        // Create a tree parser that reads information into Globals
155        MDTreeParser treeParser;
156        treeParser.initializeASTFactory(factory);
# Line 229 | Line 240 | namespace OpenMD {
240        simError();
241      }
242  
243 +    simParams->setMDfileVersion(mdFileVersion);
244      return simParams;
245    }
246    
247    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
248                                    bool loadInitCoords) {
249 <
249 >    
250      const int bufferSize = 65535;
251      char buffer[bufferSize];
252      int lineNo = 0;
253      std::string mdRawData;
254      int metaDataBlockStart = -1;
255      int metaDataBlockEnd = -1;
256 <    int i;
257 <    int mdOffset;
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) {
281   #endif
282  
283 <      std::ifstream mdFile_(mdFileName.c_str());
283 >      std::ifstream mdFile_;
284 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
285        
286        if (mdFile_.fail()) {
287          sprintf(painCave.errMsg,
# Line 276 | Line 308 | namespace OpenMD {
308          painCave.isFatal = 1;
309          simError();
310        }
311 +      
312 +      // found the correct opening string, now try to get the file
313 +      // format version number.
314  
315 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
316 +      std::string fileType = tokenizer.nextToken();
317 +      toUpper(fileType);
318 +
319 +      mdFileVersion = 0;
320 +
321 +      if (fileType == "OPENMD") {
322 +        while (tokenizer.hasMoreTokens()) {
323 +          std::string token(tokenizer.nextToken());
324 +          toUpper(token);
325 +          if (token == "VERSION") {
326 +            mdFileVersion = tokenizer.nextTokenAsInt();
327 +            break;
328 +          }
329 +        }
330 +      }
331 +            
332        //scan through the input stream and find MetaData tag        
333        while(mdFile_.getline(buffer, bufferSize)) {
334          ++lineNo;
# Line 317 | Line 369 | namespace OpenMD {
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
# Line 332 | Line 395 | namespace OpenMD {
395      std::stringstream rawMetaDataStream(mdRawData);
396  
397      //parse meta-data file
398 <    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
398 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
399 >                                   metaDataBlockStart + 1);
400      
401      //create the force field
402 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
402 >    ForceField * ff = new ForceField(simParams->getForceField());
403  
404      if (ff == NULL) {
405        sprintf(painCave.errMsg,
# Line 369 | Line 433 | namespace OpenMD {
433      }
434      
435      ff->parse(forcefieldFileName);
372 –    ff->setFortranForceOptions();
436      //create SimInfo
437      SimInfo * info = new SimInfo(ff, simParams);
438  
# Line 387 | Line 450 | namespace OpenMD {
450      //create the molecules
451      createMolecules(info);
452      
453 <    
453 >    //find the storage layout
454 >
455 >    int storageLayout = computeStorageLayout(info);
456 >
457      //allocate memory for DataStorage(circular reference, need to
458      //break it)
459 <    info->setSnapshotManager(new SimSnapshotManager(info));
459 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
460      
461      //set the global index of atoms, rigidbodies and cutoffgroups
462      //(only need to be set once, the global index will never change
# Line 413 | Line 479 | namespace OpenMD {
479      
480      if (loadInitCoords)
481        loadCoordinates(info, mdFileName);    
416 –    
482      return info;
483    }
484    
# Line 448 | Line 513 | namespace OpenMD {
513    
514   #ifdef IS_MPI
515    void SimCreator::divideMolecules(SimInfo *info) {
451 –    RealType numerator;
452 –    RealType denominator;
453 –    RealType precast;
454 –    RealType x;
455 –    RealType y;
516      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;
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,
# Line 477 | Line 530 | namespace OpenMD {
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      
486 –    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();
# Line 500 | Line 554 | namespace OpenMD {
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;
# Line 565 | Line 622 | namespace OpenMD {
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;
# Line 581 | Line 638 | namespace OpenMD {
638        }
639        
640        delete myRandom;
641 <      
641 >
642        // Spray out this nonsense to all other processors:
586 –      
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:
591 –      
648        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
649 +
650      }
651      
652      info->setMolToProcMap(molToProcMap);
# Line 612 | Line 669 | namespace OpenMD {
669   #endif
670          
671          stampId = info->getMoleculeStampId(i);
672 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
673 <                                                   stampId, i, info->getLocalIndexManager());
672 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
673 >                                                   info->getMoleculeStamp(stampId),
674 >                                                   stampId, i,
675 >                                                   info->getLocalIndexManager());
676          
677          info->addMolecule(mol);
678          
# Line 626 | Line 685 | namespace OpenMD {
685      } //end for(int i=0)  
686    }
687      
688 +  int SimCreator::computeStorageLayout(SimInfo* info) {
689 +
690 +    Globals* simParams = info->getSimParams();
691 +    int nRigidBodies = info->getNGlobalRigidBodies();
692 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
693 +    set<AtomType*>::iterator i;
694 +    bool hasDirectionalAtoms = false;
695 +    bool hasFixedCharge = false;
696 +    bool hasDipoles = false;    
697 +    bool hasQuadrupoles = false;    
698 +    bool hasPolarizable = false;    
699 +    bool hasFluctuatingCharge = false;    
700 +    bool hasMetallic = false;
701 +    int storageLayout = 0;
702 +    storageLayout |= DataStorage::dslPosition;
703 +    storageLayout |= DataStorage::dslVelocity;
704 +    storageLayout |= DataStorage::dslForce;
705 +
706 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
707 +
708 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
709 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
710 +      EAMAdapter ea = EAMAdapter( (*i) );
711 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
712 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
713 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
714 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
715 +
716 +      if (da.isDirectional()){
717 +        hasDirectionalAtoms = true;
718 +      }
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 +      }
728 +      if ( fca.isFixedCharge() ){
729 +        hasFixedCharge = true;
730 +      }
731 +      if ( fqa.isFluctuatingCharge() ){
732 +        hasFluctuatingCharge = true;
733 +      }
734 +      if ( pa.isPolarizable() ){
735 +        hasPolarizable = true;
736 +      }
737 +    }
738 +    
739 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
740 +      storageLayout |= DataStorage::dslAmat;
741 +      if(storageLayout & DataStorage::dslVelocity) {
742 +        storageLayout |= DataStorage::dslAngularMomentum;
743 +      }
744 +      if (storageLayout & DataStorage::dslForce) {
745 +        storageLayout |= DataStorage::dslTorque;
746 +      }
747 +    }
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 +    }
757 +    if (hasMetallic) {
758 +      storageLayout |= DataStorage::dslDensity;
759 +      storageLayout |= DataStorage::dslFunctional;
760 +      storageLayout |= DataStorage::dslFunctionalDerivative;
761 +    }
762 +    if (hasPolarizable) {
763 +      storageLayout |= DataStorage::dslElectricField;
764 +    }
765 +    if (hasFluctuatingCharge){
766 +      storageLayout |= DataStorage::dslFlucQPosition;
767 +      if(storageLayout & DataStorage::dslVelocity) {
768 +        storageLayout |= DataStorage::dslFlucQVelocity;
769 +      }
770 +      if (storageLayout & DataStorage::dslForce) {
771 +        storageLayout |= DataStorage::dslFlucQForce;
772 +      }
773 +    }
774 +    
775 +    // if the user has asked for them, make sure we've got the memory for the
776 +    // objects defined.
777 +
778 +    if (simParams->getOutputParticlePotential()) {
779 +      storageLayout |= DataStorage::dslParticlePot;
780 +    }
781 +
782 +    if (simParams->havePrintHeatFlux()) {
783 +      if (simParams->getPrintHeatFlux()) {
784 +        storageLayout |= DataStorage::dslParticlePot;
785 +      }
786 +    }
787 +
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 +
807    void SimCreator::setGlobalIndex(SimInfo *info) {
808      SimInfo::MoleculeIterator mi;
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 <
644 <    /**@todo fixme */
645 < #ifndef IS_MPI
833 >    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
834      
835      beginAtomIndex = 0;
836 <    beginRigidBodyIndex = 0;
836 >    // The rigid body indices begin immediately after the atom indices:
837 >    beginRigidBodyIndex = info->getNGlobalAtoms();
838      beginCutoffGroupIndex = 0;
839 <    
840 < #else
841 <    
842 <    int nproc;
843 <    int myNode;
844 <    
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)) {
839 >    beginBondIndex = 0;
840 >    beginBendIndex = 0;
841 >    beginTorsionIndex = 0;
842 >    beginInversionIndex = 0;
843 >  
844 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
845        
846 <      //local index(index in DataStorge) of atom is important
847 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
848 <        atom->setGlobalIndex(beginAtomIndex++);
846 > #ifdef IS_MPI      
847 >      if (info->getMolToProc(i) == worldRank) {
848 > #endif        
849 >        // stuff to do if I own this molecule
850 >        mol = info->getMoleculeByGlobalIndex(i);
851 >
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 >        
858 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
859 >            rb = mol->nextRigidBody(ri)) {
860 >          rb->setGlobalIndex(beginRigidBodyIndex++);
861 >        }
862 >        
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 {
888 >
889 >        // stuff to do if I don't own this molecule
890 >        
891 >        int stampId = info->getMoleculeStampId(i);
892 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
893 >
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 <      
903 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
904 <          rb = mol->nextRigidBody(ri)) {
905 <        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 <    
902 > #endif          
903 >
904 >    } //end for(int i=0)  
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          }
915          
916        }      
917      }
918 <    
918 >  
919   #ifdef IS_MPI    
920      // Since the globalGroupMembership has been zero filled and we've only
921      // poked values into the atoms we know, we can do an Allreduce
# Line 729 | Line 923 | namespace OpenMD {
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);
# Line 760 | Line 963 | namespace OpenMD {
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    
# Line 781 | Line 985 | namespace OpenMD {
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* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
993 <           integrableObject = mol->nextIntegrableObject(ioi)) {
994 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
995 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
992 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
993 >           sd = mol->nextIntegrableObject(ioi)) {
994 >        sd->setGlobalIntegrableObjectIndex(globalIO);
995 >        IOIndexToIntegrableObject[globalIO] = sd;
996          globalIO++;
997        }
998      }
999 <    
999 >      
1000      info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
1001      
1002    }
1003    
1004    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
800 –    Globals* simParams;
801 –    simParams = info->getSimParams();
1005      
803 –    
1006      DumpReader reader(info, mdFileName);
1007      int nframes = reader.getNFrames();
1008      
# Line 814 | Line 1016 | namespace OpenMD {
1016        painCave.isFatal = 1;
1017        simError();
1018      }
817 –    
1019      //copy the current snapshot to previous snapshot
1020      info->getSnapshotManager()->advance();
1021    }

Diff Legend

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