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 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
Revision 1827 by gezelter, Wed Jan 9 19:49:07 2013 UTC

# Line 44 | Line 44
44   * @file SimCreator.cpp
45   * @author tlin
46   * @date 11/03/2004
47 – * @time 13:51am
47   * @version 1.0
48   */
49   #include <exception>
# Line 56 | 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 76 | 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"
# Line 92 | Line 99 | namespace OpenMD {
99   #ifdef IS_MPI            
100        int streamSize;
101        const int masterNode = 0;
102 <      int commStatus;
102 >
103        if (worldRank == masterNode) {
104 <        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
104 >        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
105   #endif                
106          SimplePreprocessor preprocessor;
107          preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
# Line 102 | Line 109 | namespace OpenMD {
109   #ifdef IS_MPI            
110          //brocasting the stream size
111          streamSize = ppStream.str().size() +1;
112 <        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
113 <
114 <        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
108 <            
109 <                
112 >        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
113 >        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI::CHAR, masterNode);
114 >                          
115        } else {
116 +        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
117  
112 –        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
113 –
118          //get stream size
119 <        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
119 >        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
120  
121          char* buf = new char[streamSize];
122          assert(buf);
123                  
124          //receive file content
125 <        commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
125 >        MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode);
126                  
127          ppStream.str(buf);
128          delete [] buf;
125 –
129        }
130   #endif            
131        // Create a scanner that reads from the input stream
# Line 247 | Line 250 | namespace OpenMD {
250      std::string mdRawData;
251      int metaDataBlockStart = -1;
252      int metaDataBlockEnd = -1;
253 <    int i;
254 <    int mdOffset;
253 >    int i, j;
254 >    streamoff mdOffset;
255      int mdFileVersion;
256  
257 +    // Create a string for embedding the version information in the MetaData
258 +    std::string version;
259 +    version.assign("## Last run using OpenMD Version: ");
260 +    version.append(OPENMD_VERSION_MAJOR);
261 +    version.append(".");
262 +    version.append(OPENMD_VERSION_MINOR);
263 +
264 +    std::string svnrev;
265 +    //convert a macro from compiler to a string in c++
266 +    STR_DEFINE(svnrev, SVN_REV );
267 +    version.append(" Revision: ");
268 +    // If there's no SVN revision, just call this the RELEASE revision.
269 +    if (!svnrev.empty()) {
270 +      version.append(svnrev);
271 +    } else {
272 +      version.append("RELEASE");
273 +    }
274 +  
275   #ifdef IS_MPI            
276      const int masterNode = 0;
277      if (worldRank == masterNode) {
278   #endif
279  
280 <      std::ifstream mdFile_(mdFileName.c_str());
280 >      std::ifstream mdFile_;
281 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
282        
283        if (mdFile_.fail()) {
284          sprintf(painCave.errMsg,
# Line 344 | Line 366 | namespace OpenMD {
366  
367        mdRawData.clear();
368  
369 +      bool foundVersion = false;
370 +
371        for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
372          mdFile_.getline(buffer, bufferSize);
373 <        mdRawData += buffer;
373 >        std::string line = trimLeftCopy(buffer);
374 >        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version");
375 >        if (static_cast<size_t>(j) != string::npos) {
376 >          foundVersion = true;
377 >          mdRawData += version;
378 >        } else {
379 >          mdRawData += buffer;
380 >        }
381          mdRawData += "\n";
382        }
383 <
383 >      
384 >      if (!foundVersion) mdRawData += version + "\n";
385 >      
386        mdFile_.close();
387  
388   #ifdef IS_MPI
# Line 363 | Line 396 | namespace OpenMD {
396                                     metaDataBlockStart + 1);
397      
398      //create the force field
399 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
399 >    ForceField * ff = new ForceField(simParams->getForceField());
400  
401      if (ff == NULL) {
402        sprintf(painCave.errMsg,
# Line 414 | Line 447 | namespace OpenMD {
447      //create the molecules
448      createMolecules(info);
449      
450 +    //find the storage layout
451 +
452 +    int storageLayout = computeStorageLayout(info);
453 +
454      //allocate memory for DataStorage(circular reference, need to
455      //break it)
456 <    info->setSnapshotManager(new SimSnapshotManager(info));
456 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
457      
458      //set the global index of atoms, rigidbodies and cutoffgroups
459      //(only need to be set once, the global index will never change
# Line 485 | Line 522 | namespace OpenMD {
522      int nTarget;
523      int done;
524      int i;
488 –    int j;
525      int loops;
526      int which_proc;
527      int nProcessors;
# Line 493 | Line 529 | namespace OpenMD {
529      int nGlobalMols = info->getNGlobalMolecules();
530      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
531      
532 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
532 >    nProcessors = MPI::COMM_WORLD.Get_size();
533      
534      if (nProcessors > nGlobalMols) {
535        sprintf(painCave.errMsg,
# Line 531 | Line 567 | namespace OpenMD {
567        nTarget = (int)(precast + 0.5);
568        
569        for(i = 0; i < nGlobalMols; i++) {
570 +
571          done = 0;
572          loops = 0;
573          
# Line 555 | Line 592 | namespace OpenMD {
592            // and be done with it.
593            
594            if (loops > 100) {
595 +
596              sprintf(painCave.errMsg,
597 <                    "I've tried 100 times to assign molecule %d to a "
598 <                    " processor, but can't find a good spot.\n"
599 <                    "I'm assigning it at random to processor %d.\n",
600 <                    i, which_proc);
601 <            
597 >                    "There have been 100 attempts to assign molecule %d to an\n"
598 >                    "\tunderworked processor, but there's no good place to\n"
599 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
600 >                    i, which_proc);
601 >          
602              painCave.isFatal = 0;
603 +            painCave.severity = OPENMD_INFO;
604              simError();
605              
606              molToProcMap[i] = which_proc;
# Line 606 | Line 645 | namespace OpenMD {
645        }
646        
647        delete myRandom;
648 <      
648 >
649        // Spray out this nonsense to all other processors:
650 <      
612 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
650 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
651      } else {
652        
653        // Listen to your marching orders from processor 0:
654 <      
655 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
654 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
655 >
656      }
657      
658      info->setMolToProcMap(molToProcMap);
# Line 652 | Line 690 | namespace OpenMD {
690        
691      } //end for(int i=0)  
692    }
693 +    
694 +  int SimCreator::computeStorageLayout(SimInfo* info) {
695 +
696 +    Globals* simParams = info->getSimParams();
697 +    int nRigidBodies = info->getNGlobalRigidBodies();
698 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
699 +    set<AtomType*>::iterator i;
700 +    bool hasDirectionalAtoms = false;
701 +    bool hasFixedCharge = false;
702 +    bool hasDipoles = false;    
703 +    bool hasQuadrupoles = false;    
704 +    bool hasPolarizable = false;    
705 +    bool hasFluctuatingCharge = false;    
706 +    bool hasMetallic = false;
707 +    int storageLayout = 0;
708 +    storageLayout |= DataStorage::dslPosition;
709 +    storageLayout |= DataStorage::dslVelocity;
710 +    storageLayout |= DataStorage::dslForce;
711 +
712 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
713 +
714 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
715 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
716 +      EAMAdapter ea = EAMAdapter( (*i) );
717 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
718 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
719 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
720 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
721 +
722 +      if (da.isDirectional()){
723 +        hasDirectionalAtoms = true;
724 +      }
725 +      if (ma.isDipole()){
726 +        hasDipoles = true;
727 +      }
728 +      if (ma.isQuadrupole()){
729 +        hasQuadrupoles = true;
730 +      }
731 +      if (ea.isEAM() || sca.isSuttonChen()){
732 +        hasMetallic = true;
733 +      }
734 +      if ( fca.isFixedCharge() ){
735 +        hasFixedCharge = true;
736 +      }
737 +      if ( fqa.isFluctuatingCharge() ){
738 +        hasFluctuatingCharge = true;
739 +      }
740 +      if ( pa.isPolarizable() ){
741 +        hasPolarizable = true;
742 +      }
743 +    }
744      
745 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
746 +      storageLayout |= DataStorage::dslAmat;
747 +      if(storageLayout & DataStorage::dslVelocity) {
748 +        storageLayout |= DataStorage::dslAngularMomentum;
749 +      }
750 +      if (storageLayout & DataStorage::dslForce) {
751 +        storageLayout |= DataStorage::dslTorque;
752 +      }
753 +    }
754 +    if (hasDipoles) {
755 +      storageLayout |= DataStorage::dslDipole;
756 +    }
757 +    if (hasQuadrupoles) {
758 +      storageLayout |= DataStorage::dslQuadrupole;
759 +    }
760 +    if (hasFixedCharge || hasFluctuatingCharge) {
761 +      storageLayout |= DataStorage::dslSkippedCharge;
762 +    }
763 +    if (hasMetallic) {
764 +      storageLayout |= DataStorage::dslDensity;
765 +      storageLayout |= DataStorage::dslFunctional;
766 +      storageLayout |= DataStorage::dslFunctionalDerivative;
767 +    }
768 +    if (hasPolarizable) {
769 +      storageLayout |= DataStorage::dslElectricField;
770 +    }
771 +    if (hasFluctuatingCharge){
772 +      storageLayout |= DataStorage::dslFlucQPosition;
773 +      if(storageLayout & DataStorage::dslVelocity) {
774 +        storageLayout |= DataStorage::dslFlucQVelocity;
775 +      }
776 +      if (storageLayout & DataStorage::dslForce) {
777 +        storageLayout |= DataStorage::dslFlucQForce;
778 +      }
779 +    }
780 +    
781 +    // if the user has asked for them, make sure we've got the memory for the
782 +    // objects defined.
783 +
784 +    if (simParams->getOutputParticlePotential()) {
785 +      storageLayout |= DataStorage::dslParticlePot;
786 +    }
787 +
788 +    if (simParams->havePrintHeatFlux()) {
789 +      if (simParams->getPrintHeatFlux()) {
790 +        storageLayout |= DataStorage::dslParticlePot;
791 +      }
792 +    }
793 +
794 +    if (simParams->getOutputElectricField()) {
795 +      storageLayout |= DataStorage::dslElectricField;
796 +    }
797 +
798 +    if (simParams->getOutputFluctuatingCharges()) {
799 +      storageLayout |= DataStorage::dslFlucQPosition;
800 +      storageLayout |= DataStorage::dslFlucQVelocity;
801 +      storageLayout |= DataStorage::dslFlucQForce;
802 +    }
803 +
804 +    return storageLayout;
805 +  }
806 +
807    void SimCreator::setGlobalIndex(SimInfo *info) {
808      SimInfo::MoleculeIterator mi;
809      Molecule::AtomIterator ai;
# Line 667 | Line 818 | namespace OpenMD {
818      int beginRigidBodyIndex;
819      int beginCutoffGroupIndex;
820      int nGlobalAtoms = info->getNGlobalAtoms();
821 +    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
822      
823      beginAtomIndex = 0;
824 <    beginRigidBodyIndex = 0;
824 >    //rigidbody's index begins right after atom's
825 >    beginRigidBodyIndex = info->getNGlobalAtoms();
826      beginCutoffGroupIndex = 0;
827  
828      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
# Line 732 | Line 885 | namespace OpenMD {
885      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
886      // docs said we could.
887      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
888 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
889 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
888 >    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
889 >                              &tmpGroupMembership[0], nGlobalAtoms,
890 >                              MPI::INT, MPI::SUM);
891      info->setGlobalGroupMembership(tmpGroupMembership);
892   #else
893      info->setGlobalGroupMembership(globalGroupMembership);
894   #endif
895      
896      //fill molMembership
897 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
897 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
898 >                                         info->getNGlobalRigidBodies(), 0);
899      
900 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
900 >    for(mol = info->beginMolecule(mi); mol != NULL;
901 >        mol = info->nextMolecule(mi)) {
902        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
903          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
904        }
905 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
906 +           rb = mol->nextRigidBody(ri)) {
907 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
908 +      }
909      }
910      
911   #ifdef IS_MPI
912 <    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
912 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
913 >                                      info->getNGlobalRigidBodies(), 0);
914 >    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
915 >                              nGlobalAtoms + nGlobalRigidBodies,
916 >                              MPI::INT, MPI::SUM);
917      
754 –    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
755 –                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
756 –    
918      info->setGlobalMolMembership(tmpMolMembership);
919   #else
920      info->setGlobalMolMembership(globalMolMembership);
# Line 763 | Line 924 | namespace OpenMD {
924      // here the molecules are listed by their global indices.
925  
926      std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
927 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
927 >    for (mol = info->beginMolecule(mi); mol != NULL;
928 >         mol = info->nextMolecule(mi)) {
929        nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
930      }
931      
932   #ifdef IS_MPI
933      std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
934 <    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
935 <                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
934 >    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
935 >                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
936   #else
937      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
938   #endif    
# Line 784 | Line 946 | namespace OpenMD {
946      }
947      
948      std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
949 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
949 >    for (mol = info->beginMolecule(mi); mol != NULL;
950 >         mol = info->nextMolecule(mi)) {
951        int myGlobalIndex = mol->getGlobalIndex();
952        int globalIO = startingIOIndexForMol[myGlobalIndex];
953 <      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
954 <           integrableObject = mol->nextIntegrableObject(ioi)) {
955 <        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
956 <        IOIndexToIntegrableObject[globalIO] = integrableObject;
953 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
954 >           sd = mol->nextIntegrableObject(ioi)) {
955 >        sd->setGlobalIntegrableObjectIndex(globalIO);
956 >        IOIndexToIntegrableObject[globalIO] = sd;
957          globalIO++;
958        }
959      }
# Line 800 | Line 963 | namespace OpenMD {
963    }
964    
965    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
803 –    Globals* simParams;
804 –
805 –    simParams = info->getSimParams();
966      
967      DumpReader reader(info, mdFileName);
968      int nframes = reader.getNFrames();
969 <
969 >    
970      if (nframes > 0) {
971        reader.readFrame(nframes - 1);
972      } else {

Diff Legend

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