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 403 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
Revision 1287 by gezelter, Wed Sep 10 18:11:32 2008 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 +  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  
150 <      mpiEventInit();
150 >      
151 >    catch(antlr::MismatchedCharException& e) {
152 >      sprintf(painCave.errMsg,
153 >              "parser exception: %s %s:%d:%d\n",
154 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
155 >      painCave.isFatal = 1;
156 >      simError();          
157 >    }
158 >    catch(antlr::MismatchedTokenException &e) {
159 >      sprintf(painCave.errMsg,
160 >              "parser exception: %s %s:%d:%d\n",
161 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
162 >      painCave.isFatal = 1;
163 >      simError();  
164 >    }
165 >    catch(antlr::NoViableAltForCharException &e) {
166 >      sprintf(painCave.errMsg,
167 >              "parser exception: %s %s:%d:%d\n",
168 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
169 >      painCave.isFatal = 1;
170 >      simError();  
171 >    }
172 >    catch(antlr::NoViableAltException &e) {
173 >      sprintf(painCave.errMsg,
174 >              "parser exception: %s %s:%d:%d\n",
175 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
176 >      painCave.isFatal = 1;
177 >      simError();  
178 >    }
179 >      
180 >    catch(antlr::TokenStreamRecognitionException& e) {
181 >      sprintf(painCave.errMsg,
182 >              "parser exception: %s %s:%d:%d\n",
183 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
184 >      painCave.isFatal = 1;
185 >      simError();  
186 >    }
187 >        
188 >    catch(antlr::TokenStreamIOException& e) {
189 >      sprintf(painCave.errMsg,
190 >              "parser exception: %s\n",
191 >              e.getMessage().c_str());
192 >      painCave.isFatal = 1;
193 >      simError();
194 >    }
195 >        
196 >    catch(antlr::TokenStreamException& e) {
197 >      sprintf(painCave.errMsg,
198 >              "parser exception: %s\n",
199 >              e.getMessage().c_str());
200 >      painCave.isFatal = 1;
201 >      simError();
202 >    }        
203 >    catch (antlr::RecognitionException& e) {
204 >      sprintf(painCave.errMsg,
205 >              "parser exception: %s %s:%d:%d\n",
206 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
207 >      painCave.isFatal = 1;
208 >      simError();          
209 >    }
210 >    catch (antlr::CharStreamException& e) {
211 >      sprintf(painCave.errMsg,
212 >              "parser exception: %s\n",
213 >              e.getMessage().c_str());
214 >      painCave.isFatal = 1;
215 >      simError();        
216 >    }
217 >    catch (OOPSEException& e) {
218 >      sprintf(painCave.errMsg,
219 >              "%s\n",
220 >              e.getMessage().c_str());
221 >      painCave.isFatal = 1;
222 >      simError();
223 >    }
224 >    catch (std::exception& e) {
225 >      sprintf(painCave.errMsg,
226 >              "parser exception: %s\n",
227 >              e.what());
228 >      painCave.isFatal = 1;
229 >      simError();
230 >    }
231  
232 < #endif
232 >    return simParams;
233 >  }
234 >  
235 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
236 >                                  bool loadInitCoords) {
237  
238 <      yacc_BASS(mdFileName.c_str());
238 >    const int bufferSize = 65535;
239 >    char buffer[bufferSize];
240 >    int lineNo = 0;
241 >    std::string mdRawData;
242 >    int metaDataBlockStart = -1;
243 >    int metaDataBlockEnd = -1;
244 >    int i;
245 >    int mdOffset;
246  
247 < #ifdef IS_MPI
247 > #ifdef IS_MPI            
248 >    const int masterNode = 0;
249 >    if (worldRank == masterNode) {
250 > #endif
251  
252 <      throwMPIEvent(NULL);
253 <    } else {
254 <      set_interface_stamps(stamps, simParams);
255 <      mpiEventInit();
256 <      MPIcheckPoint();
257 <      mpiEventLoop();
258 <    }
252 >      std::ifstream mdFile_(mdFileName.c_str());
253 >      
254 >      if (mdFile_.fail()) {
255 >        sprintf(painCave.errMsg,
256 >                "SimCreator: Cannot open file: %s\n",
257 >                mdFileName.c_str());
258 >        painCave.isFatal = 1;
259 >        simError();
260 >      }
261  
262 < #endif
262 >      mdFile_.getline(buffer, bufferSize);
263 >      ++lineNo;
264 >      std::string line = trimLeftCopy(buffer);
265 >      i = CaseInsensitiveFind(line, "<OOPSE");
266 >      if (static_cast<size_t>(i) == string::npos) {
267 >        sprintf(painCave.errMsg,
268 >                "SimCreator: File: %s is not an OOPSE file!\n",
269 >                mdFileName.c_str());
270 >        painCave.isFatal = 1;
271 >        simError();
272 >      }
273  
274 <  }
274 >      //scan through the input stream and find MetaData tag        
275 >      while(mdFile_.getline(buffer, bufferSize)) {
276 >        ++lineNo;
277 >        
278 >        std::string line = trimLeftCopy(buffer);
279 >        if (metaDataBlockStart == -1) {
280 >          i = CaseInsensitiveFind(line, "<MetaData>");
281 >          if (i != string::npos) {
282 >            metaDataBlockStart = lineNo;
283 >            mdOffset = mdFile_.tellg();
284 >          }
285 >        } else {
286 >          i = CaseInsensitiveFind(line, "</MetaData>");
287 >          if (i != string::npos) {
288 >            metaDataBlockEnd = lineNo;
289 >          }
290 >        }
291 >      }
292  
293 <  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
294 <    
295 <    MakeStamps * stamps = new MakeStamps();
293 >      if (metaDataBlockStart == -1) {
294 >        sprintf(painCave.errMsg,
295 >                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
296 >                mdFileName.c_str());
297 >        painCave.isFatal = 1;
298 >        simError();
299 >      }
300 >      if (metaDataBlockEnd == -1) {
301 >        sprintf(painCave.errMsg,
302 >                "SimCreator: File: %s did not contain a closed MetaData block!\n",
303 >                mdFileName.c_str());
304 >        painCave.isFatal = 1;
305 >        simError();
306 >      }
307 >        
308 >      mdFile_.clear();
309 >      mdFile_.seekg(0);
310 >      mdFile_.seekg(mdOffset);
311  
312 <    Globals * simParams = new Globals();
312 >      mdRawData.clear();
313  
314 <    //parse meta-data file
315 <    parseFile(mdFileName, stamps, simParams);
314 >      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
315 >        mdFile_.getline(buffer, bufferSize);
316 >        mdRawData += buffer;
317 >        mdRawData += "\n";
318 >      }
319  
320 <    //create the force field
321 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
322 <                                                                         simParams->getForceField());
320 >      mdFile_.close();
321 >
322 > #ifdef IS_MPI
323 >    }
324 > #endif
325 >
326 >    std::stringstream rawMetaDataStream(mdRawData);
327 >
328 >    //parse meta-data file
329 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
330      
331 +    //create the force field
332 +    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
333 +
334      if (ff == NULL) {
335 <      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
336 <              simParams->getForceField());
335 >      sprintf(painCave.errMsg,
336 >              "ForceField Factory can not create %s force field\n",
337 >              simParams->getForceField().c_str());
338        painCave.isFatal = 1;
339        simError();
340      }
341 <
341 >    
342      if (simParams->haveForceFieldFileName()) {
343        ff->setForceFieldFileName(simParams->getForceFieldFileName());
344      }
345      
346      std::string forcefieldFileName;
347      forcefieldFileName = ff->getForceFieldFileName();
348 <
348 >    
349      if (simParams->haveForceFieldVariant()) {
350        //If the force field has variant, the variant force field name will be
351        //Base.variant.frc. For exampel EAM.u6.frc
352 <        
352 >      
353        std::string variant = simParams->getForceFieldVariant();
354 <
354 >      
355        std::string::size_type pos = forcefieldFileName.rfind(".frc");
356        variant = "." + variant;
357        if (pos != std::string::npos) {
# Line 139 | Line 363 | namespace oopse {
363      }
364      
365      ff->parse(forcefieldFileName);
366 <    
143 <    //extract the molecule stamps
144 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
145 <    compList(stamps, simParams, moleculeStampPairs);
146 <
366 >    ff->setFortranForceOptions();
367      //create SimInfo
368 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
368 >    SimInfo * info = new SimInfo(ff, simParams);
369  
370 <    //gather parameters (SimCreator only retrieves part of the parameters)
370 >    info->setRawMetaData(mdRawData);
371 >    
372 >    //gather parameters (SimCreator only retrieves part of the
373 >    //parameters)
374      gatherParameters(info, mdFileName);
375 <
375 >    
376      //divide the molecules and determine the global index of molecules
377   #ifdef IS_MPI
378      divideMolecules(info);
379   #endif
380 <
380 >    
381      //create the molecules
382      createMolecules(info);
383 <
384 <
385 <    //allocate memory for DataStorage(circular reference, need to break it)
383 >    
384 >    
385 >    //allocate memory for DataStorage(circular reference, need to
386 >    //break it)
387      info->setSnapshotManager(new SimSnapshotManager(info));
388      
389 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
390 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
391 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
389 >    //set the global index of atoms, rigidbodies and cutoffgroups
390 >    //(only need to be set once, the global index will never change
391 >    //again). Local indices of atoms and rigidbodies are already set
392 >    //by MoleculeCreator class which actually delegates the
393 >    //responsibility to LocalIndexManager.
394      setGlobalIndex(info);
395 <
396 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
397 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
398 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
399 <    //we can determine the beginning global indices of atoms before they get created.
395 >    
396 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
397 >    //method, at that point atoms don't have the global index yet
398 >    //(their global index are all initialized to -1).  Therefore we
399 >    //have to call addInteractionPairs explicitly here. A way to work
400 >    //around is that we can determine the beginning global indices of
401 >    //atoms before they get created.
402      SimInfo::MoleculeIterator mi;
403      Molecule* mol;
404      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
405 <      info->addExcludePairs(mol);
405 >      info->addInteractionPairs(mol);
406      }
407      
180 –
181 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
182 –    //eta, chi for NPT integrator)
408      if (loadInitCoords)
409 <      loadCoordinates(info);    
409 >      loadCoordinates(info, mdFileName);    
410      
411      return info;
412    }
413 <
413 >  
414    void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
415 <
416 <    //figure out the ouput file names
415 >    
416 >    //figure out the output file names
417      std::string prefix;
418 <
418 >    
419   #ifdef IS_MPI
420 <
420 >    
421      if (worldRank == 0) {
422   #endif // is_mpi
423        Globals * simParams = info->getSimParams();
# Line 201 | Line 426 | namespace oopse {
426        } else {
427          prefix = getPrefix(mdfile);
428        }
429 <
429 >      
430        info->setFinalConfigFileName(prefix + ".eor");
431        info->setDumpFileName(prefix + ".dump");
432        info->setStatFileName(prefix + ".stat");
433 <
433 >      info->setRestFileName(prefix + ".zang");
434 >      
435   #ifdef IS_MPI
436 <
436 >      
437      }
438 <
438 >    
439   #endif
440 <
440 >    
441    }
442 <
442 >  
443   #ifdef IS_MPI
444    void SimCreator::divideMolecules(SimInfo *info) {
445 <    double numerator;
446 <    double denominator;
447 <    double precast;
448 <    double x;
449 <    double y;
450 <    double a;
445 >    RealType numerator;
446 >    RealType denominator;
447 >    RealType precast;
448 >    RealType x;
449 >    RealType y;
450 >    RealType a;
451      int old_atoms;
452      int add_atoms;
453      int new_atoms;
# Line 237 | Line 463 | namespace oopse {
463      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
464      
465      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
466 <
466 >    
467      if (nProcessors > nGlobalMols) {
468        sprintf(painCave.errMsg,
469                "nProcessors (%d) > nMol (%d)\n"
# Line 246 | Line 472 | namespace oopse {
472                "\tusable division of atoms for force decomposition.\n"
473                "\tEither try a smaller number of processors, or run the\n"
474                "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
475 <
475 >      
476        painCave.isFatal = 1;
477        simError();
478      }
479 <
479 >    
480      int seedValue;
481      Globals * simParams = info->getSimParams();
482      SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
# Line 260 | Line 486 | namespace oopse {
486      }else {
487        myRandom = new SeqRandNumGen();
488      }  
489 <
490 <
489 >    
490 >    
491      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
492 <
492 >    
493      //initialize atomsPerProc
494      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
495 <
495 >    
496      if (worldRank == 0) {
497        numerator = info->getNGlobalAtoms();
498        denominator = nProcessors;
499        precast = numerator / denominator;
500        nTarget = (int)(precast + 0.5);
501 <
501 >      
502        for(i = 0; i < nGlobalMols; i++) {
503          done = 0;
504          loops = 0;
505 <
505 >        
506          while (!done) {
507            loops++;
508 <
508 >          
509            // Pick a processor at random
510 <
510 >          
511            which_proc = (int) (myRandom->rand() * nProcessors);
512 <
512 >          
513            //get the molecule stamp first
514            int stampId = info->getMoleculeStampId(i);
515            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
516 <
516 >          
517            // How many atoms does this processor have so far?
518            old_atoms = atomsPerProc[which_proc];
519            add_atoms = moleculeStamp->getNAtoms();
520            new_atoms = old_atoms + add_atoms;
521 <
521 >          
522            // If we've been through this loop too many times, we need
523            // to just give up and assign the molecule to this processor
524            // and be done with it.
525 <
525 >          
526            if (loops > 100) {
527              sprintf(painCave.errMsg,
528                      "I've tried 100 times to assign molecule %d to a "
529                      " processor, but can't find a good spot.\n"
530                      "I'm assigning it at random to processor %d.\n",
531                      i, which_proc);
532 <
532 >            
533              painCave.isFatal = 0;
534              simError();
535 <
535 >            
536              molToProcMap[i] = which_proc;
537              atomsPerProc[which_proc] += add_atoms;
538 <
538 >            
539              done = 1;
540              continue;
541            }
542 <
542 >          
543            // If we can add this molecule to this processor without sending
544            // it above nTarget, then go ahead and do it:
545 <
545 >          
546            if (new_atoms <= nTarget) {
547              molToProcMap[i] = which_proc;
548              atomsPerProc[which_proc] += add_atoms;
549 <
549 >            
550              done = 1;
551              continue;
552            }
553 <
553 >          
554            // The only situation left is when new_atoms > nTarget.  We
555            // want to accept this with some probability that dies off the
556            // farther we are from nTarget
557 <
557 >          
558            // roughly:  x = new_atoms - nTarget
559            //           Pacc(x) = exp(- a * x)
560            // where a = penalty / (average atoms per molecule)
561 <
562 <          x = (double)(new_atoms - nTarget);
561 >          
562 >          x = (RealType)(new_atoms - nTarget);
563            y = myRandom->rand();
564 <
564 >          
565            if (y < exp(- a * x)) {
566              molToProcMap[i] = which_proc;
567              atomsPerProc[which_proc] += add_atoms;
568 <
568 >            
569              done = 1;
570              continue;
571            } else {
# Line 347 | Line 573 | namespace oopse {
573            }
574          }
575        }
576 <
576 >      
577        delete myRandom;
578 <        
578 >      
579        // Spray out this nonsense to all other processors:
580 <
580 >      
581        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
582      } else {
583 <
583 >      
584        // Listen to your marching orders from processor 0:
585 <
585 >      
586        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
587      }
588 <
588 >    
589      info->setMolToProcMap(molToProcMap);
590      sprintf(checkPointMsg,
591              "Successfully divided the molecules among the processors.\n");
592 <    MPIcheckPoint();
592 >    errorCheckPoint();
593    }
594 <
594 >  
595   #endif
596 <
596 >  
597    void SimCreator::createMolecules(SimInfo *info) {
598      MoleculeCreator molCreator;
599      int stampId;
600 <
600 >    
601      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
602 <
602 >      
603   #ifdef IS_MPI
604 <
604 >      
605        if (info->getMolToProc(i) == worldRank) {
606   #endif
607 <
607 >        
608          stampId = info->getMoleculeStampId(i);
609          Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
610                                                     stampId, i, info->getLocalIndexManager());
611 <
611 >        
612          info->addMolecule(mol);
613 <
613 >        
614   #ifdef IS_MPI
615 <
615 >        
616        }
617 <
617 >      
618   #endif
619 <
619 >      
620      } //end for(int i=0)  
621    }
622 <
397 <  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
398 <                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
399 <    int i;
400 <    char * id;
401 <    LinkedMolStamp* extractedStamp = NULL;
402 <    MoleculeStamp * currentStamp;
403 <    Component** the_components = simParams->getComponents();
404 <    int n_components = simParams->getNComponents();
405 <
406 <    if (!simParams->haveNMol()) {
407 <      // we don't have the total number of molecules, so we assume it is
408 <      // given in each component
409 <
410 <      for(i = 0; i < n_components; i++) {
411 <        if (!the_components[i]->haveNMol()) {
412 <          // we have a problem
413 <          sprintf(painCave.errMsg,
414 <                  "SimCreator Error. No global NMol or component NMol given.\n"
415 <                  "\tCannot calculate the number of atoms.\n");
416 <
417 <          painCave.isFatal = 1;
418 <          simError();
419 <        }
420 <      
421 <        id = the_components[i]->getType();
422 <
423 <        extractedStamp = stamps->extractMolStamp(id);
424 <        if (extractedStamp == NULL) {
425 <          sprintf(painCave.errMsg,
426 <                  "SimCreator error: Component \"%s\" was not found in the "
427 <                  "list of declared molecules\n", id);
428 <
429 <          painCave.isFatal = 1;
430 <          simError();
431 <        }
432 <
433 <        currentStamp = extractedStamp->getStamp();
434 <
435 <
436 <        moleculeStampPairs.push_back(
437 <                                     std::make_pair(currentStamp, the_components[i]->getNMol()));
438 <      } //end for (i = 0; i < n_components; i++)
439 <    } else {
440 <      sprintf(painCave.errMsg, "SimSetup error.\n"
441 <              "\tSorry, the ability to specify total"
442 <              " nMols and then give molfractions in the components\n"
443 <              "\tis not currently supported."
444 <              " Please give nMol in the components.\n");
445 <
446 <      painCave.isFatal = 1;
447 <      simError();
448 <    }
449 <
450 < #ifdef IS_MPI
451 <
452 <    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
453 <    MPIcheckPoint();
454 <
455 < #endif // is_mpi
456 <
457 <  }
458 <
622 >    
623    void SimCreator::setGlobalIndex(SimInfo *info) {
624      SimInfo::MoleculeIterator mi;
625      Molecule::AtomIterator ai;
626      Molecule::RigidBodyIterator ri;
627      Molecule::CutoffGroupIterator ci;
628 +    Molecule::IntegrableObjectIterator  ioi;
629      Molecule * mol;
630      Atom * atom;
631      RigidBody * rb;
# Line 469 | Line 634 | namespace oopse {
634      int beginRigidBodyIndex;
635      int beginCutoffGroupIndex;
636      int nGlobalAtoms = info->getNGlobalAtoms();
472 –    
473 – #ifndef IS_MPI
637  
638 +    /**@todo fixme */
639 + #ifndef IS_MPI
640 +    
641      beginAtomIndex = 0;
642      beginRigidBodyIndex = 0;
643      beginCutoffGroupIndex = 0;
644 <
644 >    
645   #else
646 <
646 >    
647      int nproc;
648      int myNode;
649 <
649 >    
650      myNode = worldRank;
651      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
652 <
652 >    
653      std::vector < int > tmpAtomsInProc(nproc, 0);
654      std::vector < int > tmpRigidBodiesInProc(nproc, 0);
655      std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
656      std::vector < int > NumAtomsInProc(nproc, 0);
657      std::vector < int > NumRigidBodiesInProc(nproc, 0);
658      std::vector < int > NumCutoffGroupsInProc(nproc, 0);
659 <
659 >    
660      tmpAtomsInProc[myNode] = info->getNAtoms();
661      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
662      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
663 <
663 >    
664      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
665      MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
666                    MPI_SUM, MPI_COMM_WORLD);
# Line 502 | Line 668 | namespace oopse {
668                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
669      MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
670                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
671 <
671 >    
672      beginAtomIndex = 0;
673      beginRigidBodyIndex = 0;
674      beginCutoffGroupIndex = 0;
675 <
675 >    
676      for(int i = 0; i < myNode; i++) {
677        beginAtomIndex += NumAtomsInProc[i];
678        beginRigidBodyIndex += NumRigidBodiesInProc[i];
679        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
680      }
681 <
681 >    
682   #endif
683 <
683 >    
684      //rigidbody's index begins right after atom's
685      beginRigidBodyIndex += info->getNGlobalAtoms();
686 <
686 >    
687      for(mol = info->beginMolecule(mi); mol != NULL;
688          mol = info->nextMolecule(mi)) {
689 <
689 >      
690        //local index(index in DataStorge) of atom is important
691        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
692          atom->setGlobalIndex(beginAtomIndex++);
693        }
694 <
694 >      
695        for(rb = mol->beginRigidBody(ri); rb != NULL;
696            rb = mol->nextRigidBody(ri)) {
697          rb->setGlobalIndex(beginRigidBodyIndex++);
698        }
699 <
699 >      
700        //local index of cutoff group is trivial, it only depends on the order of travesing
701        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
702            cg = mol->nextCutoffGroup(ci)) {
703          cg->setGlobalIndex(beginCutoffGroupIndex++);
704        }
705      }
706 <
706 >    
707      //fill globalGroupMembership
708      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
709      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
710        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
711 <
711 >        
712          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
713            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
714          }
715 <
715 >        
716        }      
717      }
718 <
718 >    
719   #ifdef IS_MPI    
720      // Since the globalGroupMembership has been zero filled and we've only
721      // poked values into the atoms we know, we can do an Allreduce
# Line 563 | Line 729 | namespace oopse {
729   #else
730      info->setGlobalGroupMembership(globalGroupMembership);
731   #endif
732 <
732 >    
733      //fill molMembership
734      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
735      
736      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
737 <
737 >      
738        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
739          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
740        }
741      }
742 <
742 >    
743   #ifdef IS_MPI
744      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
745 <
745 >    
746      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
747                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
748      
# Line 585 | Line 751 | namespace oopse {
751      info->setGlobalMolMembership(globalMolMembership);
752   #endif
753  
754 <  }
754 >    // nIOPerMol holds the number of integrable objects per molecule
755 >    // here the molecules are listed by their global indices.
756  
757 <  void SimCreator::loadCoordinates(SimInfo* info) {
757 >    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
758 >    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
759 >      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
760 >    }
761 >    
762 > #ifdef IS_MPI
763 >    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
764 >    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
765 >                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
766 > #else
767 >    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
768 > #endif    
769 >
770 > std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
771 >
772 > int startingIndex = 0;
773 > for (int i = 0; i < info->getNGlobalMolecules(); i++) {
774 >  startingIOIndexForMol[i] = startingIndex;
775 >  startingIndex += numIntegrableObjectsPerMol[i];
776 > }
777 >
778 > std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
779 > for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
780 >      int myGlobalIndex = mol->getGlobalIndex();
781 >      int globalIO = startingIOIndexForMol[myGlobalIndex];
782 >      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
783 >           integrableObject = mol->nextIntegrableObject(ioi)) {
784 >            integrableObject->setGlobalIntegrableObjectIndex(globalIO);
785 >            IOIndexToIntegrableObject[globalIO] = integrableObject;
786 >            globalIO++;
787 >      }
788 >    }
789 >
790 >  info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
791 >  
792 >  }
793 >  
794 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
795      Globals* simParams;
796      simParams = info->getSimParams();
797      
798 <    if (!simParams->haveInitialConfig()) {
799 <      sprintf(painCave.errMsg,
596 <              "Cannot intialize a simulation without an initial configuration file.\n");
597 <      painCave.isFatal = 1;;
598 <      simError();
599 <    }
600 <        
601 <    DumpReader reader(info, simParams->getInitialConfig());
798 >    
799 >    DumpReader reader(info, mdFileName);
800      int nframes = reader.getNFrames();
801 <
801 >    
802      if (nframes > 0) {
803        reader.readFrame(nframes - 1);
804      } else {
805        //invalid initial coordinate file
806 <      sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
807 <              simParams->getInitialConfig());
806 >      sprintf(painCave.errMsg,
807 >              "Initial configuration file %s should at least contain one frame\n",
808 >              mdFileName.c_str());
809        painCave.isFatal = 1;
810        simError();
811      }
812 <
812 >    
813      //copy the current snapshot to previous snapshot
814      info->getSnapshotManager()->advance();
815    }
816 <
816 >  
817   } //end namespace oopse
818  
819  

Diff Legend

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