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 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1810 by gezelter, Thu Nov 8 14:23:43 2012 UTC

# Line 1 | Line 1
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
# Line 36 | Line 36
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).                        
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   /**
# Line 55 | Line 56
56   #include "brains/SimCreator.hpp"
57   #include "brains/SimSnapshotManager.hpp"
58   #include "io/DumpReader.hpp"
59 < #include "UseTheForce/ForceFieldFactory.hpp"
59 > #include "brains/ForceField.hpp"
60   #include "utils/simError.h"
61   #include "utils/StringUtils.hpp"
62   #include "math/SeqRandNumGen.hpp"
# Line 75 | Line 76
76   #include "antlr/NoViableAltForCharException.hpp"
77   #include "antlr/NoViableAltException.hpp"
78  
79 + #include "types/DirectionalAdapter.hpp"
80 + #include "types/MultipoleAdapter.hpp"
81 + #include "types/EAMAdapter.hpp"
82 + #include "types/SuttonChenAdapter.hpp"
83 + #include "types/PolarizableAdapter.hpp"
84 + #include "types/FixedChargeAdapter.hpp"
85 + #include "types/FluctuatingChargeAdapter.hpp"
86 +
87   #ifdef IS_MPI
88 + #include "mpi.h"
89   #include "math/ParallelRandNumGen.hpp"
90   #endif
91  
92   namespace OpenMD {
93    
94 <  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
94 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
95      Globals* simParams = NULL;
96      try {
97  
# Line 90 | Line 100 | namespace OpenMD {
100   #ifdef IS_MPI            
101        int streamSize;
102        const int masterNode = 0;
103 <      int commStatus;
103 >
104        if (worldRank == masterNode) {
105 < #endif
106 <                
105 >        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
106 > #endif                
107          SimplePreprocessor preprocessor;
108 <        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
108 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock,
109 >                                ppStream);
110                  
111   #ifdef IS_MPI            
112          //brocasting the stream size
113          streamSize = ppStream.str().size() +1;
114 <        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
115 <
116 <        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
106 <            
114 >        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
115 >        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
116 >                              streamSize, MPI::CHAR, masterNode);
117                  
118        } else {
119 +
120 +        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
121 +
122          //get stream size
123 <        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
123 >        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
124  
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);
129 >        MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode);
130                  
131          ppStream.str(buf);
132          delete [] buf;
# Line 229 | Line 242 | namespace OpenMD {
242        simError();
243      }
244  
245 +    simParams->setMDfileVersion(mdFileVersion);
246      return simParams;
247    }
248    
249    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
250                                    bool loadInitCoords) {
251 <
251 >    
252      const int bufferSize = 65535;
253      char buffer[bufferSize];
254      int lineNo = 0;
255      std::string mdRawData;
256      int metaDataBlockStart = -1;
257      int metaDataBlockEnd = -1;
258 <    int i;
259 <    int mdOffset;
258 >    int i, j;
259 >    streamoff mdOffset(0);
260 >    int mdFileVersion;
261 >
262 >    // Create a string for embedding the version information in the MetaData
263 >    std::string version;
264 >    version.assign("## Last run using OpenMD Version: ");
265 >    version.append(OPENMD_VERSION_MAJOR);
266 >    version.append(".");
267 >    version.append(OPENMD_VERSION_MINOR);
268  
269 +    std::string svnrev;
270 +    //convert a macro from compiler to a string in c++
271 +    STR_DEFINE(svnrev, SVN_REV );
272 +    version.append(" Revision: ");
273 +    // If there's no SVN revision, just call this the RELEASE revision.
274 +    if (!svnrev.empty()) {
275 +      version.append(svnrev);
276 +    } else {
277 +      version.append("RELEASE");
278 +    }
279 +  
280   #ifdef IS_MPI            
281      const int masterNode = 0;
282      if (worldRank == masterNode) {
283   #endif
284  
285 <      std::ifstream mdFile_(mdFileName.c_str());
285 >      std::ifstream mdFile_;
286 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
287        
288        if (mdFile_.fail()) {
289          sprintf(painCave.errMsg,
# Line 276 | Line 310 | namespace OpenMD {
310          painCave.isFatal = 1;
311          simError();
312        }
313 +      
314 +      // found the correct opening string, now try to get the file
315 +      // format version number.
316  
317 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
318 +      std::string fileType = tokenizer.nextToken();
319 +      toUpper(fileType);
320 +
321 +      mdFileVersion = 0;
322 +
323 +      if (fileType == "OPENMD") {
324 +        while (tokenizer.hasMoreTokens()) {
325 +          std::string token(tokenizer.nextToken());
326 +          toUpper(token);
327 +          if (token == "VERSION") {
328 +            mdFileVersion = tokenizer.nextTokenAsInt();
329 +            break;
330 +          }
331 +        }
332 +      }
333 +            
334        //scan through the input stream and find MetaData tag        
335        while(mdFile_.getline(buffer, bufferSize)) {
336          ++lineNo;
# Line 317 | Line 371 | namespace OpenMD {
371  
372        mdRawData.clear();
373  
374 +      bool foundVersion = false;
375 +
376        for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
377          mdFile_.getline(buffer, bufferSize);
378 <        mdRawData += buffer;
378 >        std::string line = trimLeftCopy(buffer);
379 >        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version");
380 >        if (static_cast<size_t>(j) != string::npos) {
381 >          foundVersion = true;
382 >          mdRawData += version;
383 >        } else {
384 >          mdRawData += buffer;
385 >        }
386          mdRawData += "\n";
387        }
388 <
388 >      
389 >      if (!foundVersion) mdRawData += version + "\n";
390 >      
391        mdFile_.close();
392  
393   #ifdef IS_MPI
# Line 332 | Line 397 | namespace OpenMD {
397      std::stringstream rawMetaDataStream(mdRawData);
398  
399      //parse meta-data file
400 <    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
400 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
401 >                                   metaDataBlockStart + 1);
402      
403      //create the force field
404 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
404 >    ForceField * ff = new ForceField(simParams->getForceField());
405  
406      if (ff == NULL) {
407        sprintf(painCave.errMsg,
# Line 369 | Line 435 | namespace OpenMD {
435      }
436      
437      ff->parse(forcefieldFileName);
372    ff->setFortranForceOptions();
438      //create SimInfo
439      SimInfo * info = new SimInfo(ff, simParams);
440  
# Line 387 | Line 452 | namespace OpenMD {
452      //create the molecules
453      createMolecules(info);
454      
455 <    
455 >    //find the storage layout
456 >
457 >    int storageLayout = computeStorageLayout(info);
458 >
459      //allocate memory for DataStorage(circular reference, need to
460      //break it)
461 <    info->setSnapshotManager(new SimSnapshotManager(info));
461 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
462      
463      //set the global index of atoms, rigidbodies and cutoffgroups
464      //(only need to be set once, the global index will never change
# Line 413 | Line 481 | namespace OpenMD {
481      
482      if (loadInitCoords)
483        loadCoordinates(info, mdFileName);    
416    
484      return info;
485    }
486    
# Line 460 | Line 527 | namespace OpenMD {
527      int nTarget;
528      int done;
529      int i;
463    int j;
530      int loops;
531      int which_proc;
532      int nProcessors;
# Line 468 | Line 534 | namespace OpenMD {
534      int nGlobalMols = info->getNGlobalMolecules();
535      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
536      
537 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
537 >    nProcessors = MPI::COMM_WORLD.Get_size();
538      
539      if (nProcessors > nGlobalMols) {
540        sprintf(painCave.errMsg,
# Line 506 | Line 572 | namespace OpenMD {
572        nTarget = (int)(precast + 0.5);
573        
574        for(i = 0; i < nGlobalMols; i++) {
575 +
576          done = 0;
577          loops = 0;
578          
# Line 530 | Line 597 | namespace OpenMD {
597            // and be done with it.
598            
599            if (loops > 100) {
600 +
601              sprintf(painCave.errMsg,
602 <                    "I've tried 100 times to assign molecule %d to a "
603 <                    " processor, but can't find a good spot.\n"
604 <                    "I'm assigning it at random to processor %d.\n",
602 >                    "There have been 100 attempts to assign molecule %d to an\n"
603 >                    "\tunderworked processor, but there's no good place to\n"
604 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
605                      i, which_proc);
606 <            
606 >          
607              painCave.isFatal = 0;
608 +            painCave.severity = OPENMD_INFO;
609              simError();
610              
611              molToProcMap[i] = which_proc;
# Line 581 | Line 650 | namespace OpenMD {
650        }
651        
652        delete myRandom;
653 <      
653 >
654        // Spray out this nonsense to all other processors:
655 <      
587 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
655 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
656      } else {
657        
658        // Listen to your marching orders from processor 0:
659 <      
660 <      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 625 | Line 695 | namespace OpenMD {
695        
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 hasMultipoles = false;    
708 +    bool hasPolarizable = false;    
709 +    bool hasFluctuatingCharge = false;    
710 +    bool hasMetallic = false;
711 +    int storageLayout = 0;
712 +    storageLayout |= DataStorage::dslPosition;
713 +    storageLayout |= DataStorage::dslVelocity;
714 +    storageLayout |= DataStorage::dslForce;
715 +
716 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
717 +
718 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
719 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
720 +      EAMAdapter ea = EAMAdapter( (*i) );
721 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
722 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
723 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
724 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
725 +
726 +      if (da.isDirectional()){
727 +        hasDirectionalAtoms = true;
728 +      }
729 +      if (ma.isMultipole()){
730 +        hasMultipoles = true;
731 +      }
732 +      if (ea.isEAM() || sca.isSuttonChen()){
733 +        hasMetallic = true;
734 +      }
735 +      if ( fca.isFixedCharge() ){
736 +        hasFixedCharge = true;
737 +      }
738 +      if ( fqa.isFluctuatingCharge() ){
739 +        hasFluctuatingCharge = true;
740 +      }
741 +      if ( pa.isPolarizable() ){
742 +        hasPolarizable = true;
743 +      }
744 +    }
745      
746 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
747 +      storageLayout |= DataStorage::dslAmat;
748 +      if(storageLayout & DataStorage::dslVelocity) {
749 +        storageLayout |= DataStorage::dslAngularMomentum;
750 +      }
751 +      if (storageLayout & DataStorage::dslForce) {
752 +        storageLayout |= DataStorage::dslTorque;
753 +      }
754 +    }
755 +    if (hasMultipoles) {
756 +      storageLayout |= DataStorage::dslElectroFrame;
757 +    }
758 +    if (hasFixedCharge || hasFluctuatingCharge) {
759 +      storageLayout |= DataStorage::dslSkippedCharge;
760 +    }
761 +    if (hasMetallic) {
762 +      storageLayout |= DataStorage::dslDensity;
763 +      storageLayout |= DataStorage::dslFunctional;
764 +      storageLayout |= DataStorage::dslFunctionalDerivative;
765 +    }
766 +    if (hasPolarizable) {
767 +      storageLayout |= DataStorage::dslElectricField;
768 +    }
769 +    if (hasFluctuatingCharge){
770 +      storageLayout |= DataStorage::dslFlucQPosition;
771 +      if(storageLayout & DataStorage::dslVelocity) {
772 +        storageLayout |= DataStorage::dslFlucQVelocity;
773 +      }
774 +      if (storageLayout & DataStorage::dslForce) {
775 +        storageLayout |= DataStorage::dslFlucQForce;
776 +      }
777 +    }
778 +    
779 +    // if the user has asked for them, make sure we've got the memory for the
780 +    // objects defined.
781 +
782 +    if (simParams->getOutputParticlePotential()) {
783 +      storageLayout |= DataStorage::dslParticlePot;
784 +    }
785 +
786 +    if (simParams->havePrintHeatFlux()) {
787 +      if (simParams->getPrintHeatFlux()) {
788 +        storageLayout |= DataStorage::dslParticlePot;
789 +      }
790 +    }
791 +
792 +    if (simParams->getOutputElectricField()) {
793 +      storageLayout |= DataStorage::dslElectricField;
794 +    }
795 +    if (simParams->getOutputFluctuatingCharges()) {
796 +      storageLayout |= DataStorage::dslFlucQPosition;
797 +      storageLayout |= DataStorage::dslFlucQVelocity;
798 +      storageLayout |= DataStorage::dslFlucQForce;
799 +    }
800 +
801 +    return storageLayout;
802 +  }
803 +
804    void SimCreator::setGlobalIndex(SimInfo *info) {
805      SimInfo::MoleculeIterator mi;
806      Molecule::AtomIterator ai;
# Line 640 | Line 815 | namespace OpenMD {
815      int beginRigidBodyIndex;
816      int beginCutoffGroupIndex;
817      int nGlobalAtoms = info->getNGlobalAtoms();
818 <
644 <    /**@todo fixme */
645 < #ifndef IS_MPI
818 >    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
819      
820      beginAtomIndex = 0;
648    beginRigidBodyIndex = 0;
649    beginCutoffGroupIndex = 0;
650    
651 #else
652    
653    int nproc;
654    int myNode;
655    
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    
821      //rigidbody's index begins right after atom's
822 <    beginRigidBodyIndex += info->getNGlobalAtoms();
823 <    
824 <    for(mol = info->beginMolecule(mi); mol != NULL;
825 <        mol = info->nextMolecule(mi)) {
822 >    beginRigidBodyIndex = info->getNGlobalAtoms();
823 >    beginCutoffGroupIndex = 0;
824 >
825 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
826        
827 <      //local index(index in DataStorge) of atom is important
828 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
829 <        atom->setGlobalIndex(beginAtomIndex++);
827 > #ifdef IS_MPI      
828 >      if (info->getMolToProc(i) == worldRank) {
829 > #endif        
830 >        // stuff to do if I own this molecule
831 >        mol = info->getMoleculeByGlobalIndex(i);
832 >
833 >        //local index(index in DataStorge) of atom is important
834 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
835 >          atom->setGlobalIndex(beginAtomIndex++);
836 >        }
837 >        
838 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
839 >            rb = mol->nextRigidBody(ri)) {
840 >          rb->setGlobalIndex(beginRigidBodyIndex++);
841 >        }
842 >        
843 >        //local index of cutoff group is trivial, it only depends on
844 >        //the order of travesing
845 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
846 >            cg = mol->nextCutoffGroup(ci)) {
847 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
848 >        }        
849 >        
850 > #ifdef IS_MPI        
851 >      }  else {
852 >
853 >        // stuff to do if I don't own this molecule
854 >        
855 >        int stampId = info->getMoleculeStampId(i);
856 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
857 >
858 >        beginAtomIndex += stamp->getNAtoms();
859 >        beginRigidBodyIndex += stamp->getNRigidBodies();
860 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
861        }
862 <      
863 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
864 <          rb = mol->nextRigidBody(ri)) {
865 <        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 <    
862 > #endif          
863 >
864 >    } //end for(int i=0)  
865 >
866      //fill globalGroupMembership
867      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
868      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
# Line 721 | Line 874 | namespace OpenMD {
874          
875        }      
876      }
877 <    
877 >  
878   #ifdef IS_MPI    
879      // Since the globalGroupMembership has been zero filled and we've only
880      // poked values into the atoms we know, we can do an Allreduce
# Line 729 | Line 882 | namespace OpenMD {
882      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
883      // docs said we could.
884      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
885 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
886 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
885 >    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
886 >                              &tmpGroupMembership[0], nGlobalAtoms,
887 >                              MPI::INT, MPI::SUM);
888      info->setGlobalGroupMembership(tmpGroupMembership);
889   #else
890      info->setGlobalGroupMembership(globalGroupMembership);
891   #endif
892      
893      //fill molMembership
894 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
894 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
895 >                                         info->getNGlobalRigidBodies(), 0);
896      
897 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
897 >    for(mol = info->beginMolecule(mi); mol != NULL;
898 >        mol = info->nextMolecule(mi)) {
899        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
900          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
901        }
902 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
903 +           rb = mol->nextRigidBody(ri)) {
904 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
905 +      }
906      }
907      
908   #ifdef IS_MPI
909 <    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
909 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
910 >                                      info->getNGlobalRigidBodies(), 0);
911 >    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
912 >                              nGlobalAtoms + nGlobalRigidBodies,
913 >                              MPI::INT, MPI::SUM);
914      
751    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
752                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
753    
915      info->setGlobalMolMembership(tmpMolMembership);
916   #else
917      info->setGlobalMolMembership(globalMolMembership);
# Line 760 | Line 921 | namespace OpenMD {
921      // here the molecules are listed by their global indices.
922  
923      std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
924 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
924 >    for (mol = info->beginMolecule(mi); mol != NULL;
925 >         mol = info->nextMolecule(mi)) {
926        nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
927      }
928      
929   #ifdef IS_MPI
930      std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
931 <    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
932 <                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
931 >    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
932 >                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
933   #else
934      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
935   #endif    
# Line 781 | Line 943 | namespace OpenMD {
943      }
944      
945      std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
946 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
946 >    for (mol = info->beginMolecule(mi); mol != NULL;
947 >         mol = info->nextMolecule(mi)) {
948        int myGlobalIndex = mol->getGlobalIndex();
949        int globalIO = startingIOIndexForMol[myGlobalIndex];
950 <      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
951 <           integrableObject = mol->nextIntegrableObject(ioi)) {
952 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
953 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
950 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
951 >           sd = mol->nextIntegrableObject(ioi)) {
952 >        sd->setGlobalIntegrableObjectIndex(globalIO);
953 >        IOIndexToIntegrableObject[globalIO] = sd;
954          globalIO++;
955        }
956      }
957 <    
957 >      
958      info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
959      
960    }
961    
962    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
963 <    Globals* simParams;
801 <    simParams = info->getSimParams();
802 <    
803 <    
963 >
964      DumpReader reader(info, mdFileName);
965      int nframes = reader.getNFrames();
966 <    
966 >
967      if (nframes > 0) {
968        reader.readFrame(nframes - 1);
969      } else {
# Line 814 | Line 974 | namespace OpenMD {
974        painCave.isFatal = 1;
975        simError();
976      }
817    
977      //copy the current snapshot to previous snapshot
978      info->getSnapshotManager()->advance();
979    }

Comparing trunk/src/brains/SimCreator.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1810 by gezelter, Thu Nov 8 14:23:43 2012 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines