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

Diff Legend

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