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 1810 by gezelter, Thu Nov 8 14:23:43 2012 UTC vs.
Revision 1976 by gezelter, Wed Mar 12 20:01:15 2014 UTC

# Line 1 | Line 1
1   /*
2 < * copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
5   * non-exclusive, royalty free, license to use, modify and
# 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 59 | Line 64
64   #include "brains/ForceField.hpp"
65   #include "utils/simError.h"
66   #include "utils/StringUtils.hpp"
67 + #include "utils/Revision.hpp"
68   #include "math/SeqRandNumGen.hpp"
69   #include "mdParser/MDLexer.hpp"
70   #include "mdParser/MDParser.hpp"
# Line 84 | Line 90
90   #include "types/FixedChargeAdapter.hpp"
91   #include "types/FluctuatingChargeAdapter.hpp"
92  
87 – #ifdef IS_MPI
88 – #include "mpi.h"
89 – #include "math/ParallelRandNumGen.hpp"
90 – #endif
93  
94   namespace OpenMD {
95    
# Line 102 | Line 104 | namespace OpenMD {
104        const int masterNode = 0;
105  
106        if (worldRank == masterNode) {
107 <        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
107 >        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
108 >        // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
109   #endif                
110          SimplePreprocessor preprocessor;
111 <        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock,
112 <                                ppStream);
111 >        preprocessor.preprocess(rawMetaDataStream, filename,
112 >                                startOfMetaDataBlock, ppStream);
113                  
114   #ifdef IS_MPI            
115 <        //brocasting the stream size
115 >        //broadcasting the stream size
116          streamSize = ppStream.str().size() +1;
117 <        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
118 <        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
119 <                              streamSize, MPI::CHAR, masterNode);
120 <                
117 >        MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
118 >        MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
119 >                  streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
120 >
121 >        // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
122 >        // MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
123 >        //                       streamSize, MPI::CHAR, masterNode);
124 >                          
125        } else {
126  
127 <        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
127 >        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
128 >        // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
129  
130          //get stream size
131 <        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
132 <
131 >        MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
132 >        // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
133          char* buf = new char[streamSize];
134          assert(buf);
135                  
136          //receive file content
137 <        MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode);
138 <                
137 >        MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
138 >        // MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode);
139 >
140          ppStream.str(buf);
141          delete [] buf;
133 –
142        }
143   #endif            
144        // Create a scanner that reads from the input stream
# Line 152 | Line 160 | namespace OpenMD {
160        parser.initializeASTFactory(factory);
161        parser.setASTFactory(&factory);
162        parser.mdfile();
155 –
163        // Create a tree parser that reads information into Globals
164        MDTreeParser treeParser;
165        treeParser.initializeASTFactory(factory);
# Line 256 | Line 263 | namespace OpenMD {
263      int metaDataBlockStart = -1;
264      int metaDataBlockEnd = -1;
265      int i, j;
266 <    streamoff mdOffset(0);
266 >    streamoff mdOffset;
267      int mdFileVersion;
268  
269      // Create a string for embedding the version information in the MetaData
# Line 266 | Line 273 | namespace OpenMD {
273      version.append(".");
274      version.append(OPENMD_VERSION_MINOR);
275  
276 <    std::string svnrev;
276 >    std::string svnrev(g_REVISION, strnlen(g_REVISION, 20));
277      //convert a macro from compiler to a string in c++
278 <    STR_DEFINE(svnrev, SVN_REV );
278 >    // STR_DEFINE(svnrev, SVN_REV );
279      version.append(" Revision: ");
280      // If there's no SVN revision, just call this the RELEASE revision.
281      if (!svnrev.empty()) {
# Line 515 | Line 522 | namespace OpenMD {
522    
523   #ifdef IS_MPI
524    void SimCreator::divideMolecules(SimInfo *info) {
518 –    RealType numerator;
519 –    RealType denominator;
520 –    RealType precast;
521 –    RealType x;
522 –    RealType y;
525      RealType a;
524 –    int old_atoms;
525 –    int add_atoms;
526 –    int new_atoms;
527 –    int nTarget;
528 –    int done;
529 –    int i;
530 –    int loops;
531 –    int which_proc;
526      int nProcessors;
527      std::vector<int> atomsPerProc;
528      int nGlobalMols = info->getNGlobalMolecules();
529 <    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
529 >    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an
530 >                                                    // error
531 >                                                    // condition:
532      
533 <    nProcessors = MPI::COMM_WORLD.Get_size();
533 >    MPI_Comm_size( MPI_COMM_WORLD, &nProcessors);    
534 >    //nProcessors = MPI::COMM_WORLD.Get_size();
535      
536      if (nProcessors > nGlobalMols) {
537        sprintf(painCave.errMsg,
# Line 543 | Line 540 | namespace OpenMD {
540                "\tthe number of molecules.  This will not result in a \n"
541                "\tusable division of atoms for force decomposition.\n"
542                "\tEither try a smaller number of processors, or run the\n"
543 <              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
543 >              "\tsingle-processor version of OpenMD.\n", nProcessors,
544 >              nGlobalMols);
545        
546        painCave.isFatal = 1;
547        simError();
548      }
549      
552 –    int seedValue;
550      Globals * simParams = info->getSimParams();
551 <    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
551 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel
552 >                             //random number generator
553      if (simParams->haveSeed()) {
554 <      seedValue = simParams->getSeed();
554 >      int seedValue = simParams->getSeed();
555        myRandom = new SeqRandNumGen(seedValue);
556      }else {
557        myRandom = new SeqRandNumGen();
# Line 566 | Line 564 | namespace OpenMD {
564      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
565      
566      if (worldRank == 0) {
567 <      numerator = info->getNGlobalAtoms();
568 <      denominator = nProcessors;
569 <      precast = numerator / denominator;
570 <      nTarget = (int)(precast + 0.5);
567 >      RealType numerator = info->getNGlobalAtoms();
568 >      RealType denominator = nProcessors;
569 >      RealType precast = numerator / denominator;
570 >      int nTarget = (int)(precast + 0.5);
571        
572 <      for(i = 0; i < nGlobalMols; i++) {
572 >      for(int i = 0; i < nGlobalMols; i++) {
573  
574 <        done = 0;
575 <        loops = 0;
574 >        int done = 0;
575 >        int loops = 0;
576          
577          while (!done) {
578            loops++;
579            
580            // Pick a processor at random
581            
582 <          which_proc = (int) (myRandom->rand() * nProcessors);
582 >          int which_proc = (int) (myRandom->rand() * nProcessors);
583            
584            //get the molecule stamp first
585            int stampId = info->getMoleculeStampId(i);
586            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
587            
588            // How many atoms does this processor have so far?
589 <          old_atoms = atomsPerProc[which_proc];
590 <          add_atoms = moleculeStamp->getNAtoms();
591 <          new_atoms = old_atoms + add_atoms;
589 >          int old_atoms = atomsPerProc[which_proc];
590 >          int add_atoms = moleculeStamp->getNAtoms();
591 >          int new_atoms = old_atoms + add_atoms;
592            
593            // If we've been through this loop too many times, we need
594            // to just give up and assign the molecule to this processor
# Line 634 | Line 632 | namespace OpenMD {
632            //           Pacc(x) = exp(- a * x)
633            // where a = penalty / (average atoms per molecule)
634            
635 <          x = (RealType)(new_atoms - nTarget);
636 <          y = myRandom->rand();
635 >          RealType x = (RealType)(new_atoms - nTarget);
636 >          RealType y = myRandom->rand();
637            
638            if (y < exp(- a * x)) {
639              molToProcMap[i] = which_proc;
# Line 652 | Line 650 | namespace OpenMD {
650        delete myRandom;
651  
652        // Spray out this nonsense to all other processors:
653 <      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
653 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
654 >      // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
655      } else {
656        
657        // Listen to your marching orders from processor 0:
658 <      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
658 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
659 >      // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
660  
661      }
662      
# Line 704 | Line 704 | namespace OpenMD {
704      set<AtomType*>::iterator i;
705      bool hasDirectionalAtoms = false;
706      bool hasFixedCharge = false;
707 <    bool hasMultipoles = false;    
707 >    bool hasDipoles = false;    
708 >    bool hasQuadrupoles = false;    
709      bool hasPolarizable = false;    
710      bool hasFluctuatingCharge = false;    
711      bool hasMetallic = false;
# Line 726 | Line 727 | namespace OpenMD {
727        if (da.isDirectional()){
728          hasDirectionalAtoms = true;
729        }
730 <      if (ma.isMultipole()){
731 <        hasMultipoles = true;
730 >      if (ma.isDipole()){
731 >        hasDipoles = true;
732        }
733 +      if (ma.isQuadrupole()){
734 +        hasQuadrupoles = true;
735 +      }
736        if (ea.isEAM() || sca.isSuttonChen()){
737          hasMetallic = true;
738        }
# Line 752 | Line 756 | namespace OpenMD {
756          storageLayout |= DataStorage::dslTorque;
757        }
758      }
759 <    if (hasMultipoles) {
760 <      storageLayout |= DataStorage::dslElectroFrame;
759 >    if (hasDipoles) {
760 >      storageLayout |= DataStorage::dslDipole;
761      }
762 +    if (hasQuadrupoles) {
763 +      storageLayout |= DataStorage::dslQuadrupole;
764 +    }
765      if (hasFixedCharge || hasFluctuatingCharge) {
766        storageLayout |= DataStorage::dslSkippedCharge;
767      }
# Line 789 | Line 796 | namespace OpenMD {
796        }
797      }
798  
799 <    if (simParams->getOutputElectricField()) {
799 >    if (simParams->getOutputElectricField() | simParams->haveElectricField()) {
800        storageLayout |= DataStorage::dslElectricField;
801      }
802 +
803      if (simParams->getOutputFluctuatingCharges()) {
804        storageLayout |= DataStorage::dslFlucQPosition;
805        storageLayout |= DataStorage::dslFlucQVelocity;
806        storageLayout |= DataStorage::dslFlucQForce;
807      }
808  
809 +    info->setStorageLayout(storageLayout);
810 +
811      return storageLayout;
812    }
813  
# Line 806 | Line 816 | namespace OpenMD {
816      Molecule::AtomIterator ai;
817      Molecule::RigidBodyIterator ri;
818      Molecule::CutoffGroupIterator ci;
819 +    Molecule::BondIterator boi;
820 +    Molecule::BendIterator bei;
821 +    Molecule::TorsionIterator ti;
822 +    Molecule::InversionIterator ii;
823      Molecule::IntegrableObjectIterator  ioi;
824 <    Molecule * mol;
825 <    Atom * atom;
826 <    RigidBody * rb;
827 <    CutoffGroup * cg;
824 >    Molecule* mol;
825 >    Atom* atom;
826 >    RigidBody* rb;
827 >    CutoffGroup* cg;
828 >    Bond* bond;
829 >    Bend* bend;
830 >    Torsion* torsion;
831 >    Inversion* inversion;
832      int beginAtomIndex;
833      int beginRigidBodyIndex;
834      int beginCutoffGroupIndex;
835 +    int beginBondIndex;
836 +    int beginBendIndex;
837 +    int beginTorsionIndex;
838 +    int beginInversionIndex;
839      int nGlobalAtoms = info->getNGlobalAtoms();
840      int nGlobalRigidBodies = info->getNGlobalRigidBodies();
841      
842      beginAtomIndex = 0;
843 <    //rigidbody's index begins right after atom's
843 >    // The rigid body indices begin immediately after the atom indices:
844      beginRigidBodyIndex = info->getNGlobalAtoms();
845      beginCutoffGroupIndex = 0;
846 <
846 >    beginBondIndex = 0;
847 >    beginBendIndex = 0;
848 >    beginTorsionIndex = 0;
849 >    beginInversionIndex = 0;
850 >  
851      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
852        
853   #ifdef IS_MPI      
# Line 830 | Line 856 | namespace OpenMD {
856          // stuff to do if I own this molecule
857          mol = info->getMoleculeByGlobalIndex(i);
858  
859 <        //local index(index in DataStorge) of atom is important
860 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
859 >        // The local index(index in DataStorge) of the atom is important:
860 >        for(atom = mol->beginAtom(ai); atom != NULL;
861 >            atom = mol->nextAtom(ai)) {
862            atom->setGlobalIndex(beginAtomIndex++);
863          }
864          
# Line 840 | Line 867 | namespace OpenMD {
867            rb->setGlobalIndex(beginRigidBodyIndex++);
868          }
869          
870 <        //local index of cutoff group is trivial, it only depends on
871 <        //the order of travesing
870 >        // The local index of other objects only depends on the order
871 >        // of traversal:
872          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
873              cg = mol->nextCutoffGroup(ci)) {
874            cg->setGlobalIndex(beginCutoffGroupIndex++);
875          }        
876 +        for(bond = mol->beginBond(boi); bond != NULL;
877 +            bond = mol->nextBond(boi)) {
878 +          bond->setGlobalIndex(beginBondIndex++);
879 +        }        
880 +        for(bend = mol->beginBend(bei); bend != NULL;
881 +            bend = mol->nextBend(bei)) {
882 +          bend->setGlobalIndex(beginBendIndex++);
883 +        }        
884 +        for(torsion = mol->beginTorsion(ti); torsion != NULL;
885 +            torsion = mol->nextTorsion(ti)) {
886 +          torsion->setGlobalIndex(beginTorsionIndex++);
887 +        }        
888 +        for(inversion = mol->beginInversion(ii); inversion != NULL;
889 +            inversion = mol->nextInversion(ii)) {
890 +          inversion->setGlobalIndex(beginInversionIndex++);
891 +        }        
892          
893   #ifdef IS_MPI        
894        }  else {
# Line 858 | Line 901 | namespace OpenMD {
901          beginAtomIndex += stamp->getNAtoms();
902          beginRigidBodyIndex += stamp->getNRigidBodies();
903          beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
904 +        beginBondIndex += stamp->getNBonds();
905 +        beginBendIndex += stamp->getNBends();
906 +        beginTorsionIndex += stamp->getNTorsions();
907 +        beginInversionIndex += stamp->getNInversions();
908        }
909   #endif          
910  
# Line 865 | Line 912 | namespace OpenMD {
912  
913      //fill globalGroupMembership
914      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
915 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
916 <      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
917 <        
915 >    for(mol = info->beginMolecule(mi); mol != NULL;
916 >        mol = info->nextMolecule(mi)) {        
917 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
918 >           cg = mol->nextCutoffGroup(ci)) {        
919          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
920            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
921          }
# Line 882 | Line 930 | namespace OpenMD {
930      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
931      // docs said we could.
932      std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
933 <    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
934 <                              &tmpGroupMembership[0], nGlobalAtoms,
935 <                              MPI::INT, MPI::SUM);
933 >    MPI_Allreduce(&globalGroupMembership[0],
934 >                  &tmpGroupMembership[0], nGlobalAtoms,
935 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
936 >    // MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
937 >    //                           &tmpGroupMembership[0], nGlobalAtoms,
938 >    //                           MPI::INT, MPI::SUM);
939      info->setGlobalGroupMembership(tmpGroupMembership);
940   #else
941      info->setGlobalGroupMembership(globalGroupMembership);
# Line 908 | Line 959 | namespace OpenMD {
959   #ifdef IS_MPI
960      std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
961                                        info->getNGlobalRigidBodies(), 0);
962 <    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
963 <                              nGlobalAtoms + nGlobalRigidBodies,
964 <                              MPI::INT, MPI::SUM);
962 >    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
963 >                  nGlobalAtoms + nGlobalRigidBodies,
964 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
965 >    // MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
966 >    //                           nGlobalAtoms + nGlobalRigidBodies,
967 >    //                           MPI::INT, MPI::SUM);
968      
969      info->setGlobalMolMembership(tmpMolMembership);
970   #else
# Line 928 | Line 982 | namespace OpenMD {
982      
983   #ifdef IS_MPI
984      std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
985 <    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
986 <                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
985 >    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
986 >      info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
987 >    // MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
988 >    //                           info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
989   #else
990      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
991   #endif    
# Line 960 | Line 1016 | namespace OpenMD {
1016    }
1017    
1018    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
1019 <
1019 >    
1020      DumpReader reader(info, mdFileName);
1021      int nframes = reader.getNFrames();
1022 <
1022 >    
1023      if (nframes > 0) {
1024        reader.readFrame(nframes - 1);
1025      } else {

Diff Legend

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