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

Comparing:
trunk/src/brains/SimCreator.cpp (file contents), Revision 1313 by gezelter, Wed Oct 22 20:01:49 2008 UTC vs.
branches/development/src/brains/SimCreator.cpp (file contents), Revision 1826 by gezelter, Wed Jan 9 19:41:48 2013 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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]  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   #include <exception>
# Line 55 | Line 55
55   #include "brains/SimCreator.hpp"
56   #include "brains/SimSnapshotManager.hpp"
57   #include "io/DumpReader.hpp"
58 < #include "UseTheForce/ForceFieldFactory.hpp"
58 > #include "brains/ForceField.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61   #include "math/SeqRandNumGen.hpp"
# Line 75 | Line 75
75   #include "antlr/NoViableAltForCharException.hpp"
76   #include "antlr/NoViableAltException.hpp"
77  
78 + #include "types/DirectionalAdapter.hpp"
79 + #include "types/MultipoleAdapter.hpp"
80 + #include "types/EAMAdapter.hpp"
81 + #include "types/SuttonChenAdapter.hpp"
82 + #include "types/PolarizableAdapter.hpp"
83 + #include "types/FixedChargeAdapter.hpp"
84 + #include "types/FluctuatingChargeAdapter.hpp"
85 +
86   #ifdef IS_MPI
87 + #include "mpi.h"
88   #include "math/ParallelRandNumGen.hpp"
89   #endif
90  
91 < namespace oopse {
91 > namespace OpenMD {
92    
93 <  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
93 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
94      Globals* simParams = NULL;
95      try {
96  
# Line 92 | Line 101 | namespace oopse {
101        const int masterNode = 0;
102        int commStatus;
103        if (worldRank == masterNode) {
104 < #endif
105 <                
104 >        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
105 > #endif                
106          SimplePreprocessor preprocessor;
107          preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
108                  
# Line 106 | Line 115 | namespace oopse {
115              
116                  
117        } else {
118 +
119 +        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
120 +
121          //get stream size
122          commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
123  
# Line 214 | Line 226 | namespace oopse {
226        painCave.isFatal = 1;
227        simError();        
228      }
229 <    catch (OOPSEException& e) {
229 >    catch (OpenMDException& e) {
230        sprintf(painCave.errMsg,
231                "%s\n",
232                e.getMessage().c_str());
# Line 229 | Line 241 | namespace oopse {
241        simError();
242      }
243  
244 +    simParams->setMDfileVersion(mdFileVersion);
245      return simParams;
246    }
247    
248    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
249                                    bool loadInitCoords) {
250 <
250 >    
251      const int bufferSize = 65535;
252      char buffer[bufferSize];
253      int lineNo = 0;
254      std::string mdRawData;
255      int metaDataBlockStart = -1;
256      int metaDataBlockEnd = -1;
257 <    int i;
258 <    int mdOffset;
257 >    int i, j;
258 >    streamoff mdOffset;
259 >    int mdFileVersion;
260  
261 +    // Create a string for embedding the version information in the MetaData
262 +    std::string version;
263 +    version.assign("## Last run using OpenMD Version: ");
264 +    version.append(OPENMD_VERSION_MAJOR);
265 +    version.append(".");
266 +    version.append(OPENMD_VERSION_MINOR);
267 +
268 +    std::string svnrev;
269 +    //convert a macro from compiler to a string in c++
270 +    STR_DEFINE(svnrev, SVN_REV );
271 +    version.append(" Revision: ");
272 +    // If there's no SVN revision, just call this the RELEASE revision.
273 +    if (!svnrev.empty()) {
274 +      version.append(svnrev);
275 +    } else {
276 +      version.append("RELEASE");
277 +    }
278 +  
279   #ifdef IS_MPI            
280      const int masterNode = 0;
281      if (worldRank == masterNode) {
282   #endif
283  
284 <      std::ifstream mdFile_(mdFileName.c_str());
284 >      std::ifstream mdFile_;
285 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
286        
287        if (mdFile_.fail()) {
288          sprintf(painCave.errMsg,
# Line 262 | Line 295 | namespace oopse {
295        mdFile_.getline(buffer, bufferSize);
296        ++lineNo;
297        std::string line = trimLeftCopy(buffer);
298 <      i = CaseInsensitiveFind(line, "<OOPSE");
298 >      i = CaseInsensitiveFind(line, "<OpenMD");
299        if (static_cast<size_t>(i) == string::npos) {
300 +        // try the older file strings to see if that works:
301 +        i = CaseInsensitiveFind(line, "<OOPSE");
302 +      }
303 +      
304 +      if (static_cast<size_t>(i) == string::npos) {
305 +        // still no luck!
306          sprintf(painCave.errMsg,
307 <                "SimCreator: File: %s is not an OOPSE file!\n",
307 >                "SimCreator: File: %s is not a valid OpenMD file!\n",
308                  mdFileName.c_str());
309          painCave.isFatal = 1;
310          simError();
311        }
312 +      
313 +      // found the correct opening string, now try to get the file
314 +      // format version number.
315  
316 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
317 +      std::string fileType = tokenizer.nextToken();
318 +      toUpper(fileType);
319 +
320 +      mdFileVersion = 0;
321 +
322 +      if (fileType == "OPENMD") {
323 +        while (tokenizer.hasMoreTokens()) {
324 +          std::string token(tokenizer.nextToken());
325 +          toUpper(token);
326 +          if (token == "VERSION") {
327 +            mdFileVersion = tokenizer.nextTokenAsInt();
328 +            break;
329 +          }
330 +        }
331 +      }
332 +            
333        //scan through the input stream and find MetaData tag        
334        while(mdFile_.getline(buffer, bufferSize)) {
335          ++lineNo;
# Line 311 | Line 370 | namespace oopse {
370  
371        mdRawData.clear();
372  
373 +      bool foundVersion = false;
374 +
375        for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
376          mdFile_.getline(buffer, bufferSize);
377 <        mdRawData += buffer;
377 >        std::string line = trimLeftCopy(buffer);
378 >        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version");
379 >        if (static_cast<size_t>(j) != string::npos) {
380 >          foundVersion = true;
381 >          mdRawData += version;
382 >        } else {
383 >          mdRawData += buffer;
384 >        }
385          mdRawData += "\n";
386        }
387 <
387 >      
388 >      if (!foundVersion) mdRawData += version + "\n";
389 >      
390        mdFile_.close();
391  
392   #ifdef IS_MPI
# Line 326 | Line 396 | namespace oopse {
396      std::stringstream rawMetaDataStream(mdRawData);
397  
398      //parse meta-data file
399 <    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
399 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
400 >                                   metaDataBlockStart + 1);
401      
402      //create the force field
403 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
403 >    ForceField * ff = new ForceField(simParams->getForceField());
404  
405      if (ff == NULL) {
406        sprintf(painCave.errMsg,
# Line 363 | Line 434 | namespace oopse {
434      }
435      
436      ff->parse(forcefieldFileName);
366    ff->setFortranForceOptions();
437      //create SimInfo
438      SimInfo * info = new SimInfo(ff, simParams);
439  
# Line 381 | Line 451 | namespace oopse {
451      //create the molecules
452      createMolecules(info);
453      
454 <    
454 >    //find the storage layout
455 >
456 >    int storageLayout = computeStorageLayout(info);
457 >
458      //allocate memory for DataStorage(circular reference, need to
459      //break it)
460 <    info->setSnapshotManager(new SimSnapshotManager(info));
460 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
461      
462      //set the global index of atoms, rigidbodies and cutoffgroups
463      //(only need to be set once, the global index will never change
# Line 407 | Line 480 | namespace oopse {
480      
481      if (loadInitCoords)
482        loadCoordinates(info, mdFileName);    
410    
483      return info;
484    }
485    
# Line 454 | Line 526 | namespace oopse {
526      int nTarget;
527      int done;
528      int i;
457    int j;
529      int loops;
530      int which_proc;
531      int nProcessors;
# Line 462 | Line 533 | namespace oopse {
533      int nGlobalMols = info->getNGlobalMolecules();
534      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
535      
536 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
536 >    nProcessors = MPI::COMM_WORLD.Get_size();
537      
538      if (nProcessors > nGlobalMols) {
539        sprintf(painCave.errMsg,
# Line 471 | Line 542 | namespace oopse {
542                "\tthe number of molecules.  This will not result in a \n"
543                "\tusable division of atoms for force decomposition.\n"
544                "\tEither try a smaller number of processors, or run the\n"
545 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
545 >              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
546        
547        painCave.isFatal = 1;
548        simError();
# Line 500 | Line 571 | namespace oopse {
571        nTarget = (int)(precast + 0.5);
572        
573        for(i = 0; i < nGlobalMols; i++) {
574 +
575          done = 0;
576          loops = 0;
577          
# Line 524 | Line 596 | namespace oopse {
596            // and be done with it.
597            
598            if (loops > 100) {
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",
599 >
600 >            sprintf(painCave.errMsg,
601 >                    "There have been 100 attempts to assign molecule %d to an\n"
602 >                    "\tunderworked processor, but there's no good place to\n"
603 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
604                      i, which_proc);
605 <            
605 >          
606              painCave.isFatal = 0;
607 +            painCave.severity = OPENMD_INFO;
608              simError();
609              
610              molToProcMap[i] = which_proc;
# Line 575 | Line 649 | namespace oopse {
649        }
650        
651        delete myRandom;
652 <      
652 >
653        // Spray out this nonsense to all other processors:
654 <      
581 <      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:
658 <      
659 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
658 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
659 >
660      }
661      
662      info->setMolToProcMap(molToProcMap);
# Line 606 | Line 679 | namespace oopse {
679   #endif
680          
681          stampId = info->getMoleculeStampId(i);
682 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
683 <                                                   stampId, i, info->getLocalIndexManager());
682 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
683 >                                                   info->getMoleculeStamp(stampId),
684 >                                                   stampId, i,
685 >                                                   info->getLocalIndexManager());
686          
687          info->addMolecule(mol);
688          
# Line 619 | Line 694 | namespace oopse {
694        
695      } //end for(int i=0)  
696    }
697 +    
698 +  int SimCreator::computeStorageLayout(SimInfo* info) {
699 +
700 +    Globals* simParams = info->getSimParams();
701 +    int nRigidBodies = info->getNGlobalRigidBodies();
702 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
703 +    set<AtomType*>::iterator i;
704 +    bool hasDirectionalAtoms = false;
705 +    bool hasFixedCharge = false;
706 +    bool hasDipoles = false;    
707 +    bool hasQuadrupoles = 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.isDipole()){
730 +        hasDipoles = true;
731 +      }
732 +      if (ma.isQuadrupole()){
733 +        hasQuadrupoles = true;
734 +      }
735 +      if (ea.isEAM() || sca.isSuttonChen()){
736 +        hasMetallic = true;
737 +      }
738 +      if ( fca.isFixedCharge() ){
739 +        hasFixedCharge = true;
740 +      }
741 +      if ( fqa.isFluctuatingCharge() ){
742 +        hasFluctuatingCharge = true;
743 +      }
744 +      if ( pa.isPolarizable() ){
745 +        hasPolarizable = true;
746 +      }
747 +    }
748 +    
749 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
750 +      storageLayout |= DataStorage::dslAmat;
751 +      if(storageLayout & DataStorage::dslVelocity) {
752 +        storageLayout |= DataStorage::dslAngularMomentum;
753 +      }
754 +      if (storageLayout & DataStorage::dslForce) {
755 +        storageLayout |= DataStorage::dslTorque;
756 +      }
757 +    }
758 +    if (hasDipoles) {
759 +      storageLayout |= DataStorage::dslDipole;
760 +    }
761 +    if (hasQuadrupoles) {
762 +      storageLayout |= DataStorage::dslQuadrupole;
763 +    }
764 +    if (hasFixedCharge || hasFluctuatingCharge) {
765 +      storageLayout |= DataStorage::dslSkippedCharge;
766 +    }
767 +    if (hasMetallic) {
768 +      storageLayout |= DataStorage::dslDensity;
769 +      storageLayout |= DataStorage::dslFunctional;
770 +      storageLayout |= DataStorage::dslFunctionalDerivative;
771 +    }
772 +    if (hasPolarizable) {
773 +      storageLayout |= DataStorage::dslElectricField;
774 +    }
775 +    if (hasFluctuatingCharge){
776 +      storageLayout |= DataStorage::dslFlucQPosition;
777 +      if(storageLayout & DataStorage::dslVelocity) {
778 +        storageLayout |= DataStorage::dslFlucQVelocity;
779 +      }
780 +      if (storageLayout & DataStorage::dslForce) {
781 +        storageLayout |= DataStorage::dslFlucQForce;
782 +      }
783 +    }
784      
785 +    // if the user has asked for them, make sure we've got the memory for the
786 +    // objects defined.
787 +
788 +    if (simParams->getOutputParticlePotential()) {
789 +      storageLayout |= DataStorage::dslParticlePot;
790 +    }
791 +
792 +    if (simParams->havePrintHeatFlux()) {
793 +      if (simParams->getPrintHeatFlux()) {
794 +        storageLayout |= DataStorage::dslParticlePot;
795 +      }
796 +    }
797 +
798 +    if (simParams->getOutputElectricField()) {
799 +      storageLayout |= DataStorage::dslElectricField;
800 +    }
801 +
802 +    if (simParams->getOutputFluctuatingCharges()) {
803 +      storageLayout |= DataStorage::dslFlucQPosition;
804 +      storageLayout |= DataStorage::dslFlucQVelocity;
805 +      storageLayout |= DataStorage::dslFlucQForce;
806 +    }
807 +
808 +    return storageLayout;
809 +  }
810 +
811    void SimCreator::setGlobalIndex(SimInfo *info) {
812      SimInfo::MoleculeIterator mi;
813      Molecule::AtomIterator ai;
# Line 634 | Line 822 | namespace oopse {
822      int beginRigidBodyIndex;
823      int beginCutoffGroupIndex;
824      int nGlobalAtoms = info->getNGlobalAtoms();
825 <
638 <    /**@todo fixme */
639 < #ifndef IS_MPI
825 >    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
826      
827      beginAtomIndex = 0;
642    beginRigidBodyIndex = 0;
643    beginCutoffGroupIndex = 0;
644    
645 #else
646    
647    int nproc;
648    int myNode;
649    
650    myNode = worldRank;
651    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
652    
653    std::vector < int > tmpAtomsInProc(nproc, 0);
654    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
655    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
656    std::vector < int > NumAtomsInProc(nproc, 0);
657    std::vector < int > NumRigidBodiesInProc(nproc, 0);
658    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
659    
660    tmpAtomsInProc[myNode] = info->getNAtoms();
661    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
662    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
663    
664    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
665    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
666                  MPI_SUM, MPI_COMM_WORLD);
667    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
668                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
669    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
670                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
671    
672    beginAtomIndex = 0;
673    beginRigidBodyIndex = 0;
674    beginCutoffGroupIndex = 0;
675    
676    for(int i = 0; i < myNode; i++) {
677      beginAtomIndex += NumAtomsInProc[i];
678      beginRigidBodyIndex += NumRigidBodiesInProc[i];
679      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
680    }
681    
682 #endif
683    
828      //rigidbody's index begins right after atom's
829 <    beginRigidBodyIndex += info->getNGlobalAtoms();
830 <    
831 <    for(mol = info->beginMolecule(mi); mol != NULL;
832 <        mol = info->nextMolecule(mi)) {
829 >    beginRigidBodyIndex = info->getNGlobalAtoms();
830 >    beginCutoffGroupIndex = 0;
831 >
832 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
833        
834 <      //local index(index in DataStorge) of atom is important
835 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
836 <        atom->setGlobalIndex(beginAtomIndex++);
834 > #ifdef IS_MPI      
835 >      if (info->getMolToProc(i) == worldRank) {
836 > #endif        
837 >        // stuff to do if I own this molecule
838 >        mol = info->getMoleculeByGlobalIndex(i);
839 >
840 >        //local index(index in DataStorge) of atom is important
841 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
842 >          atom->setGlobalIndex(beginAtomIndex++);
843 >        }
844 >        
845 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
846 >            rb = mol->nextRigidBody(ri)) {
847 >          rb->setGlobalIndex(beginRigidBodyIndex++);
848 >        }
849 >        
850 >        //local index of cutoff group is trivial, it only depends on
851 >        //the order of travesing
852 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
853 >            cg = mol->nextCutoffGroup(ci)) {
854 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
855 >        }        
856 >        
857 > #ifdef IS_MPI        
858 >      }  else {
859 >
860 >        // stuff to do if I don't own this molecule
861 >        
862 >        int stampId = info->getMoleculeStampId(i);
863 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
864 >
865 >        beginAtomIndex += stamp->getNAtoms();
866 >        beginRigidBodyIndex += stamp->getNRigidBodies();
867 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
868        }
869 <      
870 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
871 <          rb = mol->nextRigidBody(ri)) {
872 <        rb->setGlobalIndex(beginRigidBodyIndex++);
698 <      }
699 <      
700 <      //local index of cutoff group is trivial, it only depends on the order of travesing
701 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
702 <          cg = mol->nextCutoffGroup(ci)) {
703 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
704 <      }
705 <    }
706 <    
869 > #endif          
870 >
871 >    } //end for(int i=0)  
872 >
873      //fill globalGroupMembership
874      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
875      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
# Line 715 | Line 881 | namespace oopse {
881          
882        }      
883      }
884 <    
884 >  
885   #ifdef IS_MPI    
886      // Since the globalGroupMembership has been zero filled and we've only
887      // poked values into the atoms we know, we can do an Allreduce
# Line 723 | Line 889 | namespace oopse {
889      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
890      // docs said we could.
891      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
892 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
893 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
892 >    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
893 >                              &tmpGroupMembership[0], nGlobalAtoms,
894 >                              MPI::INT, MPI::SUM);
895      info->setGlobalGroupMembership(tmpGroupMembership);
896   #else
897      info->setGlobalGroupMembership(globalGroupMembership);
898   #endif
899      
900      //fill molMembership
901 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
901 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
902 >                                         info->getNGlobalRigidBodies(), 0);
903      
904 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
904 >    for(mol = info->beginMolecule(mi); mol != NULL;
905 >        mol = info->nextMolecule(mi)) {
906        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
907          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
908        }
909 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
910 +           rb = mol->nextRigidBody(ri)) {
911 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
912 +      }
913      }
914      
915   #ifdef IS_MPI
916 <    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
916 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
917 >                                      info->getNGlobalRigidBodies(), 0);
918 >    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
919 >                              nGlobalAtoms + nGlobalRigidBodies,
920 >                              MPI::INT, MPI::SUM);
921      
745    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
746                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
747    
922      info->setGlobalMolMembership(tmpMolMembership);
923   #else
924      info->setGlobalMolMembership(globalMolMembership);
# Line 754 | Line 928 | namespace oopse {
928      // here the molecules are listed by their global indices.
929  
930      std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
931 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
931 >    for (mol = info->beginMolecule(mi); mol != NULL;
932 >         mol = info->nextMolecule(mi)) {
933        nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
934      }
935      
936   #ifdef IS_MPI
937      std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
938 <    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
939 <                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
938 >    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
939 >                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
940   #else
941      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
942   #endif    
# Line 775 | Line 950 | namespace oopse {
950      }
951      
952      std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
953 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
953 >    for (mol = info->beginMolecule(mi); mol != NULL;
954 >         mol = info->nextMolecule(mi)) {
955        int myGlobalIndex = mol->getGlobalIndex();
956        int globalIO = startingIOIndexForMol[myGlobalIndex];
957 <      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
958 <           integrableObject = mol->nextIntegrableObject(ioi)) {
959 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
960 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
957 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
958 >           sd = mol->nextIntegrableObject(ioi)) {
959 >        sd->setGlobalIntegrableObjectIndex(globalIO);
960 >        IOIndexToIntegrableObject[globalIO] = sd;
961          globalIO++;
962        }
963      }
964 <    
964 >      
965      info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
966      
967    }
968    
969    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
794    Globals* simParams;
795    simParams = info->getSimParams();
970      
797    
971      DumpReader reader(info, mdFileName);
972      int nframes = reader.getNFrames();
973      
# Line 808 | Line 981 | namespace oopse {
981        painCave.isFatal = 1;
982        simError();
983      }
811    
984      //copy the current snapshot to previous snapshot
985      info->getSnapshotManager()->advance();
986    }
987    
988 < } //end namespace oopse
988 > } //end namespace OpenMD
989  
990  

Comparing:
trunk/src/brains/SimCreator.cpp (property svn:keywords), Revision 1313 by gezelter, Wed Oct 22 20:01:49 2008 UTC vs.
branches/development/src/brains/SimCreator.cpp (property svn:keywords), Revision 1826 by gezelter, Wed Jan 9 19:41:48 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines