ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/SimCreator.cpp
(Generate patch)

Comparing:
trunk/src/brains/SimCreator.cpp (file contents), Revision 1024 by tim, Wed Aug 30 18:42:29 2006 UTC vs.
branches/development/src/brains/SimCreator.cpp (file contents), Revision 1600 by gezelter, Wed Aug 3 20:20:37 2011 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   /**
# Line 79 | Line 79
79   #include "math/ParallelRandNumGen.hpp"
80   #endif
81  
82 < namespace oopse {
82 > namespace OpenMD {
83    
84    Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
85      Globals* simParams = NULL;
# Line 108 | Line 108 | namespace oopse {
108        } else {
109          //get stream size
110          commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
111 <                
111 >
112          char* buf = new char[streamSize];
113          assert(buf);
114                  
# Line 116 | Line 116 | namespace oopse {
116          commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
117                  
118          ppStream.str(buf);
119 <        delete buf;
119 >        delete [] buf;
120  
121        }
122   #endif            
# Line 145 | Line 145 | namespace oopse {
145        treeParser.initializeASTFactory(factory);
146        treeParser.setASTFactory(&factory);
147        simParams = treeParser.walkTree(parser.getAST());
148
148      }
149  
150        
# Line 215 | Line 214 | namespace oopse {
214        painCave.isFatal = 1;
215        simError();        
216      }
217 <    catch (OOPSEException& e) {
217 >    catch (OpenMDException& e) {
218        sprintf(painCave.errMsg,
219                "%s\n",
220                e.getMessage().c_str());
# Line 235 | Line 234 | namespace oopse {
234    
235    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
236                                    bool loadInitCoords) {
237 <
237 >    
238      const int bufferSize = 65535;
239      char buffer[bufferSize];
240      int lineNo = 0;
# Line 263 | Line 262 | namespace oopse {
262        mdFile_.getline(buffer, bufferSize);
263        ++lineNo;
264        std::string line = trimLeftCopy(buffer);
265 <      i = CaseInsensitiveFind(line, "<OOPSE");
266 <      if (i == string::npos) {
265 >      i = CaseInsensitiveFind(line, "<OpenMD");
266 >      if (static_cast<size_t>(i) == string::npos) {
267 >        // try the older file strings to see if that works:
268 >        i = CaseInsensitiveFind(line, "<OOPSE");
269 >      }
270 >      
271 >      if (static_cast<size_t>(i) == string::npos) {
272 >        // still no luck!
273          sprintf(painCave.errMsg,
274 <                "SimCreator: File: %s is not an OOPSE file!\n",
274 >                "SimCreator: File: %s is not a valid OpenMD file!\n",
275                  mdFileName.c_str());
276          painCave.isFatal = 1;
277          simError();
# Line 330 | Line 335 | namespace oopse {
335      Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
336      
337      //create the force field
338 <    ForceField * ff = ForceFieldFactory::getInstance()
339 <      ->createForceField(simParams->getForceField());
335 <    
338 >    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
339 >
340      if (ff == NULL) {
341        sprintf(painCave.errMsg,
342                "ForceField Factory can not create %s force field\n",
# Line 365 | Line 369 | namespace oopse {
369      }
370      
371      ff->parse(forcefieldFileName);
368    ff->setFortranForceOptions();
372      //create SimInfo
373      SimInfo * info = new SimInfo(ff, simParams);
374  
# Line 383 | Line 386 | namespace oopse {
386      //create the molecules
387      createMolecules(info);
388      
386    
389      //allocate memory for DataStorage(circular reference, need to
390      //break it)
391      info->setSnapshotManager(new SimSnapshotManager(info));
# Line 395 | Line 397 | namespace oopse {
397      //responsibility to LocalIndexManager.
398      setGlobalIndex(info);
399      
400 <    //Although addExcludePairs is called inside SimInfo's addMolecule
400 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
401      //method, at that point atoms don't have the global index yet
402      //(their global index are all initialized to -1).  Therefore we
403 <    //have to call addExcludePairs explicitly here. A way to work
403 >    //have to call addInteractionPairs explicitly here. A way to work
404      //around is that we can determine the beginning global indices of
405      //atoms before they get created.
406      SimInfo::MoleculeIterator mi;
407      Molecule* mol;
408      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
409 <      info->addExcludePairs(mol);
409 >      info->addInteractionPairs(mol);
410      }
411      
412      if (loadInitCoords)
413        loadCoordinates(info, mdFileName);    
412    
414      return info;
415    }
416    
# Line 473 | Line 474 | namespace oopse {
474                "\tthe number of molecules.  This will not result in a \n"
475                "\tusable division of atoms for force decomposition.\n"
476                "\tEither try a smaller number of processors, or run the\n"
477 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
477 >              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
478        
479        painCave.isFatal = 1;
480        simError();
# Line 588 | Line 589 | namespace oopse {
589        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
590      }
591      
592 +    cerr << "molToProcMap:\n";
593 +    for (int i = 0; i < molToProcMap.size(); i++) {
594 +      cerr << "m = " << i << " mtpr[m] = " << molToProcMap[i] <<"\n";
595 +    }
596 +
597      info->setMolToProcMap(molToProcMap);
598      sprintf(checkPointMsg,
599              "Successfully divided the molecules among the processors.\n");
600 <    MPIcheckPoint();
600 >    errorCheckPoint();
601    }
602    
603   #endif
# Line 608 | Line 614 | namespace oopse {
614   #endif
615          
616          stampId = info->getMoleculeStampId(i);
617 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
618 <                                                   stampId, i, info->getLocalIndexManager());
617 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
618 >                                                   info->getMoleculeStamp(stampId),
619 >                                                   stampId, i,
620 >                                                   info->getLocalIndexManager());
621          
622          info->addMolecule(mol);
623          
# Line 636 | Line 644 | namespace oopse {
644      int beginRigidBodyIndex;
645      int beginCutoffGroupIndex;
646      int nGlobalAtoms = info->getNGlobalAtoms();
639
640    /**@todo fixme */
641 #ifndef IS_MPI
647      
648      beginAtomIndex = 0;
649      beginRigidBodyIndex = 0;
650      beginCutoffGroupIndex = 0;
651 <    
652 < #else
648 <    
649 <    int nproc;
650 <    int myNode;
651 <    
652 <    myNode = worldRank;
653 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
654 <    
655 <    std::vector < int > tmpAtomsInProc(nproc, 0);
656 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
657 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
658 <    std::vector < int > NumAtomsInProc(nproc, 0);
659 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
660 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
661 <    
662 <    tmpAtomsInProc[myNode] = info->getNAtoms();
663 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
664 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
665 <    
666 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
667 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
668 <                  MPI_SUM, MPI_COMM_WORLD);
669 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
670 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
671 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
672 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
673 <    
674 <    beginAtomIndex = 0;
675 <    beginRigidBodyIndex = 0;
676 <    beginCutoffGroupIndex = 0;
677 <    
678 <    for(int i = 0; i < myNode; i++) {
679 <      beginAtomIndex += NumAtomsInProc[i];
680 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
681 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
682 <    }
683 <    
684 < #endif
685 <    
686 <    //rigidbody's index begins right after atom's
687 <    beginRigidBodyIndex += info->getNGlobalAtoms();
688 <    
689 <    for(mol = info->beginMolecule(mi); mol != NULL;
690 <        mol = info->nextMolecule(mi)) {
651 >
652 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
653        
654 <      //local index(index in DataStorge) of atom is important
655 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
656 <        atom->setGlobalIndex(beginAtomIndex++);
654 > #ifdef IS_MPI      
655 >      if (info->getMolToProc(i) == worldRank) {
656 > #endif        
657 >        // stuff to do if I own this molecule
658 >        mol = info->getMoleculeByGlobalIndex(i);
659 >
660 >        //local index(index in DataStorge) of atom is important
661 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
662 >          atom->setGlobalIndex(beginAtomIndex++);
663 >        }
664 >        
665 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
666 >            rb = mol->nextRigidBody(ri)) {
667 >          rb->setGlobalIndex(beginRigidBodyIndex++);
668 >        }
669 >        
670 >        //local index of cutoff group is trivial, it only depends on
671 >        //the order of travesing
672 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
673 >            cg = mol->nextCutoffGroup(ci)) {
674 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
675 >        }        
676 >        
677 > #ifdef IS_MPI        
678 >      }  else {
679 >
680 >        // stuff to do if I don't own this molecule
681 >        
682 >        int stampId = info->getMoleculeStampId(i);
683 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
684 >
685 >        beginAtomIndex += stamp->getNAtoms();
686 >        beginRigidBodyIndex += stamp->getNRigidBodies();
687 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
688        }
689 <      
690 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
691 <          rb = mol->nextRigidBody(ri)) {
692 <        rb->setGlobalIndex(beginRigidBodyIndex++);
700 <      }
701 <      
702 <      //local index of cutoff group is trivial, it only depends on the order of travesing
703 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
704 <          cg = mol->nextCutoffGroup(ci)) {
705 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
706 <      }
707 <    }
708 <    
689 > #endif          
690 >
691 >    } //end for(int i=0)  
692 >
693      //fill globalGroupMembership
694      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
695      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
# Line 717 | Line 701 | namespace oopse {
701          
702        }      
703      }
704 <    
704 >  
705   #ifdef IS_MPI    
706      // Since the globalGroupMembership has been zero filled and we've only
707      // poked values into the atoms we know, we can do an Allreduce
708      // to get the full globalGroupMembership array (We think).
709      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
710      // docs said we could.
711 <    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
711 >    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
712      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
713                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
714      info->setGlobalGroupMembership(tmpGroupMembership);
715 +
716 +    cerr << "ggm:\n";
717 +    for (int i = 0; i < tmpGroupMembership.size(); i++)
718 +      cerr << "i = " << i << "\t ggm(i) = " << tmpGroupMembership[i] << "\n";
719 +
720   #else
721      info->setGlobalGroupMembership(globalGroupMembership);
722   #endif
# Line 736 | Line 725 | namespace oopse {
725      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
726      
727      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
739      
728        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
729          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
730        }
731      }
732      
733   #ifdef IS_MPI
734 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
734 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
735      
736      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
737                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
# Line 769 | Line 757 | namespace oopse {
757      std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
758   #endif    
759  
760 < std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
761 <
762 < int startingIndex = 0;
763 < for (int i = 0; i < info->getNGlobalMolecules(); i++) {
764 <  startingIOIndexForMol[i] = startingIndex;
765 <  startingIndex += numIntegrableObjectsPerMol[i];
766 < }
767 <
768 < std::cerr << "nGIO = " << info->getNGlobalIntegrableObjects() << "\n";
769 < std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
782 < for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
760 >    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
761 >    
762 >    int startingIndex = 0;
763 >    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
764 >      startingIOIndexForMol[i] = startingIndex;
765 >      startingIndex += numIntegrableObjectsPerMol[i];
766 >    }
767 >    
768 >    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
769 >    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
770        int myGlobalIndex = mol->getGlobalIndex();
771        int globalIO = startingIOIndexForMol[myGlobalIndex];
785      std::cerr << "myGlobalIndex = " << myGlobalIndex << " globalIO = " << globalIO << "\n";
772        for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
773             integrableObject = mol->nextIntegrableObject(ioi)) {
774 <            integrableObject->setGlobalIntegrableObjectIndex(globalIO);
775 <            IOIndexToIntegrableObject[globalIO] = integrableObject;
776 <            globalIO++;
774 >        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
775 >        IOIndexToIntegrableObject[globalIO] = integrableObject;
776 >        globalIO++;
777        }
778      }
779 <
780 <  info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
781 <  
779 >    cerr << "ioi2io:\n";
780 >    for (int i = 0; i < IOIndexToIntegrableObject.size(); i++) {
781 >      if (IOIndexToIntegrableObject[i] != NULL) {
782 >        cerr << "i = " << i << "globalIOindex = " << IOIndexToIntegrableObject[i]->getGlobalIntegrableObjectIndex() << "\n";
783 >      }
784 >    }
785 >      
786 >    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
787 >    
788    }
789    
790    void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
791      Globals* simParams;
792 +
793      simParams = info->getSimParams();
794      
802    
795      DumpReader reader(info, mdFileName);
796      int nframes = reader.getNFrames();
797 <    
797 >
798      if (nframes > 0) {
799        reader.readFrame(nframes - 1);
800      } else {
# Line 813 | Line 805 | int startingIndex = 0;
805        painCave.isFatal = 1;
806        simError();
807      }
816    
808      //copy the current snapshot to previous snapshot
809      info->getSnapshotManager()->advance();
810    }
811    
812 < } //end namespace oopse
812 > } //end namespace OpenMD
813  
814  

Comparing:
trunk/src/brains/SimCreator.cpp (property svn:keywords), Revision 1024 by tim, Wed Aug 30 18:42:29 2006 UTC vs.
branches/development/src/brains/SimCreator.cpp (property svn:keywords), Revision 1600 by gezelter, Wed Aug 3 20:20:37 2011 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines