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 1051 by gezelter, Mon Sep 25 22:08:33 2006 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 <      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();
90 <      MPIcheckPoint();
91 <      mpiEventLoop();
92 <    }
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
105 <    parseFile(mdFileName, stamps, simParams);
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()->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      }
343 <
343 >    
344      if (simParams->haveForceFieldFileName()) {
345 +      std::cout<< simParams->getForceFieldFileName() << "\n";
346        ff->setForceFieldFileName(simParams->getForceFieldFileName());
347      }
348      
349      std::string forcefieldFileName;
350      forcefieldFileName = ff->getForceFieldFileName();
351 <
351 >    
352      if (simParams->haveForceFieldVariant()) {
353        //If the force field has variant, the variant force field name will be
354        //Base.variant.frc. For exampel EAM.u6.frc
355 <        
355 >      
356        std::string variant = simParams->getForceFieldVariant();
357 <
357 >      
358        std::string::size_type pos = forcefieldFileName.rfind(".frc");
359        variant = "." + variant;
360        if (pos != std::string::npos) {
# Line 139 | Line 366 | namespace oopse {
366      }
367      
368      ff->parse(forcefieldFileName);
369 <    
143 <    //extract the molecule stamps
144 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
145 <    compList(stamps, simParams, moleculeStampPairs);
146 <
369 >    ff->setFortranForceOptions();
370      //create SimInfo
371 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
371 >    SimInfo * info = new SimInfo(ff, simParams);
372  
373 <    //gather parameters (SimCreator only retrieves part of the parameters)
373 >    info->setRawMetaData(mdRawData);
374 >    
375 >    //gather parameters (SimCreator only retrieves part of the
376 >    //parameters)
377      gatherParameters(info, mdFileName);
378 <
378 >    
379      //divide the molecules and determine the global index of molecules
380   #ifdef IS_MPI
381      divideMolecules(info);
382   #endif
383 <
383 >    
384      //create the molecules
385      createMolecules(info);
386 <
387 <
388 <    //allocate memory for DataStorage(circular reference, need to break it)
386 >    
387 >    
388 >    //allocate memory for DataStorage(circular reference, need to
389 >    //break it)
390      info->setSnapshotManager(new SimSnapshotManager(info));
391      
392 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
393 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
394 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
392 >    //set the global index of atoms, rigidbodies and cutoffgroups
393 >    //(only need to be set once, the global index will never change
394 >    //again). Local indices of atoms and rigidbodies are already set
395 >    //by MoleculeCreator class which actually delegates the
396 >    //responsibility to LocalIndexManager.
397      setGlobalIndex(info);
398 <
399 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
400 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
401 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
402 <    //we can determine the beginning global indices of atoms before they get created.
398 >    
399 >    //Although addExcludePairs is called inside SimInfo's addMolecule
400 >    //method, at that point atoms don't have the global index yet
401 >    //(their global index are all initialized to -1).  Therefore we
402 >    //have to call addExcludePairs explicitly here. A way to work
403 >    //around is that we can determine the beginning global indices of
404 >    //atoms before they get created.
405      SimInfo::MoleculeIterator mi;
406      Molecule* mol;
407      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
408        info->addExcludePairs(mol);
409      }
410      
180 –
181 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
182 –    //eta, chi for NPT integrator)
411      if (loadInitCoords)
412 <      loadCoordinates(info);    
412 >      loadCoordinates(info, mdFileName);    
413      
414      return info;
415    }
416 <
416 >  
417    void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
418 <
419 <    //figure out the ouput file names
418 >    
419 >    //figure out the output file names
420      std::string prefix;
421 <
421 >    
422   #ifdef IS_MPI
423 <
423 >    
424      if (worldRank == 0) {
425   #endif // is_mpi
426        Globals * simParams = info->getSimParams();
# Line 201 | Line 429 | namespace oopse {
429        } else {
430          prefix = getPrefix(mdfile);
431        }
432 <
432 >      
433        info->setFinalConfigFileName(prefix + ".eor");
434        info->setDumpFileName(prefix + ".dump");
435        info->setStatFileName(prefix + ".stat");
436 <
436 >      info->setRestFileName(prefix + ".zang");
437 >      
438   #ifdef IS_MPI
439 <
439 >      
440      }
441 <
441 >    
442   #endif
443 <
443 >    
444    }
445 <
445 >  
446   #ifdef IS_MPI
447    void SimCreator::divideMolecules(SimInfo *info) {
448 <    double numerator;
449 <    double denominator;
450 <    double precast;
451 <    double x;
452 <    double y;
453 <    double a;
448 >    RealType numerator;
449 >    RealType denominator;
450 >    RealType precast;
451 >    RealType x;
452 >    RealType y;
453 >    RealType a;
454      int old_atoms;
455      int add_atoms;
456      int new_atoms;
# Line 237 | Line 466 | namespace oopse {
466      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
467      
468      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
469 <
469 >    
470      if (nProcessors > nGlobalMols) {
471        sprintf(painCave.errMsg,
472                "nProcessors (%d) > nMol (%d)\n"
# Line 246 | Line 475 | namespace oopse {
475                "\tusable division of atoms for force decomposition.\n"
476                "\tEither try a smaller number of processors, or run the\n"
477                "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
478 <
478 >      
479        painCave.isFatal = 1;
480        simError();
481      }
482 <
482 >    
483      int seedValue;
484      Globals * simParams = info->getSimParams();
485      SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
# Line 260 | Line 489 | namespace oopse {
489      }else {
490        myRandom = new SeqRandNumGen();
491      }  
492 <
493 <
492 >    
493 >    
494      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
495 <
495 >    
496      //initialize atomsPerProc
497      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
498 <
498 >    
499      if (worldRank == 0) {
500        numerator = info->getNGlobalAtoms();
501        denominator = nProcessors;
502        precast = numerator / denominator;
503        nTarget = (int)(precast + 0.5);
504 <
504 >      
505        for(i = 0; i < nGlobalMols; i++) {
506          done = 0;
507          loops = 0;
508 <
508 >        
509          while (!done) {
510            loops++;
511 <
511 >          
512            // Pick a processor at random
513 <
513 >          
514            which_proc = (int) (myRandom->rand() * nProcessors);
515 <
515 >          
516            //get the molecule stamp first
517            int stampId = info->getMoleculeStampId(i);
518            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
519 <
519 >          
520            // How many atoms does this processor have so far?
521            old_atoms = atomsPerProc[which_proc];
522            add_atoms = moleculeStamp->getNAtoms();
523            new_atoms = old_atoms + add_atoms;
524 <
524 >          
525            // If we've been through this loop too many times, we need
526            // to just give up and assign the molecule to this processor
527            // and be done with it.
528 <
528 >          
529            if (loops > 100) {
530              sprintf(painCave.errMsg,
531                      "I've tried 100 times to assign molecule %d to a "
532                      " processor, but can't find a good spot.\n"
533                      "I'm assigning it at random to processor %d.\n",
534                      i, which_proc);
535 <
535 >            
536              painCave.isFatal = 0;
537              simError();
538 <
538 >            
539              molToProcMap[i] = which_proc;
540              atomsPerProc[which_proc] += add_atoms;
541 <
541 >            
542              done = 1;
543              continue;
544            }
545 <
545 >          
546            // If we can add this molecule to this processor without sending
547            // it above nTarget, then go ahead and do it:
548 <
548 >          
549            if (new_atoms <= nTarget) {
550              molToProcMap[i] = which_proc;
551              atomsPerProc[which_proc] += add_atoms;
552 <
552 >            
553              done = 1;
554              continue;
555            }
556 <
556 >          
557            // The only situation left is when new_atoms > nTarget.  We
558            // want to accept this with some probability that dies off the
559            // farther we are from nTarget
560 <
560 >          
561            // roughly:  x = new_atoms - nTarget
562            //           Pacc(x) = exp(- a * x)
563            // where a = penalty / (average atoms per molecule)
564 <
565 <          x = (double)(new_atoms - nTarget);
564 >          
565 >          x = (RealType)(new_atoms - nTarget);
566            y = myRandom->rand();
567 <
567 >          
568            if (y < exp(- a * x)) {
569              molToProcMap[i] = which_proc;
570              atomsPerProc[which_proc] += add_atoms;
571 <
571 >            
572              done = 1;
573              continue;
574            } else {
# Line 347 | Line 576 | namespace oopse {
576            }
577          }
578        }
579 <
579 >      
580        delete myRandom;
581 <        
581 >      
582        // Spray out this nonsense to all other processors:
583 <
583 >      
584        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
585      } else {
586 <
586 >      
587        // Listen to your marching orders from processor 0:
588 <
588 >      
589        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
590      }
591 <
591 >    
592      info->setMolToProcMap(molToProcMap);
593      sprintf(checkPointMsg,
594              "Successfully divided the molecules among the processors.\n");
595      MPIcheckPoint();
596    }
597 <
597 >  
598   #endif
599 <
599 >  
600    void SimCreator::createMolecules(SimInfo *info) {
601      MoleculeCreator molCreator;
602      int stampId;
603 <
603 >    
604      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
605 <
605 >      
606   #ifdef IS_MPI
607 <
607 >      
608        if (info->getMolToProc(i) == worldRank) {
609   #endif
610 <
610 >        
611          stampId = info->getMoleculeStampId(i);
612          Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
613                                                     stampId, i, info->getLocalIndexManager());
614 <
614 >        
615          info->addMolecule(mol);
616 <
616 >        
617   #ifdef IS_MPI
618 <
618 >        
619        }
620 <
620 >      
621   #endif
622 <
622 >      
623      } //end for(int i=0)  
624    }
625 <
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 <
625 >    
626    void SimCreator::setGlobalIndex(SimInfo *info) {
627      SimInfo::MoleculeIterator mi;
628      Molecule::AtomIterator ai;
629      Molecule::RigidBodyIterator ri;
630      Molecule::CutoffGroupIterator ci;
631 +    Molecule::IntegrableObjectIterator  ioi;
632      Molecule * mol;
633      Atom * atom;
634      RigidBody * rb;
# Line 469 | Line 637 | namespace oopse {
637      int beginRigidBodyIndex;
638      int beginCutoffGroupIndex;
639      int nGlobalAtoms = info->getNGlobalAtoms();
472 –    
473 – #ifndef IS_MPI
640  
641 +    /**@todo fixme */
642 + #ifndef IS_MPI
643 +    
644      beginAtomIndex = 0;
645      beginRigidBodyIndex = 0;
646      beginCutoffGroupIndex = 0;
647 <
647 >    
648   #else
649 <
649 >    
650      int nproc;
651      int myNode;
652 <
652 >    
653      myNode = worldRank;
654      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
655 <
655 >    
656      std::vector < int > tmpAtomsInProc(nproc, 0);
657      std::vector < int > tmpRigidBodiesInProc(nproc, 0);
658      std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
659      std::vector < int > NumAtomsInProc(nproc, 0);
660      std::vector < int > NumRigidBodiesInProc(nproc, 0);
661      std::vector < int > NumCutoffGroupsInProc(nproc, 0);
662 <
662 >    
663      tmpAtomsInProc[myNode] = info->getNAtoms();
664      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
665      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
666 <
666 >    
667      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
668      MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
669                    MPI_SUM, MPI_COMM_WORLD);
# Line 502 | Line 671 | namespace oopse {
671                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
672      MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
673                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
674 <
674 >    
675      beginAtomIndex = 0;
676      beginRigidBodyIndex = 0;
677      beginCutoffGroupIndex = 0;
678 <
678 >    
679      for(int i = 0; i < myNode; i++) {
680        beginAtomIndex += NumAtomsInProc[i];
681        beginRigidBodyIndex += NumRigidBodiesInProc[i];
682        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
683      }
684 <
684 >    
685   #endif
686 <
686 >    
687      //rigidbody's index begins right after atom's
688      beginRigidBodyIndex += info->getNGlobalAtoms();
689 <
689 >    
690      for(mol = info->beginMolecule(mi); mol != NULL;
691          mol = info->nextMolecule(mi)) {
692 <
692 >      
693        //local index(index in DataStorge) of atom is important
694        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
695          atom->setGlobalIndex(beginAtomIndex++);
696        }
697 <
697 >      
698        for(rb = mol->beginRigidBody(ri); rb != NULL;
699            rb = mol->nextRigidBody(ri)) {
700          rb->setGlobalIndex(beginRigidBodyIndex++);
701        }
702 <
702 >      
703        //local index of cutoff group is trivial, it only depends on the order of travesing
704        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
705            cg = mol->nextCutoffGroup(ci)) {
706          cg->setGlobalIndex(beginCutoffGroupIndex++);
707        }
708      }
709 <
709 >    
710      //fill globalGroupMembership
711      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
712      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
713        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
714 <
714 >        
715          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
716            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
717          }
718 <
718 >        
719        }      
720      }
721 <
721 >    
722   #ifdef IS_MPI    
723      // Since the globalGroupMembership has been zero filled and we've only
724      // poked values into the atoms we know, we can do an Allreduce
# Line 563 | Line 732 | namespace oopse {
732   #else
733      info->setGlobalGroupMembership(globalGroupMembership);
734   #endif
735 <
735 >    
736      //fill molMembership
737      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
738      
739      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
740 <
740 >      
741        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
742          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
743        }
744      }
745 <
745 >    
746   #ifdef IS_MPI
747      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
748 <
748 >    
749      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
750                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
751      
# Line 585 | Line 754 | namespace oopse {
754      info->setGlobalMolMembership(globalMolMembership);
755   #endif
756  
757 <  }
757 >    // nIOPerMol holds the number of integrable objects per molecule
758 >    // here the molecules are listed by their global indices.
759  
760 <  void SimCreator::loadCoordinates(SimInfo* info) {
760 >    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
761 >    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
762 >      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
763 >    }
764 >    
765 > #ifdef IS_MPI
766 >    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
767 >    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
768 >                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
769 > #else
770 >    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
771 > #endif    
772 >
773 > std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
774 >
775 > int startingIndex = 0;
776 > for (int i = 0; i < info->getNGlobalMolecules(); i++) {
777 >  startingIOIndexForMol[i] = startingIndex;
778 >  startingIndex += numIntegrableObjectsPerMol[i];
779 > }
780 >
781 > std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
782 > for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
783 >      int myGlobalIndex = mol->getGlobalIndex();
784 >      int globalIO = startingIOIndexForMol[myGlobalIndex];
785 >      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
786 >           integrableObject = mol->nextIntegrableObject(ioi)) {
787 >            integrableObject->setGlobalIntegrableObjectIndex(globalIO);
788 >            IOIndexToIntegrableObject[globalIO] = integrableObject;
789 >            globalIO++;
790 >      }
791 >    }
792 >
793 >  info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
794 >  
795 >  }
796 >  
797 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
798      Globals* simParams;
799      simParams = info->getSimParams();
800      
801 <    if (!simParams->haveInitialConfig()) {
802 <      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());
801 >    
802 >    DumpReader reader(info, mdFileName);
803      int nframes = reader.getNFrames();
804 <
804 >    
805      if (nframes > 0) {
806        reader.readFrame(nframes - 1);
807      } else {
808        //invalid initial coordinate file
809 <      sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
810 <              simParams->getInitialConfig());
809 >      sprintf(painCave.errMsg,
810 >              "Initial configuration file %s should at least contain one frame\n",
811 >              mdFileName.c_str());
812        painCave.isFatal = 1;
813        simError();
814      }
815 <
815 >    
816      //copy the current snapshot to previous snapshot
817      info->getSnapshotManager()->advance();
818    }
819 <
819 >  
820   } //end namespace oopse
821  
822  

Diff Legend

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