src/brains/SimCreator.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
/**
 * @file SimCreator.cpp
 * @author tlin
 * @date 11/03/2004
 * @time 13:51am
 * @version 1.0
 */

#include <iostream>
#include <sstream>
#include <string>

#include "brains/MoleculeCreator.hpp"
#include "brains/SimCreator.hpp"
#include "brains/SimSnapshotManager.hpp"
#include "io/DumpReader.hpp"
#include "UseTheForce/ForceFieldFactory.hpp"
#include "utils/simError.h"
#include "utils/StringUtils.hpp"
#include "math/SeqRandNumGen.hpp"
#include "mdParser/MDLexer.hpp"
#include "mdParser/MDParser.hpp"
#include "mdParser/MDTreeParser.hpp"
#include "mdParser/SimplePreprocessor.hpp"


#ifdef IS_MPI
#include "math/ParallelRandNumGen.hpp"
#endif

namespace oopse {
  
Globals* SimCreator::parseFile(const std::string mdFileName){
        Globals* simParams = NULL;
        try {

            // Create a preprocessor that preprocesses md file into an ostringstream
            std::stringstream ppStream;
#ifdef IS_MPI            
            int streamSize;
            const int masterNode = 0;
            int commStatus;
            if (worldRank == masterNode) {
#endif 
                
                SimplePreprocessor preprocessor;
                preprocessor.preprocess(mdFileName, ppStream);
                
#ifdef IS_MPI            
                //brocasting the stream size
                streamSize = ppStream.str().size() +1;
                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                   

                commStatus = MPI_Bcast(ppStream.str().c_str(), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); 
            
                
            } else {
                //get stream size
                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);   
                
                  char* buf = new char[streamSize];
                  assert(buf);
                
                  //receive file content
                  commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); 
                
                  ppStream.str(buf);
                  delete buf;

            }
#endif            
            // Create a scanner that reads from the input stream
            MDLexer lexer(ppStream);
            lexer.setFilename(mdFileName);
            lexer.initDeferredLineCount();
    
            // Create a parser that reads from the scanner
            MDParser parser(lexer);
            parser.setFilename(mdFileName);

            // Create an observer that synchorizes file name change
            FilenameObserver observer;
            observer.setLexer(&lexer);
            observer.setParser(&parser);
            lexer.setObserver(&observer);
    
            antlr::ASTFactory factory;
            parser.initializeASTFactory(factory);
            parser.setASTFactory(&factory);
            parser.mdfile();

            // Create a tree parser that reads information into Globals
            MDTreeParser treeParser;
            treeParser.initializeASTFactory(factory);
            treeParser.setASTFactory(&factory);
             simParams = treeParser.walkTree(parser.getAST());

        }
        catch (exception& e) {
            cerr << "parser exception: " << e.what() << endl;
        }

        return simParams;
  }
  
  SimInfo*  SimCreator::createSim(const std::string & mdFileName, 
                                  bool loadInitCoords) {

    //parse meta-data file
    Globals* simParams = parseFile(mdFileName);
    
    //create the force field
    ForceField * ff = ForceFieldFactory::getInstance()
      ->createForceField(simParams->getForceField());
    
    if (ff == NULL) {
      sprintf(painCave.errMsg, 
              "ForceField Factory can not create %s force field\n",
              simParams->getForceField().c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    if (simParams->haveForceFieldFileName()) {
      ff->setForceFieldFileName(simParams->getForceFieldFileName());
    }
    
    std::string forcefieldFileName;
    forcefieldFileName = ff->getForceFieldFileName();
    
    if (simParams->haveForceFieldVariant()) {
      //If the force field has variant, the variant force field name will be
      //Base.variant.frc. For exampel EAM.u6.frc
      
      std::string variant = simParams->getForceFieldVariant();
      
      std::string::size_type pos = forcefieldFileName.rfind(".frc");
      variant = "." + variant;
      if (pos != std::string::npos) {
        forcefieldFileName.insert(pos, variant);
      } else {
        //If the default force field file name does not containt .frc suffix, just append the .variant
        forcefieldFileName.append(variant);
      }
    } 
    
    ff->parse(forcefieldFileName);
        
    //create SimInfo
    SimInfo * info = new SimInfo(ff, simParams);
     
    //gather parameters (SimCreator only retrieves part of the parameters)
    gatherParameters(info, mdFileName);
    
    //divide the molecules and determine the global index of molecules
#ifdef IS_MPI
    divideMolecules(info);
#endif 
    
    //create the molecules
    createMolecules(info);
    
    
    //allocate memory for DataStorage(circular reference, need to break it)
    info->setSnapshotManager(new SimSnapshotManager(info));
    
    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
    //global index will never change again). Local indices of atoms and rigidbodies are already set by 
    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager. 
    setGlobalIndex(info);
    
    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
    //atoms don't have the global index yet  (their global index are all initialized to -1).
    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
    //we can determine the beginning global indices of atoms before they get created.
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
      info->addExcludePairs(mol);
    }
    
    if (loadInitCoords)
      loadCoordinates(info);    
    
    return info;
  }
  
  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
    
    //figure out the output file names
    std::string prefix;
    
#ifdef IS_MPI
    
    if (worldRank == 0) {
#endif // is_mpi
      Globals * simParams = info->getSimParams();
      if (simParams->haveFinalConfig()) {
        prefix = getPrefix(simParams->getFinalConfig());
      } else {
        prefix = getPrefix(mdfile);
      }
      
      info->setFinalConfigFileName(prefix + ".eor");
      info->setDumpFileName(prefix + ".dump");
      info->setStatFileName(prefix + ".stat");
      info->setRestFileName(prefix + ".zang");
      
#ifdef IS_MPI
      
    }
    
#endif
    
  }
  
