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

Comparing branches/development/src/brains/SimCreator.cpp (file contents):
Revision 1597 by gezelter, Tue Jul 26 15:49:24 2011 UTC vs.
Revision 1711 by gezelter, Sat May 19 02:58:35 2012 UTC

# 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 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 92 | Line 102 | namespace OpenMD {
102        const int masterNode = 0;
103        int commStatus;
104        if (worldRank == masterNode) {
105 < #endif
106 <                
105 >        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
106 > #endif                
107          SimplePreprocessor preprocessor;
108          preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
109                  
# Line 106 | Line 116 | namespace OpenMD {
116              
117                  
118        } else {
119 +
120 +        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
121 +
122          //get stream size
123          commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
124  
# Line 229 | Line 242 | namespace OpenMD {
242        simError();
243      }
244  
245 +    simParams->setMDfileVersion(mdFileVersion);
246      return simParams;
247    }
248    
# Line 243 | Line 257 | namespace OpenMD {
257      int metaDataBlockEnd = -1;
258      int i;
259      int mdOffset;
260 +    int mdFileVersion;
261  
262   #ifdef IS_MPI            
263      const int masterNode = 0;
# Line 276 | Line 291 | namespace OpenMD {
291          painCave.isFatal = 1;
292          simError();
293        }
294 +      
295 +      // found the correct opening string, now try to get the file
296 +      // format version number.
297  
298 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
299 +      std::string fileType = tokenizer.nextToken();
300 +      toUpper(fileType);
301 +
302 +      mdFileVersion = 0;
303 +
304 +      if (fileType == "OPENMD") {
305 +        while (tokenizer.hasMoreTokens()) {
306 +          std::string token(tokenizer.nextToken());
307 +          toUpper(token);
308 +          if (token == "VERSION") {
309 +            mdFileVersion = tokenizer.nextTokenAsInt();
310 +            break;
311 +          }
312 +        }
313 +      }
314 +            
315        //scan through the input stream and find MetaData tag        
316        while(mdFile_.getline(buffer, bufferSize)) {
317          ++lineNo;
# Line 332 | Line 367 | namespace OpenMD {
367      std::stringstream rawMetaDataStream(mdRawData);
368  
369      //parse meta-data file
370 <    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
370 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
371 >                                   metaDataBlockStart + 1);
372      
373      //create the force field
374      ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
# Line 386 | Line 422 | namespace OpenMD {
422      //create the molecules
423      createMolecules(info);
424      
425 +    //find the storage layout
426 +
427 +    int storageLayout = computeStorageLayout(info);
428 +
429 +    cerr << "computed Storage Layout = " << storageLayout << "\n";
430 +
431      //allocate memory for DataStorage(circular reference, need to
432      //break it)
433 <    info->setSnapshotManager(new SimSnapshotManager(info));
433 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
434      
435      //set the global index of atoms, rigidbodies and cutoffgroups
436      //(only need to be set once, the global index will never change
# Line 589 | Line 631 | namespace OpenMD {
631        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
632      }
633      
592 –    cerr << "molToProcMap:\n";
593 –    for (int i = 0; i < molToProcMap.size(); i++) {
594 –      cerr << "m = " << i << " mtpr[m] = " << molToProcMap[i] <<"\n";
595 –    }
596 –
634      info->setMolToProcMap(molToProcMap);
635      sprintf(checkPointMsg,
636              "Successfully divided the molecules among the processors.\n");
# Line 614 | Line 651 | namespace OpenMD {
651   #endif
652          
653          stampId = info->getMoleculeStampId(i);
654 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
655 <                                                   stampId, i, info->getLocalIndexManager());
654 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
655 >                                                   info->getMoleculeStamp(stampId),
656 >                                                   stampId, i,
657 >                                                   info->getLocalIndexManager());
658          
659          info->addMolecule(mol);
660          
# Line 627 | Line 666 | namespace OpenMD {
666        
667      } //end for(int i=0)  
668    }
669 +    
670 +  int SimCreator::computeStorageLayout(SimInfo* info) {
671 +
672 +    Globals* simParams = info->getSimParams();
673 +    int nRigidBodies = info->getNGlobalRigidBodies();
674 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
675 +    set<AtomType*>::iterator i;
676 +    bool hasDirectionalAtoms = false;
677 +    bool hasFixedCharge = false;
678 +    bool hasMultipoles = false;    
679 +    bool hasPolarizable = false;    
680 +    bool hasFluctuatingCharge = false;    
681 +    bool hasMetallic = false;
682 +    int storageLayout = 0;
683 +    storageLayout |= DataStorage::dslPosition;
684 +    storageLayout |= DataStorage::dslVelocity;
685 +    storageLayout |= DataStorage::dslForce;
686 +
687 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
688 +
689 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
690 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
691 +      EAMAdapter ea = EAMAdapter( (*i) );
692 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
693 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
694 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
695 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
696 +
697 +      if (da.isDirectional()){
698 +        hasDirectionalAtoms = true;
699 +      }
700 +      if (ma.isMultipole()){
701 +        hasMultipoles = true;
702 +      }
703 +      if (ea.isEAM() || sca.isSuttonChen()){
704 +        hasMetallic = true;
705 +      }
706 +      if ( fca.isFixedCharge() ){
707 +        hasFixedCharge = true;
708 +      }
709 +      if ( fqa.isFluctuatingCharge() ){
710 +        hasFluctuatingCharge = true;
711 +      }
712 +      if ( pa.isPolarizable() ){
713 +        hasPolarizable = true;
714 +      }
715 +    }
716      
717 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
718 +      storageLayout |= DataStorage::dslAmat;
719 +      if(storageLayout & DataStorage::dslVelocity) {
720 +        storageLayout |= DataStorage::dslAngularMomentum;
721 +      }
722 +      if (storageLayout & DataStorage::dslForce) {
723 +        storageLayout |= DataStorage::dslTorque;
724 +      }
725 +    }
726 +    if (hasMultipoles) {
727 +      storageLayout |= DataStorage::dslElectroFrame;
728 +    }
729 +    if (hasFixedCharge || hasFluctuatingCharge) {
730 +      storageLayout |= DataStorage::dslSkippedCharge;
731 +    }
732 +    if (hasMetallic) {
733 +      storageLayout |= DataStorage::dslDensity;
734 +      storageLayout |= DataStorage::dslFunctional;
735 +      storageLayout |= DataStorage::dslFunctionalDerivative;
736 +    }
737 +    if (hasPolarizable) {
738 +      storageLayout |= DataStorage::dslElectricField;
739 +    }
740 +    if (hasFluctuatingCharge){
741 +      storageLayout |= DataStorage::dslFlucQPosition;
742 +      if(storageLayout & DataStorage::dslVelocity) {
743 +        storageLayout |= DataStorage::dslFlucQVelocity;
744 +      }
745 +      if (storageLayout & DataStorage::dslForce) {
746 +        storageLayout |= DataStorage::dslFlucQForce;
747 +      }
748 +    }
749 +
750 +    if (simParams->getOutputParticlePotential()) {
751 +      storageLayout |= DataStorage::dslParticlePot;
752 +    }
753 +
754 +    return storageLayout;
755 +  }
756 +
757    void SimCreator::setGlobalIndex(SimInfo *info) {
758      SimInfo::MoleculeIterator mi;
759      Molecule::AtomIterator ai;
# Line 642 | Line 768 | namespace OpenMD {
768      int beginRigidBodyIndex;
769      int beginCutoffGroupIndex;
770      int nGlobalAtoms = info->getNGlobalAtoms();
645 –
646 –    /**@todo fixme */
647 – #ifndef IS_MPI
771      
772      beginAtomIndex = 0;
773      beginRigidBodyIndex = 0;
774      beginCutoffGroupIndex = 0;
775 <    
776 < #else
654 <    
655 <    int nproc;
656 <    int myNode;
657 <    
658 <    myNode = worldRank;
659 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
660 <    
661 <    std::vector < int > tmpAtomsInProc(nproc, 0);
662 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
663 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
664 <    std::vector < int > NumAtomsInProc(nproc, 0);
665 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
666 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
667 <    
668 <    tmpAtomsInProc[myNode] = info->getNAtoms();
669 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
670 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
671 <    
672 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
673 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
674 <                  MPI_SUM, MPI_COMM_WORLD);
675 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
676 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
677 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
678 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
679 <    
680 <    beginAtomIndex = 0;
681 <    beginRigidBodyIndex = 0;
682 <    beginCutoffGroupIndex = 0;
683 <    
684 <    for(int i = 0; i < myNode; i++) {
685 <      beginAtomIndex += NumAtomsInProc[i];
686 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
687 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
688 <    }
689 <    
690 < #endif
691 <    
692 <    //rigidbody's index begins right after atom's
693 <    beginRigidBodyIndex += info->getNGlobalAtoms();
694 <    
695 <    for(mol = info->beginMolecule(mi); mol != NULL;
696 <        mol = info->nextMolecule(mi)) {
775 >
776 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
777        
778 <      //local index(index in DataStorge) of atom is important
779 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
780 <        atom->setGlobalIndex(beginAtomIndex++);
778 > #ifdef IS_MPI      
779 >      if (info->getMolToProc(i) == worldRank) {
780 > #endif        
781 >        // stuff to do if I own this molecule
782 >        mol = info->getMoleculeByGlobalIndex(i);
783 >
784 >        //local index(index in DataStorge) of atom is important
785 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
786 >          atom->setGlobalIndex(beginAtomIndex++);
787 >        }
788 >        
789 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
790 >            rb = mol->nextRigidBody(ri)) {
791 >          rb->setGlobalIndex(beginRigidBodyIndex++);
792 >        }
793 >        
794 >        //local index of cutoff group is trivial, it only depends on
795 >        //the order of travesing
796 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
797 >            cg = mol->nextCutoffGroup(ci)) {
798 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
799 >        }        
800 >        
801 > #ifdef IS_MPI        
802 >      }  else {
803 >
804 >        // stuff to do if I don't own this molecule
805 >        
806 >        int stampId = info->getMoleculeStampId(i);
807 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
808 >
809 >        beginAtomIndex += stamp->getNAtoms();
810 >        beginRigidBodyIndex += stamp->getNRigidBodies();
811 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
812        }
813 <      
814 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
815 <          rb = mol->nextRigidBody(ri)) {
816 <        rb->setGlobalIndex(beginRigidBodyIndex++);
706 <      }
707 <      
708 <      //local index of cutoff group is trivial, it only depends on the order of travesing
709 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
710 <          cg = mol->nextCutoffGroup(ci)) {
711 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
712 <      }
713 <    }
714 <    
813 > #endif          
814 >
815 >    } //end for(int i=0)  
816 >
817      //fill globalGroupMembership
818 <    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), -1);
818 >    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
819      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
820        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
821          
# Line 732 | Line 834 | namespace OpenMD {
834      // docs said we could.
835      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
836      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
837 <                  MPI_INT, MPI_MAX, MPI_COMM_WORLD);
837 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
838      info->setGlobalGroupMembership(tmpGroupMembership);
737 –
738 –    cerr << "ggm:\n";
739 –    for (int i = 0; i < tmpGroupMembership.size(); i++)
740 –      cerr << "i = " << i << "\t ggm(i) = " << tmpGroupMembership[i] << "\n";
741 –
839   #else
840      info->setGlobalGroupMembership(globalGroupMembership);
841   #endif
# Line 798 | Line 895 | namespace OpenMD {
895          globalIO++;
896        }
897      }
801 –    cerr << "ioi2io:\n";
802 –    for (int i = 0; i < IOIndexToIntegrableObject.size(); i++) {
803 –      if (IOIndexToIntegrableObject[i] != NULL) {
804 –        cerr << "i = " << i << "globalIOindex = " << IOIndexToIntegrableObject[i]->getGlobalIntegrableObjectIndex() << "\n";
805 –      }
806 –    }
898        
899      info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
900      

Diff Legend

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