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 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1787 by gezelter, Wed Aug 29 18:13:11 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 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 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    
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;
# Line 242 | Line 256 | namespace OpenMD {
256      int metaDataBlockStart = -1;
257      int metaDataBlockEnd = -1;
258      int i;
259 <    int mdOffset;
259 >    streamoff 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());
374 >    ForceField * ff = new ForceField(simParams->getForceField());
375  
376      if (ff == NULL) {
377        sprintf(painCave.errMsg,
# Line 369 | Line 405 | namespace OpenMD {
405      }
406      
407      ff->parse(forcefieldFileName);
372 –    ff->setFortranForceOptions();
408      //create SimInfo
409      SimInfo * info = new SimInfo(ff, simParams);
410  
# Line 387 | Line 422 | namespace OpenMD {
422      //create the molecules
423      createMolecules(info);
424      
425 <    
425 >    //find the storage layout
426 >
427 >    int storageLayout = computeStorageLayout(info);
428 >
429      //allocate memory for DataStorage(circular reference, need to
430      //break it)
431 <    info->setSnapshotManager(new SimSnapshotManager(info));
431 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
432      
433      //set the global index of atoms, rigidbodies and cutoffgroups
434      //(only need to be set once, the global index will never change
# Line 413 | Line 451 | namespace OpenMD {
451      
452      if (loadInitCoords)
453        loadCoordinates(info, mdFileName);    
416 –    
454      return info;
455    }
456    
# Line 612 | Line 649 | namespace OpenMD {
649   #endif
650          
651          stampId = info->getMoleculeStampId(i);
652 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
653 <                                                   stampId, i, info->getLocalIndexManager());
652 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
653 >                                                   info->getMoleculeStamp(stampId),
654 >                                                   stampId, i,
655 >                                                   info->getLocalIndexManager());
656          
657          info->addMolecule(mol);
658          
# Line 626 | Line 665 | namespace OpenMD {
665      } //end for(int i=0)  
666    }
667      
668 +  int SimCreator::computeStorageLayout(SimInfo* info) {
669 +
670 +    Globals* simParams = info->getSimParams();
671 +    int nRigidBodies = info->getNGlobalRigidBodies();
672 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
673 +    set<AtomType*>::iterator i;
674 +    bool hasDirectionalAtoms = false;
675 +    bool hasFixedCharge = false;
676 +    bool hasMultipoles = false;    
677 +    bool hasPolarizable = false;    
678 +    bool hasFluctuatingCharge = false;    
679 +    bool hasMetallic = false;
680 +    int storageLayout = 0;
681 +    storageLayout |= DataStorage::dslPosition;
682 +    storageLayout |= DataStorage::dslVelocity;
683 +    storageLayout |= DataStorage::dslForce;
684 +
685 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
686 +
687 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
688 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
689 +      EAMAdapter ea = EAMAdapter( (*i) );
690 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
691 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
692 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
693 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
694 +
695 +      if (da.isDirectional()){
696 +        hasDirectionalAtoms = true;
697 +      }
698 +      if (ma.isMultipole()){
699 +        hasMultipoles = true;
700 +      }
701 +      if (ea.isEAM() || sca.isSuttonChen()){
702 +        hasMetallic = true;
703 +      }
704 +      if ( fca.isFixedCharge() ){
705 +        hasFixedCharge = true;
706 +      }
707 +      if ( fqa.isFluctuatingCharge() ){
708 +        hasFluctuatingCharge = true;
709 +      }
710 +      if ( pa.isPolarizable() ){
711 +        hasPolarizable = true;
712 +      }
713 +    }
714 +    
715 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
716 +      storageLayout |= DataStorage::dslAmat;
717 +      if(storageLayout & DataStorage::dslVelocity) {
718 +        storageLayout |= DataStorage::dslAngularMomentum;
719 +      }
720 +      if (storageLayout & DataStorage::dslForce) {
721 +        storageLayout |= DataStorage::dslTorque;
722 +      }
723 +    }
724 +    if (hasFixedCharge || hasFluctuatingCharge) {
725 +      storageLayout |= DataStorage::dslSkippedCharge;
726 +    }
727 +    if (hasMetallic) {
728 +      storageLayout |= DataStorage::dslDensity;
729 +      storageLayout |= DataStorage::dslFunctional;
730 +      storageLayout |= DataStorage::dslFunctionalDerivative;
731 +    }
732 +    if (hasPolarizable) {
733 +      storageLayout |= DataStorage::dslElectricField;
734 +    }
735 +    if (hasFluctuatingCharge){
736 +      storageLayout |= DataStorage::dslFlucQPosition;
737 +      if(storageLayout & DataStorage::dslVelocity) {
738 +        storageLayout |= DataStorage::dslFlucQVelocity;
739 +      }
740 +      if (storageLayout & DataStorage::dslForce) {
741 +        storageLayout |= DataStorage::dslFlucQForce;
742 +      }
743 +    }
744 +    
745 +    // if the user has asked for them, make sure we've got the memory for the
746 +    // objects defined.
747 +
748 +    if (simParams->getOutputParticlePotential()) {
749 +      storageLayout |= DataStorage::dslParticlePot;
750 +    }
751 +
752 +    if (simParams->havePrintHeatFlux()) {
753 +      if (simParams->getPrintHeatFlux()) {
754 +        storageLayout |= DataStorage::dslParticlePot;
755 +      }
756 +    }
757 +
758 +    if (simParams->getOutputElectricField()) {
759 +      storageLayout |= DataStorage::dslElectricField;
760 +    }
761 +
762 +    if (simParams->getOutputFluctuatingCharges()) {
763 +      storageLayout |= DataStorage::dslFlucQPosition;
764 +      storageLayout |= DataStorage::dslFlucQVelocity;
765 +      storageLayout |= DataStorage::dslFlucQForce;
766 +    }
767 +
768 +    return storageLayout;
769 +  }
770 +
771    void SimCreator::setGlobalIndex(SimInfo *info) {
772      SimInfo::MoleculeIterator mi;
773      Molecule::AtomIterator ai;
# Line 640 | Line 782 | namespace OpenMD {
782      int beginRigidBodyIndex;
783      int beginCutoffGroupIndex;
784      int nGlobalAtoms = info->getNGlobalAtoms();
643 –
644 –    /**@todo fixme */
645 – #ifndef IS_MPI
785      
786      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 –    
787      //rigidbody's index begins right after atom's
788 <    beginRigidBodyIndex += info->getNGlobalAtoms();
789 <    
790 <    for(mol = info->beginMolecule(mi); mol != NULL;
791 <        mol = info->nextMolecule(mi)) {
788 >    beginRigidBodyIndex = info->getNGlobalAtoms();
789 >    beginCutoffGroupIndex = 0;
790 >
791 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
792        
793 <      //local index(index in DataStorge) of atom is important
794 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
795 <        atom->setGlobalIndex(beginAtomIndex++);
793 > #ifdef IS_MPI      
794 >      if (info->getMolToProc(i) == worldRank) {
795 > #endif        
796 >        // stuff to do if I own this molecule
797 >        mol = info->getMoleculeByGlobalIndex(i);
798 >
799 >        //local index(index in DataStorge) of atom is important
800 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
801 >          atom->setGlobalIndex(beginAtomIndex++);
802 >        }
803 >        
804 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
805 >            rb = mol->nextRigidBody(ri)) {
806 >          rb->setGlobalIndex(beginRigidBodyIndex++);
807 >        }
808 >        
809 >        //local index of cutoff group is trivial, it only depends on
810 >        //the order of travesing
811 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
812 >            cg = mol->nextCutoffGroup(ci)) {
813 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
814 >        }        
815 >        
816 > #ifdef IS_MPI        
817 >      }  else {
818 >
819 >        // stuff to do if I don't own this molecule
820 >        
821 >        int stampId = info->getMoleculeStampId(i);
822 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
823 >
824 >        beginAtomIndex += stamp->getNAtoms();
825 >        beginRigidBodyIndex += stamp->getNRigidBodies();
826 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
827        }
828 <      
829 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
830 <          rb = mol->nextRigidBody(ri)) {
831 <        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 <    
828 > #endif          
829 >
830 >    } //end for(int i=0)  
831 >
832      //fill globalGroupMembership
833      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
834      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
# Line 721 | Line 840 | namespace OpenMD {
840          
841        }      
842      }
843 <    
843 >  
844   #ifdef IS_MPI    
845      // Since the globalGroupMembership has been zero filled and we've only
846      // poked values into the atoms we know, we can do an Allreduce
# Line 784 | Line 903 | namespace OpenMD {
903      for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
904        int myGlobalIndex = mol->getGlobalIndex();
905        int globalIO = startingIOIndexForMol[myGlobalIndex];
906 <      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
907 <           integrableObject = mol->nextIntegrableObject(ioi)) {
908 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
909 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
906 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
907 >           sd = mol->nextIntegrableObject(ioi)) {
908 >        sd->setGlobalIntegrableObjectIndex(globalIO);
909 >        IOIndexToIntegrableObject[globalIO] = sd;
910          globalIO++;
911        }
912      }
913 <    
913 >      
914      info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
915      
916    }
917    
918    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
919      Globals* simParams;
920 +
921      simParams = info->getSimParams();
922      
803 –    
923      DumpReader reader(info, mdFileName);
924      int nframes = reader.getNFrames();
925 <    
925 >
926      if (nframes > 0) {
927        reader.readFrame(nframes - 1);
928      } else {
# Line 814 | Line 933 | namespace OpenMD {
933        painCave.isFatal = 1;
934        simError();
935      }
817 –    
936      //copy the current snapshot to previous snapshot
937      info->getSnapshotManager()->advance();
938    }

Diff Legend

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