#ifdef IS_MPI
  void SimCreator::divideMolecules(SimInfo *info) {
    double numerator;
    double denominator;
    double precast;
    double x;
    double y;
    double a;
    int old_atoms;
    int add_atoms;
    int new_atoms;
    int nTarget;
    int done;
    int i;
    int j;
    int loops;
    int which_proc;
    int nProcessors;
    std::vector<int> atomsPerProc;
    int nGlobalMols = info->getNGlobalMolecules();
    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
    
    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
    
    if (nProcessors > nGlobalMols) {
      sprintf(painCave.errMsg,
              "nProcessors (%d) > nMol (%d)\n"
              "\tThe number of processors is larger than\n"
              "\tthe number of molecules.  This will not result in a \n"
              "\tusable division of atoms for force decomposition.\n"
              "\tEither try a smaller number of processors, or run the\n"
              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
      
      painCave.isFatal = 1;
      simError();
    }
    
    int seedValue;
    Globals * simParams = info->getSimParams();
    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
    if (simParams->haveSeed()) {
      seedValue = simParams->getSeed();
      myRandom = new SeqRandNumGen(seedValue);
    }else {
      myRandom = new SeqRandNumGen();
    }   
    
    
    a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
    
    //initialize atomsPerProc
    atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
    
    if (worldRank == 0) {
      numerator = info->getNGlobalAtoms();
      denominator = nProcessors;
      precast = numerator / denominator;
      nTarget = (int)(precast + 0.5);
      
      for(i = 0; i < nGlobalMols; i++) {
        done = 0;
        loops = 0;
        
        while (!done) {
          loops++;
          
          // Pick a processor at random
          
          which_proc = (int) (myRandom->rand() * nProcessors);
          
          //get the molecule stamp first
          int stampId = info->getMoleculeStampId(i);
          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
          
          // How many atoms does this processor have so far?
          old_atoms = atomsPerProc[which_proc];
          add_atoms = moleculeStamp->getNAtoms();
          new_atoms = old_atoms + add_atoms;
          
          // If we've been through this loop too many times, we need
          // to just give up and assign the molecule to this processor
          // and be done with it. 
          
          if (loops > 100) {
            sprintf(painCave.errMsg,
                    "I've tried 100 times to assign molecule %d to a "
                    " processor, but can't find a good spot.\n"
                    "I'm assigning it at random to processor %d.\n",
                    i, which_proc);
            
            painCave.isFatal = 0;
            simError();
            
            molToProcMap[i] = which_proc;
            atomsPerProc[which_proc] += add_atoms;
            
            done = 1;
            continue;
          }
          
          // If we can add this molecule to this processor without sending
          // it above nTarget, then go ahead and do it:
          
          if (new_atoms <= nTarget) {
            molToProcMap[i] = which_proc;
            atomsPerProc[which_proc] += add_atoms;
            
            done = 1;
            continue;
          }
          
          // The only situation left is when new_atoms > nTarget.  We
          // want to accept this with some probability that dies off the
          // farther we are from nTarget
          
          // roughly:  x = new_atoms - nTarget
          //           Pacc(x) = exp(- a * x)
          // where a = penalty / (average atoms per molecule)
          
          x = (double)(new_atoms - nTarget);
          y = myRandom->rand();
          
          if (y < exp(- a * x)) {
            molToProcMap[i] = which_proc;
            atomsPerProc[which_proc] += add_atoms;
            
            done = 1;
            continue;
          } else {
            continue;
          }
        }
      }
      
      delete myRandom;
      
      // Spray out this nonsense to all other processors:
      
      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    } else {
      
      // Listen to your marching orders from processor 0:
      
      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    }
    
    info->setMolToProcMap(molToProcMap);
    sprintf(checkPointMsg,
            "Successfully divided the molecules among the processors.\n");
    MPIcheckPoint();
  }
  
#endif
  
  void SimCreator::createMolecules(SimInfo *info) {
    MoleculeCreator molCreator;
    int stampId;
    
    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
      
#ifdef IS_MPI
      
      if (info->getMolToProc(i) == worldRank) {
#endif
        
        stampId = info->getMoleculeStampId(i);
        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
                                                   stampId, i, info->getLocalIndexManager());
        
        info->addMolecule(mol);
        
#ifdef IS_MPI
        
      }
      
#endif
      
    } //end for(int i=0)   
  }
    
  void SimCreator::setGlobalIndex(SimInfo *info) {
    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator ri;
    Molecule::CutoffGroupIterator ci;
    Molecule * mol;
    Atom * atom;
    RigidBody * rb;
    CutoffGroup * cg;
    int beginAtomIndex;
    int beginRigidBodyIndex;
    int beginCutoffGroupIndex;
    int nGlobalAtoms = info->getNGlobalAtoms();
    
#ifndef IS_MPI
    
    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;
    
#else
    
    int nproc;
    int myNode;
    
    myNode = worldRank;
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    
    std::vector < int > tmpAtomsInProc(nproc, 0);
    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
    std::vector < int > NumAtomsInProc(nproc, 0);
    std::vector < int > NumRigidBodiesInProc(nproc, 0);
    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
    
    tmpAtomsInProc[myNode] = info->getNAtoms();
    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
    
    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
                  MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    
    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;
    
    for(int i = 0; i < myNode; i++) {
      beginAtomIndex += NumAtomsInProc[i];
      beginRigidBodyIndex += NumRigidBodiesInProc[i];
      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
    }
    
#endif
    
    //rigidbody's index begins right after atom's
    beginRigidBodyIndex += info->getNGlobalAtoms();
    
    for(mol = info->beginMolecule(mi); mol != NULL;
        mol = info->nextMolecule(mi)) {
      
      //local index(index in DataStorge) of atom is important
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
        atom->setGlobalIndex(beginAtomIndex++);
      }
      
      for(rb = mol->beginRigidBody(ri); rb != NULL;
          rb = mol->nextRigidBody(ri)) {
        rb->setGlobalIndex(beginRigidBodyIndex++);
      }
      
      //local index of cutoff group is trivial, it only depends on the order of travesing
      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
          cg = mol->nextCutoffGroup(ci)) {
        cg->setGlobalIndex(beginCutoffGroupIndex++);
      }
    }
    
    //fill globalGroupMembership
    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
        
        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
        }
        
      }       
    }
    
