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 1793 by gezelter, Fri Aug 31 21:16:10 2012 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

# Line 35 | Line 35
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# 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 112 | Line 111 | namespace OpenMD {
111          //brocasting the stream size
112          streamSize = ppStream.str().size() +1;
113          MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
114 <        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
115 <                              streamSize, MPI::CHAR, masterNode);
117 <                
114 >        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI::CHAR, masterNode);
115 >                          
116        } else {
119 –
117          MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
118  
119          //get stream size
# Line 130 | Line 127 | namespace OpenMD {
127                  
128          ppStream.str(buf);
129          delete [] buf;
133 –
130        }
131   #endif            
132        // Create a scanner that reads from the input stream
# Line 255 | Line 251 | namespace OpenMD {
251      std::string mdRawData;
252      int metaDataBlockStart = -1;
253      int metaDataBlockEnd = -1;
254 <    int i;
255 <    streamoff mdOffset(0);
254 >    int i, j;
255 >    streamoff mdOffset;
256      int mdFileVersion;
257  
258 +    // Create a string for embedding the version information in the MetaData
259 +    std::string version;
260 +    version.assign("## Last run using OpenMD Version: ");
261 +    version.append(OPENMD_VERSION_MAJOR);
262 +    version.append(".");
263 +    version.append(OPENMD_VERSION_MINOR);
264  
265 +    std::string svnrev;
266 +    //convert a macro from compiler to a string in c++
267 +    STR_DEFINE(svnrev, SVN_REV );
268 +    version.append(" Revision: ");
269 +    // If there's no SVN revision, just call this the RELEASE revision.
270 +    if (!svnrev.empty()) {
271 +      version.append(svnrev);
272 +    } else {
273 +      version.append("RELEASE");
274 +    }
275 +  
276   #ifdef IS_MPI            
277      const int masterNode = 0;
278      if (worldRank == masterNode) {
# Line 354 | Line 367 | namespace OpenMD {
367  
368        mdRawData.clear();
369  
370 +      bool foundVersion = false;
371 +
372        for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
373          mdFile_.getline(buffer, bufferSize);
374 <        mdRawData += buffer;
374 >        std::string line = trimLeftCopy(buffer);
375 >        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version");
376 >        if (static_cast<size_t>(j) != string::npos) {
377 >          foundVersion = true;
378 >          mdRawData += version;
379 >        } else {
380 >          mdRawData += buffer;
381 >        }
382          mdRawData += "\n";
383        }
384 <
384 >      
385 >      if (!foundVersion) mdRawData += version + "\n";
386 >      
387        mdFile_.close();
388  
389   #ifdef IS_MPI
# Line 487 | Line 511 | namespace OpenMD {
511    
512   #ifdef IS_MPI
513    void SimCreator::divideMolecules(SimInfo *info) {
490 –    RealType numerator;
491 –    RealType denominator;
492 –    RealType precast;
493 –    RealType x;
494 –    RealType y;
514      RealType a;
496 –    int old_atoms;
497 –    int add_atoms;
498 –    int new_atoms;
499 –    int nTarget;
500 –    int done;
501 –    int i;
502 –    int loops;
503 –    int which_proc;
515      int nProcessors;
516      std::vector<int> atomsPerProc;
517      int nGlobalMols = info->getNGlobalMolecules();
518 <    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
518 >    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an
519 >                                                    // error
520 >                                                    // condition:
521      
522 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
522 >    nProcessors = MPI::COMM_WORLD.Get_size();
523      
524      if (nProcessors > nGlobalMols) {
525        sprintf(painCave.errMsg,
# Line 515 | Line 528 | namespace OpenMD {
528                "\tthe number of molecules.  This will not result in a \n"
529                "\tusable division of atoms for force decomposition.\n"
530                "\tEither try a smaller number of processors, or run the\n"
531 <              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
531 >              "\tsingle-processor version of OpenMD.\n", nProcessors,
532 >              nGlobalMols);
533        
534        painCave.isFatal = 1;
535        simError();
536      }
537      
524 –    int seedValue;
538      Globals * simParams = info->getSimParams();
539 <    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
539 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel
540 >                             //random number generator
541      if (simParams->haveSeed()) {
542 <      seedValue = simParams->getSeed();
542 >      int seedValue = simParams->getSeed();
543        myRandom = new SeqRandNumGen(seedValue);
544      }else {
545        myRandom = new SeqRandNumGen();
# Line 538 | Line 552 | namespace OpenMD {
552      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
553      
554      if (worldRank == 0) {
555 <      numerator = info->getNGlobalAtoms();
556 <      denominator = nProcessors;
557 <      precast = numerator / denominator;
558 <      nTarget = (int)(precast + 0.5);
555 >      RealType numerator = info->getNGlobalAtoms();
556 >      RealType denominator = nProcessors;
557 >      RealType precast = numerator / denominator;
558 >      int nTarget = (int)(precast + 0.5);
559        
560 <      for(i = 0; i < nGlobalMols; i++) {
561 <        done = 0;
562 <        loops = 0;
560 >      for(int i = 0; i < nGlobalMols; i++) {
561 >
562 >        int done = 0;
563 >        int loops = 0;
564          
565          while (!done) {
566            loops++;
567            
568            // Pick a processor at random
569            
570 <          which_proc = (int) (myRandom->rand() * nProcessors);
570 >          int which_proc = (int) (myRandom->rand() * nProcessors);
571            
572            //get the molecule stamp first
573            int stampId = info->getMoleculeStampId(i);
574            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
575            
576            // How many atoms does this processor have so far?
577 <          old_atoms = atomsPerProc[which_proc];
578 <          add_atoms = moleculeStamp->getNAtoms();
579 <          new_atoms = old_atoms + add_atoms;
577 >          int old_atoms = atomsPerProc[which_proc];
578 >          int add_atoms = moleculeStamp->getNAtoms();
579 >          int new_atoms = old_atoms + add_atoms;
580            
581            // If we've been through this loop too many times, we need
582            // to just give up and assign the molecule to this processor
583            // and be done with it.
584            
585            if (loops > 100) {
586 +
587              sprintf(painCave.errMsg,
588 <                    "I've tried 100 times to assign molecule %d to a "
589 <                    " processor, but can't find a good spot.\n"
590 <                    "I'm assigning it at random to processor %d.\n",
588 >                    "There have been 100 attempts to assign molecule %d to an\n"
589 >                    "\tunderworked processor, but there's no good place to\n"
590 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
591                      i, which_proc);
592 <            
592 >          
593              painCave.isFatal = 0;
594 +            painCave.severity = OPENMD_INFO;
595              simError();
596              
597              molToProcMap[i] = which_proc;
# Line 603 | Line 620 | namespace OpenMD {
620            //           Pacc(x) = exp(- a * x)
621            // where a = penalty / (average atoms per molecule)
622            
623 <          x = (RealType)(new_atoms - nTarget);
624 <          y = myRandom->rand();
623 >          RealType x = (RealType)(new_atoms - nTarget);
624 >          RealType y = myRandom->rand();
625            
626            if (y < exp(- a * x)) {
627              molToProcMap[i] = which_proc;
# Line 619 | Line 636 | namespace OpenMD {
636        }
637        
638        delete myRandom;
639 <      
639 >
640        // Spray out this nonsense to all other processors:
641 <      
625 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
641 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
642      } else {
643        
644        // Listen to your marching orders from processor 0:
645 <      
646 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
645 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
646 >
647      }
648      
649      info->setMolToProcMap(molToProcMap);
# Line 674 | Line 690 | namespace OpenMD {
690      set<AtomType*>::iterator i;
691      bool hasDirectionalAtoms = false;
692      bool hasFixedCharge = false;
693 <    bool hasMultipoles = false;    
693 >    bool hasDipoles = false;    
694 >    bool hasQuadrupoles = false;    
695      bool hasPolarizable = false;    
696      bool hasFluctuatingCharge = false;    
697      bool hasMetallic = false;
# Line 696 | Line 713 | namespace OpenMD {
713        if (da.isDirectional()){
714          hasDirectionalAtoms = true;
715        }
716 <      if (ma.isMultipole()){
717 <        hasMultipoles = true;
716 >      if (ma.isDipole()){
717 >        hasDipoles = true;
718        }
719 +      if (ma.isQuadrupole()){
720 +        hasQuadrupoles = true;
721 +      }
722        if (ea.isEAM() || sca.isSuttonChen()){
723          hasMetallic = true;
724        }
# Line 722 | Line 742 | namespace OpenMD {
742          storageLayout |= DataStorage::dslTorque;
743        }
744      }
745 <    if (hasMultipoles) {
746 <      storageLayout |= DataStorage::dslElectroFrame;
745 >    if (hasDipoles) {
746 >      storageLayout |= DataStorage::dslDipole;
747      }
748 +    if (hasQuadrupoles) {
749 +      storageLayout |= DataStorage::dslQuadrupole;
750 +    }
751      if (hasFixedCharge || hasFluctuatingCharge) {
752        storageLayout |= DataStorage::dslSkippedCharge;
753      }
# Line 759 | Line 782 | namespace OpenMD {
782        }
783      }
784  
785 <    if (simParams->getOutputElectricField()) {
785 >    if (simParams->getOutputElectricField() | simParams->haveElectricField()) {
786        storageLayout |= DataStorage::dslElectricField;
787      }
788 +
789      if (simParams->getOutputFluctuatingCharges()) {
790        storageLayout |= DataStorage::dslFlucQPosition;
791        storageLayout |= DataStorage::dslFlucQVelocity;
792        storageLayout |= DataStorage::dslFlucQForce;
793      }
794  
795 +    info->setStorageLayout(storageLayout);
796 +
797      return storageLayout;
798    }
799  
# Line 785 | Line 811 | namespace OpenMD {
811      int beginRigidBodyIndex;
812      int beginCutoffGroupIndex;
813      int nGlobalAtoms = info->getNGlobalAtoms();
814 +    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
815      
816      beginAtomIndex = 0;
817      //rigidbody's index begins right after atom's
# Line 851 | Line 878 | namespace OpenMD {
878      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
879      // docs said we could.
880      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
881 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
882 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
881 >    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
882 >                              &tmpGroupMembership[0], nGlobalAtoms,
883 >                              MPI::INT, MPI::SUM);
884      info->setGlobalGroupMembership(tmpGroupMembership);
885   #else
886      info->setGlobalGroupMembership(globalGroupMembership);
887   #endif
888      
889      //fill molMembership
890 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
890 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
891 >                                         info->getNGlobalRigidBodies(), 0);
892      
893 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
893 >    for(mol = info->beginMolecule(mi); mol != NULL;
894 >        mol = info->nextMolecule(mi)) {
895        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
896          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
897        }
898 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
899 +           rb = mol->nextRigidBody(ri)) {
900 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
901 +      }
902      }
903      
904   #ifdef IS_MPI
905 <    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
905 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
906 >                                      info->getNGlobalRigidBodies(), 0);
907 >    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
908 >                              nGlobalAtoms + nGlobalRigidBodies,
909 >                              MPI::INT, MPI::SUM);
910      
873 –    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
874 –                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
875 –    
911      info->setGlobalMolMembership(tmpMolMembership);
912   #else
913      info->setGlobalMolMembership(globalMolMembership);
# Line 882 | Line 917 | namespace OpenMD {
917      // here the molecules are listed by their global indices.
918  
919      std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
920 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
920 >    for (mol = info->beginMolecule(mi); mol != NULL;
921 >         mol = info->nextMolecule(mi)) {
922        nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
923      }
924      
925   #ifdef IS_MPI
926      std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
927 <    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
928 <                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
927 >    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
928 >                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
929   #else
930      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
931   #endif    
# Line 903 | Line 939 | namespace OpenMD {
939      }
940      
941      std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
942 <    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
942 >    for (mol = info->beginMolecule(mi); mol != NULL;
943 >         mol = info->nextMolecule(mi)) {
944        int myGlobalIndex = mol->getGlobalIndex();
945        int globalIO = startingIOIndexForMol[myGlobalIndex];
946        for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
# Line 919 | Line 956 | namespace OpenMD {
956    }
957    
958    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
959 <
959 >    
960      DumpReader reader(info, mdFileName);
961      int nframes = reader.getNFrames();
962 <
962 >    
963      if (nframes > 0) {
964        reader.readFrame(nframes - 1);
965      } else {

Diff Legend

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