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 285 by tim, Fri Feb 4 05:26:30 2005 UTC vs.
Revision 1241 by gezelter, Fri Apr 25 15:14:47 2008 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# 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  
50 – #include <sprng.h>
51 –
54   #include "brains/MoleculeCreator.hpp"
55   #include "brains/SimCreator.hpp"
56   #include "brains/SimSnapshotManager.hpp"
57   #include "io/DumpReader.hpp"
56 – #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
79 < #include "io/mpiBASS.h"
62 < #include "math/randomSPRNG.hpp"
79 > #include "math/ParallelRandNumGen.hpp"
80   #endif
81  
82   namespace oopse {
83 +  
84 +  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
85 +    Globals* simParams = NULL;
86 +    try {
87  
88 < void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){
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 < #ifdef IS_MPI
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 <    if (worldRank == 0) {
122 < #endif // is_mpi
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 >      // Create a parser that reads from the scanner
129 >      MDParser parser(lexer);
130 >      parser.setFilename(filename);
131  
132 <        simParams->initalize();
133 <        set_interface_stamps(stamps, simParams);
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 >      antlr::ASTFactory factory;
139 >      parser.initializeASTFactory(factory);
140 >      parser.setASTFactory(&factory);
141 >      parser.mdfile();
142  
143 < #ifdef IS_MPI
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 <        mpiEventInit();
149 >    }
150  
151 < #endif
151 >      
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 >    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 >    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 <        yacc_BASS(mdFileName.c_str());
233 >    return simParams;
234 >  }
235 >  
236 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
237 >                                  bool loadInitCoords) {
238  
239 < #ifdef IS_MPI
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 <        throwMPIEvent(NULL);
249 <    } else {
250 <        set_interface_stamps(stamps, simParams);
251 <        mpiEventInit();
91 <        MPIcheckPoint();
92 <        mpiEventLoop();
93 <    }
248 > #ifdef IS_MPI            
249 >    const int masterNode = 0;
250 >    if (worldRank == masterNode) {
251 > #endif
252  
253 < #endif
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 < }
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 < SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
276 <    
277 <    MakeStamps * stamps = new MakeStamps();
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 <    Globals * simParams = new Globals();
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 <    //parse meta-data file
314 <    parseFile(mdFileName, stamps, simParams);
315 <
316 <    //create the force field
317 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
318 <                          simParams->getForceField());
319 <    
320 <    if (ff == NULL) {
321 <        sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
322 <                simParams->getForceField());
323 <        painCave.isFatal = 1;
116 <        simError();
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 +    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
331 +    
332 +    //create the force field
333 +    ForceField * ff = ForceFieldFactory::getInstance()
334 +      ->createForceField(simParams->getForceField());
335 +    
336 +    if (ff == NULL) {
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 +    }
343 +    
344      if (simParams->haveForceFieldFileName()) {
345 <        ff->setForceFieldFileName(simParams->getForceFieldFileName());
345 >      ff->setForceFieldFileName(simParams->getForceFieldFileName());
346      }
347      
348      std::string forcefieldFileName;
349      forcefieldFileName = ff->getForceFieldFileName();
350 <
350 >    
351      if (simParams->haveForceFieldVariant()) {
352 <        //If the force field has variant, the variant force field name will be
353 <        //Base.variant.frc. For exampel EAM.u6.frc
354 <        
355 <        std::string variant = simParams->getForceFieldVariant();
356 <
357 <        std::string::size_type pos = forcefieldFileName.rfind(".frc");
358 <        variant = "." + variant;
359 <        if (pos != std::string::npos) {
360 <            forcefieldFileName.insert(pos, variant);
361 <        } else {
362 <            //If the default force field file name does not containt .frc suffix, just append the .variant
363 <            forcefieldFileName.append(variant);
364 <        }
352 >      //If the force field has variant, the variant force field name will be
353 >      //Base.variant.frc. For exampel EAM.u6.frc
354 >      
355 >      std::string variant = simParams->getForceFieldVariant();
356 >      
357 >      std::string::size_type pos = forcefieldFileName.rfind(".frc");
358 >      variant = "." + variant;
359 >      if (pos != std::string::npos) {
360 >        forcefieldFileName.insert(pos, variant);
361 >      } else {
362 >        //If the default force field file name does not containt .frc suffix, just append the .variant
363 >        forcefieldFileName.append(variant);
364 >      }
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);
370 >    SimInfo * info = new SimInfo(ff, simParams);
371  
372 <    //gather parameters (SimCreator only retrieves part of the parameters)
372 >    info->setRawMetaData(mdRawData);
373 >    
374 >    //gather parameters (SimCreator only retrieves part of the
375 >    //parameters)
376      gatherParameters(info, mdFileName);
377 <
377 >    
378      //divide the molecules and determine the global index of molecules
379   #ifdef IS_MPI
380      divideMolecules(info);
381   #endif
382 <
382 >    
383      //create the molecules
384      createMolecules(info);
385 <
386 <
387 <    //allocate memory for DataStorage(circular reference, need to break it)
385 >    
386 >    
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.
397 >    
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);
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 <    //setup seed for random number generator
193 <    int seedValue;
194 <    Globals * simParams = info->getSimParams();
195 <
196 <    if (simParams->haveSeed()) {
197 <        seedValue = simParams->getSeed();
198 <
199 <        if (seedValue < 100000000 ) {
200 <            sprintf(painCave.errMsg,
201 <                    "Seed for sprng library should contain at least 9 digits\n"
202 <                        "OOPSE will generate a seed for user\n");
203 <
204 <            painCave.isFatal = 0;
205 <            simError();
206 <
207 <            //using seed generated by system instead of invalid seed set by user
208 <
209 < #ifndef IS_MPI
210 <
211 <            seedValue = make_sprng_seed();
212 <
213 < #else
214 <
215 <            if (worldRank == 0) {
216 <                seedValue = make_sprng_seed();
217 <            }
218 <
219 <            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
220 <
221 < #endif
222 <
223 <        } //end if (seedValue /1000000000 == 0)
224 <    } else {
225 <
226 < #ifndef IS_MPI
227 <
228 <        seedValue = make_sprng_seed();
229 <
230 < #else
231 <
232 <        if (worldRank == 0) {
233 <            seedValue = make_sprng_seed();
234 <        }
235 <
236 <        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
237 <
238 < #endif
239 <
240 <    } //end of simParams->haveSeed()
241 <
242 <    info->setSeed(seedValue);
243 <
244 <
245 <    //figure out the ouput file names
414 >  }
415 >  
416 >  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
417 >    
418 >    //figure out the output file names
419      std::string prefix;
420 <
420 >    
421   #ifdef IS_MPI
422 <
422 >    
423      if (worldRank == 0) {
424   #endif // is_mpi
425 <
426 <        if (simParams->haveFinalConfig()) {
427 <            prefix = getPrefix(simParams->getFinalConfig());
428 <        } else {
429 <            prefix = getPrefix(mdfile);
430 <        }
431 <
432 <        info->setFinalConfigFileName(prefix + ".eor");
433 <        info->setDumpFileName(prefix + ".dump");
434 <        info->setStatFileName(prefix + ".stat");
435 <
425 >      Globals * simParams = info->getSimParams();
426 >      if (simParams->haveFinalConfig()) {
427 >        prefix = getPrefix(simParams->getFinalConfig());
428 >      } else {
429 >        prefix = getPrefix(mdfile);
430 >      }
431 >      
432 >      info->setFinalConfigFileName(prefix + ".eor");
433 >      info->setDumpFileName(prefix + ".dump");
434 >      info->setStatFileName(prefix + ".stat");
435 >      info->setRestFileName(prefix + ".zang");
436 >      
437   #ifdef IS_MPI
438 <
438 >      
439      }
440 <
440 >    
441   #endif
442 <
443 < }
444 <
442 >    
443 >  }
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;
446 >  void SimCreator::divideMolecules(SimInfo *info) {
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 287 | Line 461 | void SimCreator::divideMolecules(SimInfo *info) {
461      int which_proc;
462      int nProcessors;
463      std::vector<int> atomsPerProc;
290 –    randomSPRNG myRandom(info->getSeed());
464      int nGlobalMols = info->getNGlobalMolecules();
465      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
466      
467      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
468 <
468 >    
469      if (nProcessors > nGlobalMols) {
470 <        sprintf(painCave.errMsg,
471 <                "nProcessors (%d) > nMol (%d)\n"
472 <                    "\tThe number of processors is larger than\n"
473 <                    "\tthe number of molecules.  This will not result in a \n"
474 <                    "\tusable division of atoms for force decomposition.\n"
475 <                    "\tEither try a smaller number of processors, or run the\n"
476 <                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
477 <
478 <        painCave.isFatal = 1;
479 <        simError();
480 <    }
481 <
470 >      sprintf(painCave.errMsg,
471 >              "nProcessors (%d) > nMol (%d)\n"
472 >              "\tThe number of processors is larger than\n"
473 >              "\tthe number of molecules.  This will not result in a \n"
474 >              "\tusable division of atoms for force decomposition.\n"
475 >              "\tEither try a smaller number of processors, or run the\n"
476 >              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
477 >      
478 >      painCave.isFatal = 1;
479 >      simError();
480 >    }
481 >    
482 >    int seedValue;
483 >    Globals * simParams = info->getSimParams();
484 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
485 >    if (simParams->haveSeed()) {
486 >      seedValue = simParams->getSeed();
487 >      myRandom = new SeqRandNumGen(seedValue);
488 >    }else {
489 >      myRandom = new SeqRandNumGen();
490 >    }  
491 >    
492 >    
493      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
494 <
494 >    
495      //initialize atomsPerProc
496      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
497 <
497 >    
498      if (worldRank == 0) {
499 <        numerator = info->getNGlobalAtoms();
500 <        denominator = nProcessors;
501 <        precast = numerator / denominator;
502 <        nTarget = (int)(precast + 0.5);
503 <
504 <        for(i = 0; i < nGlobalMols; i++) {
505 <            done = 0;
506 <            loops = 0;
507 <
508 <            while (!done) {
509 <                loops++;
510 <
511 <                // Pick a processor at random
512 <
513 <                which_proc = (int) (myRandom.getRandom() * nProcessors);
514 <
515 <                //get the molecule stamp first
516 <                int stampId = info->getMoleculeStampId(i);
517 <                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
518 <
519 <                // How many atoms does this processor have so far?
520 <                old_atoms = atomsPerProc[which_proc];
521 <                add_atoms = moleculeStamp->getNAtoms();
522 <                new_atoms = old_atoms + add_atoms;
523 <
524 <                // If we've been through this loop too many times, we need
525 <                // to just give up and assign the molecule to this processor
526 <                // and be done with it.
527 <
528 <                if (loops > 100) {
529 <                    sprintf(painCave.errMsg,
530 <                            "I've tried 100 times to assign molecule %d to a "
531 <                                " processor, but can't find a good spot.\n"
532 <                                "I'm assigning it at random to processor %d.\n",
533 <                            i, which_proc);
534 <
535 <                    painCave.isFatal = 0;
536 <                    simError();
537 <
538 <                    molToProcMap[i] = which_proc;
539 <                    atomsPerProc[which_proc] += add_atoms;
540 <
541 <                    done = 1;
542 <                    continue;
543 <                }
544 <
545 <                // If we can add this molecule to this processor without sending
546 <                // it above nTarget, then go ahead and do it:
547 <
548 <                if (new_atoms <= nTarget) {
549 <                    molToProcMap[i] = which_proc;
550 <                    atomsPerProc[which_proc] += add_atoms;
551 <
552 <                    done = 1;
553 <                    continue;
554 <                }
555 <
556 <                // The only situation left is when new_atoms > nTarget.  We
557 <                // want to accept this with some probability that dies off the
558 <                // farther we are from nTarget
559 <
560 <                // roughly:  x = new_atoms - nTarget
561 <                //           Pacc(x) = exp(- a * x)
562 <                // where a = penalty / (average atoms per molecule)
563 <
564 <                x = (double)(new_atoms - nTarget);
565 <                y = myRandom.getRandom();
566 <
567 <                if (y < exp(- a * x)) {
568 <                    molToProcMap[i] = which_proc;
569 <                    atomsPerProc[which_proc] += add_atoms;
570 <
571 <                    done = 1;
572 <                    continue;
573 <                } else {
574 <                    continue;
575 <                }
392 <            }
499 >      numerator = info->getNGlobalAtoms();
500 >      denominator = nProcessors;
501 >      precast = numerator / denominator;
502 >      nTarget = (int)(precast + 0.5);
503 >      
504 >      for(i = 0; i < nGlobalMols; i++) {
505 >        done = 0;
506 >        loops = 0;
507 >        
508 >        while (!done) {
509 >          loops++;
510 >          
511 >          // Pick a processor at random
512 >          
513 >          which_proc = (int) (myRandom->rand() * nProcessors);
514 >          
515 >          //get the molecule stamp first
516 >          int stampId = info->getMoleculeStampId(i);
517 >          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
518 >          
519 >          // How many atoms does this processor have so far?
520 >          old_atoms = atomsPerProc[which_proc];
521 >          add_atoms = moleculeStamp->getNAtoms();
522 >          new_atoms = old_atoms + add_atoms;
523 >          
524 >          // If we've been through this loop too many times, we need
525 >          // to just give up and assign the molecule to this processor
526 >          // and be done with it.
527 >          
528 >          if (loops > 100) {
529 >            sprintf(painCave.errMsg,
530 >                    "I've tried 100 times to assign molecule %d to a "
531 >                    " processor, but can't find a good spot.\n"
532 >                    "I'm assigning it at random to processor %d.\n",
533 >                    i, which_proc);
534 >            
535 >            painCave.isFatal = 0;
536 >            simError();
537 >            
538 >            molToProcMap[i] = which_proc;
539 >            atomsPerProc[which_proc] += add_atoms;
540 >            
541 >            done = 1;
542 >            continue;
543 >          }
544 >          
545 >          // If we can add this molecule to this processor without sending
546 >          // it above nTarget, then go ahead and do it:
547 >          
548 >          if (new_atoms <= nTarget) {
549 >            molToProcMap[i] = which_proc;
550 >            atomsPerProc[which_proc] += add_atoms;
551 >            
552 >            done = 1;
553 >            continue;
554 >          }
555 >          
556 >          // The only situation left is when new_atoms > nTarget.  We
557 >          // want to accept this with some probability that dies off the
558 >          // farther we are from nTarget
559 >          
560 >          // roughly:  x = new_atoms - nTarget
561 >          //           Pacc(x) = exp(- a * x)
562 >          // where a = penalty / (average atoms per molecule)
563 >          
564 >          x = (RealType)(new_atoms - nTarget);
565 >          y = myRandom->rand();
566 >          
567 >          if (y < exp(- a * x)) {
568 >            molToProcMap[i] = which_proc;
569 >            atomsPerProc[which_proc] += add_atoms;
570 >            
571 >            done = 1;
572 >            continue;
573 >          } else {
574 >            continue;
575 >          }
576          }
577 <
578 <        // Spray out this nonsense to all other processors:
579 <
580 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
577 >      }
578 >      
579 >      delete myRandom;
580 >      
581 >      // Spray out this nonsense to all other processors:
582 >      
583 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
584      } else {
585 <
586 <        // Listen to your marching orders from processor 0:
587 <
588 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
585 >      
586 >      // Listen to your marching orders from processor 0:
587 >      
588 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
589      }
590 <
590 >    
591      info->setMolToProcMap(molToProcMap);
592      sprintf(checkPointMsg,
593              "Successfully divided the molecules among the processors.\n");
594 <    MPIcheckPoint();
595 < }
596 <
594 >    errorCheckPoint();
595 >  }
596 >  
597   #endif
598 <
599 < void SimCreator::createMolecules(SimInfo *info) {
598 >  
599 >  void SimCreator::createMolecules(SimInfo *info) {
600      MoleculeCreator molCreator;
601      int stampId;
602 <
602 >    
603      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
604 <
604 >      
605   #ifdef IS_MPI
606 <
607 <        if (info->getMolToProc(i) == worldRank) {
606 >      
607 >      if (info->getMolToProc(i) == worldRank) {
608   #endif
609 <
610 <            stampId = info->getMoleculeStampId(i);
611 <            Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
612 <                                                                                    stampId, i, info->getLocalIndexManager());
613 <
614 <            info->addMolecule(mol);
615 <
609 >        
610 >        stampId = info->getMoleculeStampId(i);
611 >        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
612 >                                                   stampId, i, info->getLocalIndexManager());
613 >        
614 >        info->addMolecule(mol);
615 >        
616   #ifdef IS_MPI
617 <
618 <        }
619 <
617 >        
618 >      }
619 >      
620   #endif
621 <
621 >      
622      } //end for(int i=0)  
623 < }
624 <
625 < void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
440 <                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
441 <    int i;
442 <    char * id;
443 <    MoleculeStamp * currentStamp;
444 <    Component** the_components = simParams->getComponents();
445 <    int n_components = simParams->getNComponents();
446 <
447 <    if (!simParams->haveNMol()) {
448 <        // we don't have the total number of molecules, so we assume it is
449 <        // given in each component
450 <
451 <        for(i = 0; i < n_components; i++) {
452 <            if (!the_components[i]->haveNMol()) {
453 <                // we have a problem
454 <                sprintf(painCave.errMsg,
455 <                        "SimCreator Error. No global NMol or component NMol given.\n"
456 <                            "\tCannot calculate the number of atoms.\n");
457 <
458 <                painCave.isFatal = 1;
459 <                simError();
460 <            }
461 <
462 <            id = the_components[i]->getType();
463 <            currentStamp = (stamps->extractMolStamp(id))->getStamp();
464 <
465 <            if (currentStamp == NULL) {
466 <                sprintf(painCave.errMsg,
467 <                        "SimCreator error: Component \"%s\" was not found in the "
468 <                            "list of declared molecules\n", id);
469 <
470 <                painCave.isFatal = 1;
471 <                simError();
472 <            }
473 <
474 <            moleculeStampPairs.push_back(
475 <                std::make_pair(currentStamp, the_components[i]->getNMol()));
476 <        } //end for (i = 0; i < n_components; i++)
477 <    } else {
478 <        sprintf(painCave.errMsg, "SimSetup error.\n"
479 <                                     "\tSorry, the ability to specify total"
480 <                                     " nMols and then give molfractions in the components\n"
481 <                                     "\tis not currently supported."
482 <                                     " Please give nMol in the components.\n");
483 <
484 <        painCave.isFatal = 1;
485 <        simError();
486 <    }
487 <
488 < #ifdef IS_MPI
489 <
490 <    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
491 <    MPIcheckPoint();
492 <
493 < #endif // is_mpi
494 <
495 < }
496 <
497 < void SimCreator::setGlobalIndex(SimInfo *info) {
623 >  }
624 >    
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 507 | Line 636 | void SimCreator::setGlobalIndex(SimInfo *info) {
636      int beginRigidBodyIndex;
637      int beginCutoffGroupIndex;
638      int nGlobalAtoms = info->getNGlobalAtoms();
510 –    
511 – #ifndef IS_MPI
639  
640 +    /**@todo fixme */
641 + #ifndef IS_MPI
642 +    
643      beginAtomIndex = 0;
644      beginRigidBodyIndex = 0;
645      beginCutoffGroupIndex = 0;
646 <
646 >    
647   #else
648 <
648 >    
649      int nproc;
650      int myNode;
651 <
651 >    
652      myNode = worldRank;
653      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
654 <
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 <
661 >    
662      tmpAtomsInProc[myNode] = info->getNAtoms();
663      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
664      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
665 <
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);
# Line 540 | Line 670 | void SimCreator::setGlobalIndex(SimInfo *info) {
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 <
673 >    
674      beginAtomIndex = 0;
675      beginRigidBodyIndex = 0;
676      beginCutoffGroupIndex = 0;
677 <
677 >    
678      for(int i = 0; i < myNode; i++) {
679 <        beginAtomIndex += NumAtomsInProc[i];
680 <        beginRigidBodyIndex += NumRigidBodiesInProc[i];
681 <        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
679 >      beginAtomIndex += NumAtomsInProc[i];
680 >      beginRigidBodyIndex += NumRigidBodiesInProc[i];
681 >      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
682      }
683 <
683 >    
684 > #endif
685 >    
686      //rigidbody's index begins right after atom's
687      beginRigidBodyIndex += info->getNGlobalAtoms();
688 < #endif
557 <
688 >    
689      for(mol = info->beginMolecule(mi); mol != NULL;
690          mol = info->nextMolecule(mi)) {
691 <
692 <        //local index(index in DataStorge) of atom is important
693 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
694 <            atom->setGlobalIndex(beginAtomIndex++);
695 <        }
696 <
697 <        for(rb = mol->beginRigidBody(ri); rb != NULL;
698 <            rb = mol->nextRigidBody(ri)) {
699 <            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 <        }
691 >      
692 >      //local index(index in DataStorge) of atom is important
693 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
694 >        atom->setGlobalIndex(beginAtomIndex++);
695 >      }
696 >      
697 >      for(rb = mol->beginRigidBody(ri); rb != NULL;
698 >          rb = mol->nextRigidBody(ri)) {
699 >        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 <
708 >    
709      //fill globalGroupMembership
710      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
711      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
712 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
713 <
714 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
715 <                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
716 <            }
717 <
718 <        }      
712 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
713 >        
714 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
715 >          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
716 >        }
717 >        
718 >      }      
719      }
720 <
720 >    
721   #ifdef IS_MPI    
722      // Since the globalGroupMembership has been zero filled and we've only
723      // poked values into the atoms we know, we can do an Allreduce
# Line 596 | Line 727 | void SimCreator::setGlobalIndex(SimInfo *info) {
727      std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
728      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
729                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
730 <     info->setGlobalGroupMembership(tmpGroupMembership);
730 >    info->setGlobalGroupMembership(tmpGroupMembership);
731   #else
732      info->setGlobalGroupMembership(globalGroupMembership);
733   #endif
734 <
734 >    
735      //fill molMembership
736      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
737      
738      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
739 <
740 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
741 <            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
742 <        }
739 >      
740 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
741 >        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
742 >      }
743      }
744 <
744 >    
745   #ifdef IS_MPI
746      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
747 <
747 >    
748      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
749                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
750      
# Line 622 | Line 753 | void SimCreator::setGlobalIndex(SimInfo *info) {
753      info->setGlobalMolMembership(globalMolMembership);
754   #endif
755  
756 < }
756 >    // nIOPerMol holds the number of integrable objects per molecule
757 >    // here the molecules are listed by their global indices.
758  
759 < void SimCreator::loadCoordinates(SimInfo* info) {
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, const std::string& mdFileName) {
797      Globals* simParams;
798      simParams = info->getSimParams();
799      
800 <    if (!simParams->haveInitialConfig()) {
801 <        sprintf(painCave.errMsg,
633 <                "Cannot intialize a simulation without an initial configuration file.\n");
634 <        painCave.isFatal = 1;;
635 <        simError();
636 <    }
637 <        
638 <    DumpReader reader(info, simParams->getInitialConfig());
800 >    
801 >    DumpReader reader(info, mdFileName);
802      int nframes = reader.getNFrames();
803 <
803 >    
804      if (nframes > 0) {
805 <        reader.readFrame(nframes - 1);
805 >      reader.readFrame(nframes - 1);
806      } else {
807 <        //invalid initial coordinate file
808 <        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
809 <                simParams->getInitialConfig());
810 <        painCave.isFatal = 1;
811 <        simError();
807 >      //invalid initial coordinate file
808 >      sprintf(painCave.errMsg,
809 >              "Initial configuration file %s should at least contain one frame\n",
810 >              mdFileName.c_str());
811 >      painCave.isFatal = 1;
812 >      simError();
813      }
814 <
814 >    
815      //copy the current snapshot to previous snapshot
816      info->getSnapshotManager()->advance();
817 < }
818 <
817 >  }
818 >  
819   } //end namespace oopse
820  
821  

Diff Legend

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