#ifdef IS_MPI    
    // Since the globalGroupMembership has been zero filled and we've only
    // poked values into the atoms we know, we can do an Allreduce
    // to get the full globalGroupMembership array (We think).
    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
    // docs said we could.
    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    info->setGlobalGroupMembership(tmpGroupMembership);
#else
    info->setGlobalGroupMembership(globalGroupMembership);
#endif
    
    //fill molMembership
    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
    
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
      
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
      }
    }
    
#ifdef IS_MPI
    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
    
    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    
    info->setGlobalMolMembership(tmpMolMembership);
#else
    info->setGlobalMolMembership(globalMolMembership);
#endif
    
  }
  
  void SimCreator::loadCoordinates(SimInfo* info) {
    Globals* simParams;
    simParams = info->getSimParams();
    
    if (!simParams->haveInitialConfig()) {
      sprintf(painCave.errMsg,
              "Cannot intialize a simulation without an initial configuration file.\n");
      painCave.isFatal = 1;;
      simError();
    }
    
    DumpReader reader(info, simParams->getInitialConfig());
    int nframes = reader.getNFrames();
    
    if (nframes > 0) {
      reader.readFrame(nframes - 1);
    } else {
      //invalid initial coordinate file
      sprintf(painCave.errMsg, 
              "Initial configuration file %s should at least contain one frame\n",
              simParams->getInitialConfig().c_str());
      painCave.isFatal = 1;
      simError();
    }
    
    //copy the current snapshot to previous snapshot
    info->getSnapshotManager()->advance();
  }
  
} //end namespace oopse


