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 417 by chrisfen, Thu Mar 10 15:10:24 2005 UTC vs.
Revision 1136 by chuckv, Sat May 26 17:53:52 2007 UTC

# Line 46 | Line 46
46   * @time 13:51am
47   * @version 1.0
48   */
49 + #include <exception>
50 + #include <iostream>
51 + #include <sstream>
52 + #include <string>
53  
54   #include "brains/MoleculeCreator.hpp"
55   #include "brains/SimCreator.hpp"
56   #include "brains/SimSnapshotManager.hpp"
57   #include "io/DumpReader.hpp"
54 – #include "io/parse_me.h"
58   #include "UseTheForce/ForceFieldFactory.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61   #include "math/SeqRandNumGen.hpp"
62 + #include "mdParser/MDLexer.hpp"
63 + #include "mdParser/MDParser.hpp"
64 + #include "mdParser/MDTreeParser.hpp"
65 + #include "mdParser/SimplePreprocessor.hpp"
66 + #include "antlr/ANTLRException.hpp"
67 + #include "antlr/TokenStreamRecognitionException.hpp"
68 + #include "antlr/TokenStreamIOException.hpp"
69 + #include "antlr/TokenStreamException.hpp"
70 + #include "antlr/RecognitionException.hpp"
71 + #include "antlr/CharStreamException.hpp"
72 +
73 + #include "antlr/MismatchedCharException.hpp"
74 + #include "antlr/MismatchedTokenException.hpp"
75 + #include "antlr/NoViableAltForCharException.hpp"
76 + #include "antlr/NoViableAltException.hpp"
77 +
78   #ifdef IS_MPI
60 – #include "io/mpiBASS.h"
79   #include "math/ParallelRandNumGen.hpp"
80   #endif
81  
82   namespace oopse {
83    
84 <  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps,
85 <                             Globals* simParams){
84 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
85 >    Globals* simParams = NULL;
86 >    try {
87 >
88 >      // Create a preprocessor that preprocesses md file into an ostringstream
89 >      std::stringstream ppStream;
90 > #ifdef IS_MPI            
91 >      int streamSize;
92 >      const int masterNode = 0;
93 >      int commStatus;
94 >      if (worldRank == masterNode) {
95 > #endif
96 >                
97 >        SimplePreprocessor preprocessor;
98 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
99 >                
100 > #ifdef IS_MPI            
101 >        //brocasting the stream size
102 >        streamSize = ppStream.str().size() +1;
103 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
104 >
105 >        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
106 >            
107 >                
108 >      } else {
109 >        //get stream size
110 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
111 >                
112 >        char* buf = new char[streamSize];
113 >        assert(buf);
114 >                
115 >        //receive file content
116 >        commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
117 >                
118 >        ppStream.str(buf);
119 >        delete buf;
120 >
121 >      }
122 > #endif            
123 >      // Create a scanner that reads from the input stream
124 >      MDLexer lexer(ppStream);
125 >      lexer.setFilename(filename);
126 >      lexer.initDeferredLineCount();
127      
128 < #ifdef IS_MPI
128 >      // Create a parser that reads from the scanner
129 >      MDParser parser(lexer);
130 >      parser.setFilename(filename);
131 >
132 >      // Create an observer that synchorizes file name change
133 >      FilenameObserver observer;
134 >      observer.setLexer(&lexer);
135 >      observer.setParser(&parser);
136 >      lexer.setObserver(&observer);
137      
138 <    if (worldRank == 0) {
139 < #endif // is_mpi
138 >      antlr::ASTFactory factory;
139 >      parser.initializeASTFactory(factory);
140 >      parser.setASTFactory(&factory);
141 >      parser.mdfile();
142 >
143 >      // Create a tree parser that reads information into Globals
144 >      MDTreeParser treeParser;
145 >      treeParser.initializeASTFactory(factory);
146 >      treeParser.setASTFactory(&factory);
147 >      simParams = treeParser.walkTree(parser.getAST());
148 >
149 >    }
150 >
151        
152 <      simParams->initalize();
153 <      set_interface_stamps(stamps, simParams);
152 >    catch(antlr::MismatchedCharException& e) {
153 >      sprintf(painCave.errMsg,
154 >              "parser exception: %s %s:%d:%d\n",
155 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
156 >      painCave.isFatal = 1;
157 >      simError();          
158 >    }
159 >    catch(antlr::MismatchedTokenException &e) {
160 >      sprintf(painCave.errMsg,
161 >              "parser exception: %s %s:%d:%d\n",
162 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
163 >      painCave.isFatal = 1;
164 >      simError();  
165 >    }
166 >    catch(antlr::NoViableAltForCharException &e) {
167 >      sprintf(painCave.errMsg,
168 >              "parser exception: %s %s:%d:%d\n",
169 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
170 >      painCave.isFatal = 1;
171 >      simError();  
172 >    }
173 >    catch(antlr::NoViableAltException &e) {
174 >      sprintf(painCave.errMsg,
175 >              "parser exception: %s %s:%d:%d\n",
176 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
177 >      painCave.isFatal = 1;
178 >      simError();  
179 >    }
180        
181 < #ifdef IS_MPI
182 <      
183 <      mpiEventInit();
184 <      
185 < #endif
186 <      
83 <      yacc_BASS(mdFileName.c_str());
84 <      
85 < #ifdef IS_MPI
86 <      
87 <      throwMPIEvent(NULL);
88 <    } else {
89 <      set_interface_stamps(stamps, simParams);
90 <      mpiEventInit();
91 <      MPIcheckPoint();
92 <      mpiEventLoop();
181 >    catch(antlr::TokenStreamRecognitionException& e) {
182 >      sprintf(painCave.errMsg,
183 >              "parser exception: %s %s:%d:%d\n",
184 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
185 >      painCave.isFatal = 1;
186 >      simError();  
187      }
188 <    
189 < #endif
190 <    
188 >        
189 >    catch(antlr::TokenStreamIOException& e) {
190 >      sprintf(painCave.errMsg,
191 >              "parser exception: %s\n",
192 >              e.getMessage().c_str());
193 >      painCave.isFatal = 1;
194 >      simError();
195 >    }
196 >        
197 >    catch(antlr::TokenStreamException& e) {
198 >      sprintf(painCave.errMsg,
199 >              "parser exception: %s\n",
200 >              e.getMessage().c_str());
201 >      painCave.isFatal = 1;
202 >      simError();
203 >    }        
204 >    catch (antlr::RecognitionException& e) {
205 >      sprintf(painCave.errMsg,
206 >              "parser exception: %s %s:%d:%d\n",
207 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
208 >      painCave.isFatal = 1;
209 >      simError();          
210 >    }
211 >    catch (antlr::CharStreamException& e) {
212 >      sprintf(painCave.errMsg,
213 >              "parser exception: %s\n",
214 >              e.getMessage().c_str());
215 >      painCave.isFatal = 1;
216 >      simError();        
217 >    }
218 >    catch (OOPSEException& e) {
219 >      sprintf(painCave.errMsg,
220 >              "%s\n",
221 >              e.getMessage().c_str());
222 >      painCave.isFatal = 1;
223 >      simError();
224 >    }
225 >    catch (std::exception& e) {
226 >      sprintf(painCave.errMsg,
227 >              "parser exception: %s\n",
228 >              e.what());
229 >      painCave.isFatal = 1;
230 >      simError();
231 >    }
232 >
233 >    return simParams;
234    }
235    
236 <  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
237 <    
238 <    MakeStamps * stamps = new MakeStamps();
239 <    
240 <    Globals * simParams = new Globals();
241 <    
236 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
237 >                                  bool loadInitCoords) {
238 >
239 >    const int bufferSize = 65535;
240 >    char buffer[bufferSize];
241 >    int lineNo = 0;
242 >    std::string mdRawData;
243 >    int metaDataBlockStart = -1;
244 >    int metaDataBlockEnd = -1;
245 >    int i;
246 >    int mdOffset;
247 >
248 > #ifdef IS_MPI            
249 >    const int masterNode = 0;
250 >    if (worldRank == masterNode) {
251 > #endif
252 >
253 >      std::ifstream mdFile_(mdFileName.c_str());
254 >      
255 >      if (mdFile_.fail()) {
256 >        sprintf(painCave.errMsg,
257 >                "SimCreator: Cannot open file: %s\n",
258 >                mdFileName.c_str());
259 >        painCave.isFatal = 1;
260 >        simError();
261 >      }
262 >
263 >      mdFile_.getline(buffer, bufferSize);
264 >      ++lineNo;
265 >      std::string line = trimLeftCopy(buffer);
266 >      i = CaseInsensitiveFind(line, "<OOPSE");
267 >      if (i == string::npos) {
268 >        sprintf(painCave.errMsg,
269 >                "SimCreator: File: %s is not an OOPSE file!\n",
270 >                mdFileName.c_str());
271 >        painCave.isFatal = 1;
272 >        simError();
273 >      }
274 >
275 >      //scan through the input stream and find MetaData tag        
276 >      while(mdFile_.getline(buffer, bufferSize)) {
277 >        ++lineNo;
278 >        
279 >        std::string line = trimLeftCopy(buffer);
280 >        if (metaDataBlockStart == -1) {
281 >          i = CaseInsensitiveFind(line, "<MetaData>");
282 >          if (i != string::npos) {
283 >            metaDataBlockStart = lineNo;
284 >            mdOffset = mdFile_.tellg();
285 >          }
286 >        } else {
287 >          i = CaseInsensitiveFind(line, "</MetaData>");
288 >          if (i != string::npos) {
289 >            metaDataBlockEnd = lineNo;
290 >          }
291 >        }
292 >      }
293 >
294 >      if (metaDataBlockStart == -1) {
295 >        sprintf(painCave.errMsg,
296 >                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
297 >                mdFileName.c_str());
298 >        painCave.isFatal = 1;
299 >        simError();
300 >      }
301 >      if (metaDataBlockEnd == -1) {
302 >        sprintf(painCave.errMsg,
303 >                "SimCreator: File: %s did not contain a closed MetaData block!\n",
304 >                mdFileName.c_str());
305 >        painCave.isFatal = 1;
306 >        simError();
307 >      }
308 >        
309 >      mdFile_.clear();
310 >      mdFile_.seekg(0);
311 >      mdFile_.seekg(mdOffset);
312 >
313 >      mdRawData.clear();
314 >
315 >      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
316 >        mdFile_.getline(buffer, bufferSize);
317 >        mdRawData += buffer;
318 >        mdRawData += "\n";
319 >      }
320 >
321 >      mdFile_.close();
322 >
323 > #ifdef IS_MPI
324 >    }
325 > #endif
326 >
327 >    std::stringstream rawMetaDataStream(mdRawData);
328 >
329      //parse meta-data file
330 <    parseFile(mdFileName, stamps, simParams);
330 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
331      
332      //create the force field
333 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
334 <                                                                         simParams->getForceField());
333 >    ForceField * ff = ForceFieldFactory::getInstance()
334 >      ->createForceField(simParams->getForceField());
335      
336      if (ff == NULL) {
337 <      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
338 <              simParams->getForceField());
337 >      sprintf(painCave.errMsg,
338 >              "ForceField Factory can not create %s force field\n",
339 >              simParams->getForceField().c_str());
340        painCave.isFatal = 1;
341        simError();
342      }
# Line 140 | Line 365 | namespace oopse {
365      }
366      
367      ff->parse(forcefieldFileName);
368 <    
144 <    //extract the molecule stamps
145 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
146 <    compList(stamps, simParams, moleculeStampPairs);
147 <    
368 >    ff->setFortranForceOptions();
369      //create SimInfo
370 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
371 <    
372 <    //gather parameters (SimCreator only retrieves part of the parameters)
370 >    SimInfo * info = new SimInfo(ff, simParams);
371 >
372 >    info->setRawMetaData(mdRawData);
373 >    
374 >    //gather parameters (SimCreator only retrieves part of the
375 >    //parameters)
376      gatherParameters(info, mdFileName);
377      
378      //divide the molecules and determine the global index of molecules
# Line 160 | Line 384 | namespace oopse {
384      createMolecules(info);
385      
386      
387 <    //allocate memory for DataStorage(circular reference, need to break it)
387 >    //allocate memory for DataStorage(circular reference, need to
388 >    //break it)
389      info->setSnapshotManager(new SimSnapshotManager(info));
390      
391 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
392 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
393 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
391 >    //set the global index of atoms, rigidbodies and cutoffgroups
392 >    //(only need to be set once, the global index will never change
393 >    //again). Local indices of atoms and rigidbodies are already set
394 >    //by MoleculeCreator class which actually delegates the
395 >    //responsibility to LocalIndexManager.
396      setGlobalIndex(info);
397      
398 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
399 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
400 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
401 <    //we can determine the beginning global indices of atoms before they get created.
398 >    //Although addExcludePairs is called inside SimInfo's addMolecule
399 >    //method, at that point atoms don't have the global index yet
400 >    //(their global index are all initialized to -1).  Therefore we
401 >    //have to call addExcludePairs explicitly here. A way to work
402 >    //around is that we can determine the beginning global indices of
403 >    //atoms before they get created.
404      SimInfo::MoleculeIterator mi;
405      Molecule* mol;
406      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
407        info->addExcludePairs(mol);
408      }
409      
181 –    
182 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
183 –    //eta, chi for NPT integrator)
410      if (loadInitCoords)
411 <      loadCoordinates(info);    
411 >      loadCoordinates(info, mdFileName);    
412      
413      return info;
414    }
415    
416    void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
417      
418 <    //figure out the ouput file names
418 >    //figure out the output file names
419      std::string prefix;
420      
421   #ifdef IS_MPI
# Line 218 | Line 444 | namespace oopse {
444    
445   #ifdef IS_MPI
446    void SimCreator::divideMolecules(SimInfo *info) {
447 <    double numerator;
448 <    double denominator;
449 <    double precast;
450 <    double x;
451 <    double y;
452 <    double a;
447 >    RealType numerator;
448 >    RealType denominator;
449 >    RealType precast;
450 >    RealType x;
451 >    RealType y;
452 >    RealType a;
453      int old_atoms;
454      int add_atoms;
455      int new_atoms;
# Line 335 | Line 561 | namespace oopse {
561            //           Pacc(x) = exp(- a * x)
562            // where a = penalty / (average atoms per molecule)
563            
564 <          x = (double)(new_atoms - nTarget);
564 >          x = (RealType)(new_atoms - nTarget);
565            y = myRandom->rand();
566            
567            if (y < exp(- a * x)) {
# Line 395 | Line 621 | namespace oopse {
621        
622      } //end for(int i=0)  
623    }
398 –  
399 –  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
400 –                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
401 –    int i;
402 –    char * id;
403 –    LinkedMolStamp* extractedStamp = NULL;
404 –    MoleculeStamp * currentStamp;
405 –    Component** the_components = simParams->getComponents();
406 –    int n_components = simParams->getNComponents();
624      
408 –    if (!simParams->haveNMol()) {
409 –      // we don't have the total number of molecules, so we assume it is
410 –      // given in each component
411 –      
412 –      for(i = 0; i < n_components; i++) {
413 –        if (!the_components[i]->haveNMol()) {
414 –          // we have a problem
415 –          sprintf(painCave.errMsg,
416 –                  "SimCreator Error. No global NMol or component NMol given.\n"
417 –                  "\tCannot calculate the number of atoms.\n");
418 –          
419 –          painCave.isFatal = 1;
420 –          simError();
421 –        }
422 –        
423 –        id = the_components[i]->getType();
424 –        
425 –        extractedStamp = stamps->extractMolStamp(id);
426 –        if (extractedStamp == NULL) {
427 –          sprintf(painCave.errMsg,
428 –                  "SimCreator error: Component \"%s\" was not found in the "
429 –                  "list of declared molecules\n", id);
430 –          
431 –          painCave.isFatal = 1;
432 –          simError();
433 –        }
434 –        
435 –        currentStamp = extractedStamp->getStamp();
436 –        
437 –        
438 –        moleculeStampPairs.push_back(
439 –                                     std::make_pair(currentStamp, the_components[i]->getNMol()));
440 –      } //end for (i = 0; i < n_components; i++)
441 –    } else {
442 –      sprintf(painCave.errMsg, "SimSetup error.\n"
443 –              "\tSorry, the ability to specify total"
444 –              " nMols and then give molfractions in the components\n"
445 –              "\tis not currently supported."
446 –              " Please give nMol in the components.\n");
447 –      
448 –      painCave.isFatal = 1;
449 –      simError();
450 –    }
451 –    
452 – #ifdef IS_MPI
453 –    
454 –    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
455 –    MPIcheckPoint();
456 –    
457 – #endif // is_mpi
458 –    
459 –  }
460 –  
625    void SimCreator::setGlobalIndex(SimInfo *info) {
626      SimInfo::MoleculeIterator mi;
627      Molecule::AtomIterator ai;
628      Molecule::RigidBodyIterator ri;
629      Molecule::CutoffGroupIterator ci;
630 +    Molecule::IntegrableObjectIterator  ioi;
631      Molecule * mol;
632      Atom * atom;
633      RigidBody * rb;
# Line 471 | Line 636 | namespace oopse {
636      int beginRigidBodyIndex;
637      int beginCutoffGroupIndex;
638      int nGlobalAtoms = info->getNGlobalAtoms();
639 <    
639 >
640 >    /**@todo fixme */
641   #ifndef IS_MPI
642      
643      beginAtomIndex = 0;
# Line 586 | Line 752 | namespace oopse {
752   #else
753      info->setGlobalMolMembership(globalMolMembership);
754   #endif
755 +
756 +    // nIOPerMol holds the number of integrable objects per molecule
757 +    // here the molecules are listed by their global indices.
758 +
759 +    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
760 +    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
761 +      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
762 +    }
763      
764 + #ifdef IS_MPI
765 +    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
766 +    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
767 +                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
768 + #else
769 +    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
770 + #endif    
771 +
772 + std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
773 +
774 + int startingIndex = 0;
775 + for (int i = 0; i < info->getNGlobalMolecules(); i++) {
776 +  startingIOIndexForMol[i] = startingIndex;
777 +  startingIndex += numIntegrableObjectsPerMol[i];
778 + }
779 +
780 + std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
781 + for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
782 +      int myGlobalIndex = mol->getGlobalIndex();
783 +      int globalIO = startingIOIndexForMol[myGlobalIndex];
784 +      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
785 +           integrableObject = mol->nextIntegrableObject(ioi)) {
786 +            integrableObject->setGlobalIntegrableObjectIndex(globalIO);
787 +            IOIndexToIntegrableObject[globalIO] = integrableObject;
788 +            globalIO++;
789 +      }
790 +    }
791 +
792 +  info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
793 +  
794    }
795    
796 <  void SimCreator::loadCoordinates(SimInfo* info) {
796 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
797      Globals* simParams;
798      simParams = info->getSimParams();
799      
596 –    if (!simParams->haveInitialConfig()) {
597 –      sprintf(painCave.errMsg,
598 –              "Cannot intialize a simulation without an initial configuration file.\n");
599 –      painCave.isFatal = 1;;
600 –      simError();
601 –    }
800      
801 <    DumpReader reader(info, simParams->getInitialConfig());
801 >    DumpReader reader(info, mdFileName);
802      int nframes = reader.getNFrames();
803      
804      if (nframes > 0) {
# Line 609 | Line 807 | namespace oopse {
807        //invalid initial coordinate file
808        sprintf(painCave.errMsg,
809                "Initial configuration file %s should at least contain one frame\n",
810 <              simParams->getInitialConfig());
810 >              mdFileName.c_str());
811        painCave.isFatal = 1;
812        simError();
813      }

Diff Legend

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