Revision:	770
Committed:	Fri Dec 2 15:38:03 2005 UTC (19 years, 5 months ago) by tim
File size:	19946 byte(s)
Log Message:	End of the Link --> List Return of the Oject-Oriented replace yacc/lex parser with antlr parser
#	User	Rev	Content
1	gezelter	403	/*
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
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
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42			/**
43			* @file SimCreator.cpp
44			* @author tlin
45			* @date 11/03/2004
46			* @time 13:51am
47			* @version 1.0
48			*/
49
50	tim	770	#include <iostream>
51			#include <sstream>
52			#include <string>
53
54	gezelter	403	#include "brains/MoleculeCreator.hpp"
55			#include "brains/SimCreator.hpp"
56			#include "brains/SimSnapshotManager.hpp"
57			#include "io/DumpReader.hpp"
58			#include "UseTheForce/ForceFieldFactory.hpp"
59			#include "utils/simError.h"
60			#include "utils/StringUtils.hpp"
61			#include "math/SeqRandNumGen.hpp"
62	tim	770	#include "mdParser/MDLexer.hpp"
63			#include "mdParser/MDParser.hpp"
64			#include "mdParser/MDTreeParser.hpp"
65			#include "mdParser/SimplePreprocessor.hpp"
66
67
68	gezelter	403	#ifdef IS_MPI
69			#include "math/ParallelRandNumGen.hpp"
70			#endif
71
72			namespace oopse {
73	chrisfen	417
74	tim	770	Globals* SimCreator::parseFile(const std::string mdFileName){
75			Globals* simParams = NULL;
76			try {
77
78			// Create a preprocessor that preprocesses md file into an ostringstream
79			std::stringstream ppStream;
80			#ifdef IS_MPI
81			int streamSize;
82			const int masterNode = 0;
83			int commStatus;
84			if (worldRank == masterNode) {
85			#endif
86
87			SimplePreprocessor preprocessor;
88			preprocessor.preprocess(mdFileName, ppStream);
89
90			#ifdef IS_MPI
91			//brocasting the stream size
92			streamSize = ppStream.str().size() +1;
93			commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);
94
95			commStatus = MPI_Bcast(ppStream.str().c_str(), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
96
97
98			} else {
99			//get stream size
100			commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);
101
102			char* buf = new char[streamSize];
103			assert(buf);
104
105			//receive file content
106			commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
107
108			ppStream.str(buf);
109			delete buf;
110
111			}
112			#endif
113			// Create a scanner that reads from the input stream
114			MDLexer lexer(ppStream);
115			lexer.setFilename(mdFileName);
116			lexer.initDeferredLineCount();
117	chrisfen	417
118	tim	770	// Create a parser that reads from the scanner
119			MDParser parser(lexer);
120			parser.setFilename(mdFileName);
121
122			// Create an observer that synchorizes file name change
123			FilenameObserver observer;
124			observer.setLexer(&lexer);
125			observer.setParser(&parser);
126			lexer.setObserver(&observer);
127	chrisfen	417
128	tim	770	antlr::ASTFactory factory;
129			parser.initializeASTFactory(factory);
130			parser.setASTFactory(&factory);
131			parser.mdfile();
132
133			// Create a tree parser that reads information into Globals
134			MDTreeParser treeParser;
135			treeParser.initializeASTFactory(factory);
136			treeParser.setASTFactory(&factory);
137			simParams = treeParser.walkTree(parser.getAST());
138
139			}
140			catch (exception& e) {
141			cerr << "parser exception: " << e.what() << endl;
142			}
143
144			return simParams;
145	gezelter	403	}
146	chrisfen	417
147	chrisfen	514	SimInfo* SimCreator::createSim(const std::string & mdFileName,
148			bool loadInitCoords) {
149	tim	770
150	gezelter	403	//parse meta-data file
151	tim	770	Globals* simParams = parseFile(mdFileName);
152	chrisfen	417
153	gezelter	403	//create the force field
154	chrisfen	514	ForceField * ff = ForceFieldFactory::getInstance()
155			->createForceField(simParams->getForceField());
156	gezelter	403
157			if (ff == NULL) {
158	chrisfen	514	sprintf(painCave.errMsg,
159			"ForceField Factory can not create %s force field\n",
160	tim	665	simParams->getForceField().c_str());
161	gezelter	403	painCave.isFatal = 1;
162			simError();
163			}
164	chrisfen	417
165	gezelter	403	if (simParams->haveForceFieldFileName()) {
166			ff->setForceFieldFileName(simParams->getForceFieldFileName());
167			}
168
169			std::string forcefieldFileName;
170			forcefieldFileName = ff->getForceFieldFileName();
171	chrisfen	417
172	gezelter	403	if (simParams->haveForceFieldVariant()) {
173			//If the force field has variant, the variant force field name will be
174			//Base.variant.frc. For exampel EAM.u6.frc
175	chrisfen	417
176	gezelter	403	std::string variant = simParams->getForceFieldVariant();
177	chrisfen	417
178	gezelter	403	std::string::size_type pos = forcefieldFileName.rfind(".frc");
179			variant = "." + variant;
180			if (pos != std::string::npos) {
181			forcefieldFileName.insert(pos, variant);
182			} else {
183			//If the default force field file name does not containt .frc suffix, just append the .variant
184			forcefieldFileName.append(variant);
185			}
186			}
187
188			ff->parse(forcefieldFileName);
189	tim	770
190	gezelter	403	//create SimInfo
191	tim	770	SimInfo * info = new SimInfo(ff, simParams);
192	tim	490
193	gezelter	403	//gather parameters (SimCreator only retrieves part of the parameters)
194			gatherParameters(info, mdFileName);
195	chrisfen	417
196	gezelter	403	//divide the molecules and determine the global index of molecules
197			#ifdef IS_MPI
198			divideMolecules(info);
199			#endif
200	chrisfen	417
201	gezelter	403	//create the molecules
202			createMolecules(info);
203	chrisfen	417
204
205	gezelter	403	//allocate memory for DataStorage(circular reference, need to break it)
206			info->setSnapshotManager(new SimSnapshotManager(info));
207
208			//set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
209			//global index will never change again). Local indices of atoms and rigidbodies are already set by
210			//MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
211			setGlobalIndex(info);
212	chrisfen	417
213	gezelter	403	//Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
214			//atoms don't have the global index yet (their global index are all initialized to -1).
215			//Therefore we have to call addExcludePairs explicitly here. A way to work around is that
216			//we can determine the beginning global indices of atoms before they get created.
217			SimInfo::MoleculeIterator mi;
218			Molecule* mol;
219			for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
220			info->addExcludePairs(mol);
221			}
222
223			if (loadInitCoords)
224			loadCoordinates(info);
225
226			return info;
227			}
228	chrisfen	417
229	gezelter	403	void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
230	chrisfen	417
231	tim	749	//figure out the output file names
232	gezelter	403	std::string prefix;
233	chrisfen	417
234	gezelter	403	#ifdef IS_MPI
235	chrisfen	417
236	gezelter	403	if (worldRank == 0) {
237			#endif // is_mpi
238			Globals * simParams = info->getSimParams();
239			if (simParams->haveFinalConfig()) {
240			prefix = getPrefix(simParams->getFinalConfig());
241			} else {
242			prefix = getPrefix(mdfile);
243			}
244	chrisfen	417
245	gezelter	403	info->setFinalConfigFileName(prefix + ".eor");
246			info->setDumpFileName(prefix + ".dump");
247			info->setStatFileName(prefix + ".stat");
248	chrisfen	417	info->setRestFileName(prefix + ".zang");
249
250	gezelter	403	#ifdef IS_MPI
251	chrisfen	417
252	gezelter	403	}
253	chrisfen	417
254	gezelter	403	#endif
255	chrisfen	417
256	gezelter	403	}
257	chrisfen	417
258	gezelter	403	#ifdef IS_MPI
259			void SimCreator::divideMolecules(SimInfo *info) {
260			double numerator;
261			double denominator;
262			double precast;
263			double x;
264			double y;
265			double a;
266			int old_atoms;
267			int add_atoms;
268			int new_atoms;
269			int nTarget;
270			int done;
271			int i;
272			int j;
273			int loops;
274			int which_proc;
275			int nProcessors;
276			std::vector<int> atomsPerProc;
277			int nGlobalMols = info->getNGlobalMolecules();
278			std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
279
280			MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
281	chrisfen	417
282	gezelter	403	if (nProcessors > nGlobalMols) {
283			sprintf(painCave.errMsg,
284			"nProcessors (%d) > nMol (%d)\n"
285			"\tThe number of processors is larger than\n"
286			"\tthe number of molecules. This will not result in a \n"
287			"\tusable division of atoms for force decomposition.\n"
288			"\tEither try a smaller number of processors, or run the\n"
289			"\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
290	chrisfen	417
291	gezelter	403	painCave.isFatal = 1;
292			simError();
293			}
294	chrisfen	417
295	gezelter	403	int seedValue;
296			Globals * simParams = info->getSimParams();
297			SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
298			if (simParams->haveSeed()) {
299			seedValue = simParams->getSeed();
300			myRandom = new SeqRandNumGen(seedValue);
301			}else {
302			myRandom = new SeqRandNumGen();
303			}
304	chrisfen	417
305
306	gezelter	403	a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
307	chrisfen	417
308	gezelter	403	//initialize atomsPerProc
309			atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
310	chrisfen	417
311	gezelter	403	if (worldRank == 0) {
312			numerator = info->getNGlobalAtoms();
313			denominator = nProcessors;
314			precast = numerator / denominator;
315			nTarget = (int)(precast + 0.5);
316	chrisfen	417
317	gezelter	403	for(i = 0; i < nGlobalMols; i++) {
318			done = 0;
319			loops = 0;
320	chrisfen	417
321	gezelter	403	while (!done) {
322			loops++;
323	chrisfen	417
324	gezelter	403	// Pick a processor at random
325	chrisfen	417
326	gezelter	403	which_proc = (int) (myRandom->rand() * nProcessors);
327	chrisfen	417
328	gezelter	403	//get the molecule stamp first
329			int stampId = info->getMoleculeStampId(i);
330			MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
331	chrisfen	417
332	gezelter	403	// How many atoms does this processor have so far?
333			old_atoms = atomsPerProc[which_proc];
334			add_atoms = moleculeStamp->getNAtoms();
335			new_atoms = old_atoms + add_atoms;
336	chrisfen	417
337	gezelter	403	// If we've been through this loop too many times, we need
338			// to just give up and assign the molecule to this processor
339			// and be done with it.
340	chrisfen	417
341	gezelter	403	if (loops > 100) {
342			sprintf(painCave.errMsg,
343			"I've tried 100 times to assign molecule %d to a "
344			" processor, but can't find a good spot.\n"
345			"I'm assigning it at random to processor %d.\n",
346			i, which_proc);
347	chrisfen	417
348	gezelter	403	painCave.isFatal = 0;
349			simError();
350	chrisfen	417
351	gezelter	403	molToProcMap[i] = which_proc;
352			atomsPerProc[which_proc] += add_atoms;
353	chrisfen	417
354	gezelter	403	done = 1;
355			continue;
356			}
357	chrisfen	417
358	gezelter	403	// If we can add this molecule to this processor without sending
359			// it above nTarget, then go ahead and do it:
360	chrisfen	417
361	gezelter	403	if (new_atoms <= nTarget) {
362			molToProcMap[i] = which_proc;
363			atomsPerProc[which_proc] += add_atoms;
364	chrisfen	417
365	gezelter	403	done = 1;
366			continue;
367			}
368	chrisfen	417
369	gezelter	403	// The only situation left is when new_atoms > nTarget. We
370			// want to accept this with some probability that dies off the
371			// farther we are from nTarget
372	chrisfen	417
373	gezelter	403	// roughly: x = new_atoms - nTarget
374			// Pacc(x) = exp(- a * x)
375			// where a = penalty / (average atoms per molecule)
376	chrisfen	417
377	gezelter	403	x = (double)(new_atoms - nTarget);
378			y = myRandom->rand();
379	chrisfen	417
380	gezelter	403	if (y < exp(- a * x)) {
381			molToProcMap[i] = which_proc;
382			atomsPerProc[which_proc] += add_atoms;
383	chrisfen	417
384	gezelter	403	done = 1;
385			continue;
386			} else {
387			continue;
388			}
389			}
390			}
391	chrisfen	417
392	gezelter	403	delete myRandom;
393	chrisfen	417
394	gezelter	403	// Spray out this nonsense to all other processors:
395	chrisfen	417
396	gezelter	403	MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
397			} else {
398	chrisfen	417
399	gezelter	403	// Listen to your marching orders from processor 0:
400	chrisfen	417
401	gezelter	403	MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
402			}
403	chrisfen	417
404	gezelter	403	info->setMolToProcMap(molToProcMap);
405			sprintf(checkPointMsg,
406			"Successfully divided the molecules among the processors.\n");
407			MPIcheckPoint();
408			}
409	chrisfen	417
410	gezelter	403	#endif
411	chrisfen	417
412	gezelter	403	void SimCreator::createMolecules(SimInfo *info) {
413			MoleculeCreator molCreator;
414			int stampId;
415	chrisfen	417
416	gezelter	403	for(int i = 0; i < info->getNGlobalMolecules(); i++) {
417	chrisfen	417
418	gezelter	403	#ifdef IS_MPI
419	chrisfen	417
420	gezelter	403	if (info->getMolToProc(i) == worldRank) {
421			#endif
422	chrisfen	417
423	gezelter	403	stampId = info->getMoleculeStampId(i);
424			Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
425			stampId, i, info->getLocalIndexManager());
426	chrisfen	417
427	gezelter	403	info->addMolecule(mol);
428	chrisfen	417
429	gezelter	403	#ifdef IS_MPI
430	chrisfen	417
431	gezelter	403	}
432	chrisfen	417
433	gezelter	403	#endif
434	chrisfen	417
435	gezelter	403	} //end for(int i=0)
436			}
437	chrisfen	417
438	gezelter	403	void SimCreator::setGlobalIndex(SimInfo *info) {
439			SimInfo::MoleculeIterator mi;
440			Molecule::AtomIterator ai;
441			Molecule::RigidBodyIterator ri;
442			Molecule::CutoffGroupIterator ci;
443			Molecule * mol;
444			Atom * atom;
445			RigidBody * rb;
446			CutoffGroup * cg;
447			int beginAtomIndex;
448			int beginRigidBodyIndex;
449			int beginCutoffGroupIndex;
450			int nGlobalAtoms = info->getNGlobalAtoms();
451
452			#ifndef IS_MPI
453	chrisfen	417
454	gezelter	403	beginAtomIndex = 0;
455			beginRigidBodyIndex = 0;
456			beginCutoffGroupIndex = 0;
457	chrisfen	417
458	gezelter	403	#else
459	chrisfen	417
460	gezelter	403	int nproc;
461			int myNode;
462	chrisfen	417
463	gezelter	403	myNode = worldRank;
464			MPI_Comm_size(MPI_COMM_WORLD, &nproc);
465	chrisfen	417
466	gezelter	403	std::vector < int > tmpAtomsInProc(nproc, 0);
467			std::vector < int > tmpRigidBodiesInProc(nproc, 0);
468			std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
469			std::vector < int > NumAtomsInProc(nproc, 0);
470			std::vector < int > NumRigidBodiesInProc(nproc, 0);
471			std::vector < int > NumCutoffGroupsInProc(nproc, 0);
472	chrisfen	417
473	gezelter	403	tmpAtomsInProc[myNode] = info->getNAtoms();
474			tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
475			tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
476	chrisfen	417
477	gezelter	403	//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
478			MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
479			MPI_SUM, MPI_COMM_WORLD);
480			MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
481			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
482			MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
483			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
484	chrisfen	417
485	gezelter	403	beginAtomIndex = 0;
486			beginRigidBodyIndex = 0;
487			beginCutoffGroupIndex = 0;
488	chrisfen	417
489	gezelter	403	for(int i = 0; i < myNode; i++) {
490			beginAtomIndex += NumAtomsInProc[i];
491			beginRigidBodyIndex += NumRigidBodiesInProc[i];
492			beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
493			}
494	chrisfen	417
495	gezelter	403	#endif
496	chrisfen	417
497	gezelter	403	//rigidbody's index begins right after atom's
498			beginRigidBodyIndex += info->getNGlobalAtoms();
499	chrisfen	417
500	gezelter	403	for(mol = info->beginMolecule(mi); mol != NULL;
501			mol = info->nextMolecule(mi)) {
502	chrisfen	417
503	gezelter	403	//local index(index in DataStorge) of atom is important
504			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
505			atom->setGlobalIndex(beginAtomIndex++);
506			}
507	chrisfen	417
508	gezelter	403	for(rb = mol->beginRigidBody(ri); rb != NULL;
509			rb = mol->nextRigidBody(ri)) {
510			rb->setGlobalIndex(beginRigidBodyIndex++);
511			}
512	chrisfen	417
513	gezelter	403	//local index of cutoff group is trivial, it only depends on the order of travesing
514			for(cg = mol->beginCutoffGroup(ci); cg != NULL;
515			cg = mol->nextCutoffGroup(ci)) {
516			cg->setGlobalIndex(beginCutoffGroupIndex++);
517			}
518			}
519	chrisfen	417
520	gezelter	403	//fill globalGroupMembership
521			std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
522			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
523			for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
524	chrisfen	417
525	gezelter	403	for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
526			globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
527			}
528	chrisfen	417
529	gezelter	403	}
530			}
531	chrisfen	417
532	gezelter	403	#ifdef IS_MPI
533			// Since the globalGroupMembership has been zero filled and we've only
534			// poked values into the atoms we know, we can do an Allreduce
535			// to get the full globalGroupMembership array (We think).
536			// This would be prettier if we could use MPI_IN_PLACE like the MPI-2
537			// docs said we could.
538			std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
539			MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
540			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
541			info->setGlobalGroupMembership(tmpGroupMembership);
542			#else
543			info->setGlobalGroupMembership(globalGroupMembership);
544			#endif
545	chrisfen	417
546	gezelter	403	//fill molMembership
547			std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
548
549			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
550	chrisfen	417
551	gezelter	403	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
552			globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
553			}
554			}
555	chrisfen	417
556	gezelter	403	#ifdef IS_MPI
557			std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
558	chrisfen	417
559	gezelter	403	MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
560			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
561
562			info->setGlobalMolMembership(tmpMolMembership);
563			#else
564			info->setGlobalMolMembership(globalMolMembership);
565			#endif
566	chrisfen	417
567	gezelter	403	}
568	chrisfen	417
569	gezelter	403	void SimCreator::loadCoordinates(SimInfo* info) {
570			Globals* simParams;
571			simParams = info->getSimParams();
572
573			if (!simParams->haveInitialConfig()) {
574			sprintf(painCave.errMsg,
575			"Cannot intialize a simulation without an initial configuration file.\n");
576			painCave.isFatal = 1;;
577			simError();
578			}
579	chrisfen	417
580	gezelter	403	DumpReader reader(info, simParams->getInitialConfig());
581			int nframes = reader.getNFrames();
582	chrisfen	417
583	gezelter	403	if (nframes > 0) {
584			reader.readFrame(nframes - 1);
585			} else {
586			//invalid initial coordinate file
587	chrisfen	417	sprintf(painCave.errMsg,
588			"Initial configuration file %s should at least contain one frame\n",
589	tim	665	simParams->getInitialConfig().c_str());
590	gezelter	403	painCave.isFatal = 1;
591			simError();
592			}
593	chrisfen	417
594	gezelter	403	//copy the current snapshot to previous snapshot
595			info->getSnapshotManager()->advance();
596			}
597	chrisfen	417
598	gezelter	403	} //end namespace oopse
599
600