ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/brains/SimCreator.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/brains/SimCreator.cpp (file contents):
Revision 2087 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
Revision 2544 by tim, Wed Jan 11 19:01:20 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(const std::string mdFileName){
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(mdFileName, 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(mdFileName);
126 >            lexer.initDeferredLineCount();
127 >    
128 >            // Create a parser that reads from the scanner
129 >            MDParser parser(lexer);
130 >            parser.setFilename(mdFileName);
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
152 <
153 <      yacc_BASS(mdFileName.c_str());
154 <
155 < #ifdef IS_MPI
156 <
157 <      throwMPIEvent(NULL);
158 <    } else {
159 <      set_interface_stamps(stamps, simParams);
160 <      mpiEventInit();
161 <      MPIcheckPoint();
162 <      mpiEventLoop();
163 <    }
164 <
165 < #endif
166 <
167 <  }
168 <
169 <  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
170 <    
171 <    MakeStamps * stamps = new MakeStamps();
172 <
173 <    Globals * simParams = new Globals();
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 <    //parse meta-data file
234 <    parseFile(mdFileName, stamps, simParams);
233 >        return simParams;
234 >  }
235 >  
236 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
237 >                                  bool loadInitCoords) {
238  
239 +    //parse meta-data file
240 +    Globals* simParams = parseFile(mdFileName);
241 +    
242      //create the force field
243 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
244 <                                                                         simParams->getForceField());
243 >    ForceField * ff = ForceFieldFactory::getInstance()
244 >      ->createForceField(simParams->getForceField());
245      
246      if (ff == NULL) {
247 <      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
248 <              simParams->getForceField());
247 >      sprintf(painCave.errMsg,
248 >              "ForceField Factory can not create %s force field\n",
249 >              simParams->getForceField().c_str());
250        painCave.isFatal = 1;
251        simError();
252      }
253 <
253 >    
254      if (simParams->haveForceFieldFileName()) {
255        ff->setForceFieldFileName(simParams->getForceFieldFileName());
256      }
257      
258      std::string forcefieldFileName;
259      forcefieldFileName = ff->getForceFieldFileName();
260 <
260 >    
261      if (simParams->haveForceFieldVariant()) {
262        //If the force field has variant, the variant force field name will be
263        //Base.variant.frc. For exampel EAM.u6.frc
264 <        
264 >      
265        std::string variant = simParams->getForceFieldVariant();
266 <
266 >      
267        std::string::size_type pos = forcefieldFileName.rfind(".frc");
268        variant = "." + variant;
269        if (pos != std::string::npos) {
# Line 139 | Line 275 | namespace oopse {
275      }
276      
277      ff->parse(forcefieldFileName);
278 <    
143 <    //extract the molecule stamps
144 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
145 <    compList(stamps, simParams, moleculeStampPairs);
146 <
278 >    ff->setFortranForceOptions();
279      //create SimInfo
280 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
281 <
280 >    SimInfo * info = new SimInfo(ff, simParams);
281 >    
282      //gather parameters (SimCreator only retrieves part of the parameters)
283      gatherParameters(info, mdFileName);
284 <
284 >    
285      //divide the molecules and determine the global index of molecules
286   #ifdef IS_MPI
287      divideMolecules(info);
288   #endif
289 <
289 >    
290      //create the molecules
291      createMolecules(info);
292 <
293 <
292 >    
293 >    
294      //allocate memory for DataStorage(circular reference, need to break it)
295      info->setSnapshotManager(new SimSnapshotManager(info));
296      
# Line 166 | Line 298 | namespace oopse {
298      //global index will never change again). Local indices of atoms and rigidbodies are already set by
299      //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
300      setGlobalIndex(info);
301 <
301 >    
302      //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
303      //atoms don't have the global index yet  (their global index are all initialized to -1).
304      //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
# Line 177 | Line 309 | namespace oopse {
309        info->addExcludePairs(mol);
310      }
311      
180 –
181 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
182 –    //eta, chi for NPT integrator)
312      if (loadInitCoords)
313        loadCoordinates(info);    
314      
315      return info;
316    }
317 <
317 >  
318    void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
319 <
320 <    //figure out the ouput file names
319 >    
320 >    //figure out the output file names
321      std::string prefix;
322 <
322 >    
323   #ifdef IS_MPI
324 <
324 >    
325      if (worldRank == 0) {
326   #endif // is_mpi
327        Globals * simParams = info->getSimParams();
# Line 201 | Line 330 | namespace oopse {
330        } else {
331          prefix = getPrefix(mdfile);
332        }
333 <
333 >      
334        info->setFinalConfigFileName(prefix + ".eor");
335        info->setDumpFileName(prefix + ".dump");
336        info->setStatFileName(prefix + ".stat");
337 <
337 >      info->setRestFileName(prefix + ".zang");
338 >      
339   #ifdef IS_MPI
340 <
340 >      
341      }
342 <
342 >    
343   #endif
344 <
344 >    
345    }
346 <
346 >  
347   #ifdef IS_MPI
348    void SimCreator::divideMolecules(SimInfo *info) {
349      double numerator;
# Line 237 | Line 367 | namespace oopse {
367      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
368      
369      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
370 <
370 >    
371      if (nProcessors > nGlobalMols) {
372        sprintf(painCave.errMsg,
373                "nProcessors (%d) > nMol (%d)\n"
# Line 246 | Line 376 | namespace oopse {
376                "\tusable division of atoms for force decomposition.\n"
377                "\tEither try a smaller number of processors, or run the\n"
378                "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
379 <
379 >      
380        painCave.isFatal = 1;
381        simError();
382      }
383 <
383 >    
384      int seedValue;
385      Globals * simParams = info->getSimParams();
386      SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
# Line 260 | Line 390 | namespace oopse {
390      }else {
391        myRandom = new SeqRandNumGen();
392      }  
393 <
394 <
393 >    
394 >    
395      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
396 <
396 >    
397      //initialize atomsPerProc
398      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
399 <
399 >    
400      if (worldRank == 0) {
401        numerator = info->getNGlobalAtoms();
402        denominator = nProcessors;
403        precast = numerator / denominator;
404        nTarget = (int)(precast + 0.5);
405 <
405 >      
406        for(i = 0; i < nGlobalMols; i++) {
407          done = 0;
408          loops = 0;
409 <
409 >        
410          while (!done) {
411            loops++;
412 <
412 >          
413            // Pick a processor at random
414 <
414 >          
415            which_proc = (int) (myRandom->rand() * nProcessors);
416 <
416 >          
417            //get the molecule stamp first
418            int stampId = info->getMoleculeStampId(i);
419            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
420 <
420 >          
421            // How many atoms does this processor have so far?
422            old_atoms = atomsPerProc[which_proc];
423            add_atoms = moleculeStamp->getNAtoms();
424            new_atoms = old_atoms + add_atoms;
425 <
425 >          
426            // If we've been through this loop too many times, we need
427            // to just give up and assign the molecule to this processor
428            // and be done with it.
429 <
429 >          
430            if (loops > 100) {
431              sprintf(painCave.errMsg,
432                      "I've tried 100 times to assign molecule %d to a "
433                      " processor, but can't find a good spot.\n"
434                      "I'm assigning it at random to processor %d.\n",
435                      i, which_proc);
436 <
436 >            
437              painCave.isFatal = 0;
438              simError();
439 <
439 >            
440              molToProcMap[i] = which_proc;
441              atomsPerProc[which_proc] += add_atoms;
442 <
442 >            
443              done = 1;
444              continue;
445            }
446 <
446 >          
447            // If we can add this molecule to this processor without sending
448            // it above nTarget, then go ahead and do it:
449 <
449 >          
450            if (new_atoms <= nTarget) {
451              molToProcMap[i] = which_proc;
452              atomsPerProc[which_proc] += add_atoms;
453 <
453 >            
454              done = 1;
455              continue;
456            }
457 <
457 >          
458            // The only situation left is when new_atoms > nTarget.  We
459            // want to accept this with some probability that dies off the
460            // farther we are from nTarget
461 <
461 >          
462            // roughly:  x = new_atoms - nTarget
463            //           Pacc(x) = exp(- a * x)
464            // where a = penalty / (average atoms per molecule)
465 <
465 >          
466            x = (double)(new_atoms - nTarget);
467            y = myRandom->rand();
468 <
468 >          
469            if (y < exp(- a * x)) {
470              molToProcMap[i] = which_proc;
471              atomsPerProc[which_proc] += add_atoms;
472 <
472 >            
473              done = 1;
474              continue;
475            } else {
# Line 347 | Line 477 | namespace oopse {
477            }
478          }
479        }
480 <
480 >      
481        delete myRandom;
482 <        
482 >      
483        // Spray out this nonsense to all other processors:
484 <
484 >      
485        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
486      } else {
487 <
487 >      
488        // Listen to your marching orders from processor 0:
489 <
489 >      
490        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
491      }
492 <
492 >    
493      info->setMolToProcMap(molToProcMap);
494      sprintf(checkPointMsg,
495              "Successfully divided the molecules among the processors.\n");
496      MPIcheckPoint();
497    }
498 <
498 >  
499   #endif
500 <
500 >  
501    void SimCreator::createMolecules(SimInfo *info) {
502      MoleculeCreator molCreator;
503      int stampId;
504 <
504 >    
505      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
506 <
506 >      
507   #ifdef IS_MPI
508 <
508 >      
509        if (info->getMolToProc(i) == worldRank) {
510   #endif
511 <
511 >        
512          stampId = info->getMoleculeStampId(i);
513          Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
514                                                     stampId, i, info->getLocalIndexManager());
515 <
515 >        
516          info->addMolecule(mol);
517 <
517 >        
518   #ifdef IS_MPI
519 <
519 >        
520        }
521 <
521 >      
522   #endif
523 <
523 >      
524      } //end for(int i=0)  
525    }
526 <
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 <
526 >    
527    void SimCreator::setGlobalIndex(SimInfo *info) {
528      SimInfo::MoleculeIterator mi;
529      Molecule::AtomIterator ai;
# Line 471 | Line 539 | namespace oopse {
539      int nGlobalAtoms = info->getNGlobalAtoms();
540      
541   #ifndef IS_MPI
542 <
542 >    
543      beginAtomIndex = 0;
544      beginRigidBodyIndex = 0;
545      beginCutoffGroupIndex = 0;
546 <
546 >    
547   #else
548 <
548 >    
549      int nproc;
550      int myNode;
551 <
551 >    
552      myNode = worldRank;
553      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
554 <
554 >    
555      std::vector < int > tmpAtomsInProc(nproc, 0);
556      std::vector < int > tmpRigidBodiesInProc(nproc, 0);
557      std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
558      std::vector < int > NumAtomsInProc(nproc, 0);
559      std::vector < int > NumRigidBodiesInProc(nproc, 0);
560      std::vector < int > NumCutoffGroupsInProc(nproc, 0);
561 <
561 >    
562      tmpAtomsInProc[myNode] = info->getNAtoms();
563      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
564      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
565 <
565 >    
566      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
567      MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
568                    MPI_SUM, MPI_COMM_WORLD);
# Line 502 | Line 570 | namespace oopse {
570                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
571      MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
572                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
573 <
573 >    
574      beginAtomIndex = 0;
575      beginRigidBodyIndex = 0;
576      beginCutoffGroupIndex = 0;
577 <
577 >    
578      for(int i = 0; i < myNode; i++) {
579        beginAtomIndex += NumAtomsInProc[i];
580        beginRigidBodyIndex += NumRigidBodiesInProc[i];
581        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
582      }
583 <
583 >    
584   #endif
585 <
585 >    
586      //rigidbody's index begins right after atom's
587      beginRigidBodyIndex += info->getNGlobalAtoms();
588 <
588 >    
589      for(mol = info->beginMolecule(mi); mol != NULL;
590          mol = info->nextMolecule(mi)) {
591 <
591 >      
592        //local index(index in DataStorge) of atom is important
593        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
594          atom->setGlobalIndex(beginAtomIndex++);
595        }
596 <
596 >      
597        for(rb = mol->beginRigidBody(ri); rb != NULL;
598            rb = mol->nextRigidBody(ri)) {
599          rb->setGlobalIndex(beginRigidBodyIndex++);
600        }
601 <
601 >      
602        //local index of cutoff group is trivial, it only depends on the order of travesing
603        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
604            cg = mol->nextCutoffGroup(ci)) {
605          cg->setGlobalIndex(beginCutoffGroupIndex++);
606        }
607      }
608 <
608 >    
609      //fill globalGroupMembership
610      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
611      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
612        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
613 <
613 >        
614          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
615            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
616          }
617 <
617 >        
618        }      
619      }
620 <
620 >    
621   #ifdef IS_MPI    
622      // Since the globalGroupMembership has been zero filled and we've only
623      // poked values into the atoms we know, we can do an Allreduce
# Line 563 | Line 631 | namespace oopse {
631   #else
632      info->setGlobalGroupMembership(globalGroupMembership);
633   #endif
634 <
634 >    
635      //fill molMembership
636      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
637      
638      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
639 <
639 >      
640        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
641          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
642        }
643      }
644 <
644 >    
645   #ifdef IS_MPI
646      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
647 <
647 >    
648      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
649                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
650      
# Line 584 | Line 652 | namespace oopse {
652   #else
653      info->setGlobalMolMembership(globalMolMembership);
654   #endif
655 <
655 >    
656    }
657 <
657 >  
658    void SimCreator::loadCoordinates(SimInfo* info) {
659      Globals* simParams;
660      simParams = info->getSimParams();
# Line 597 | Line 665 | namespace oopse {
665        painCave.isFatal = 1;;
666        simError();
667      }
668 <        
668 >    
669      DumpReader reader(info, simParams->getInitialConfig());
670      int nframes = reader.getNFrames();
671 <
671 >    
672      if (nframes > 0) {
673        reader.readFrame(nframes - 1);
674      } else {
675        //invalid initial coordinate file
676 <      sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
677 <              simParams->getInitialConfig());
676 >      sprintf(painCave.errMsg,
677 >              "Initial configuration file %s should at least contain one frame\n",
678 >              simParams->getInitialConfig().c_str());
679        painCave.isFatal = 1;
680        simError();
681      }
682 <
682 >    
683      //copy the current snapshot to previous snapshot
684      info->getSnapshotManager()->advance();
685    }
686 <
686 >  
687   } //end namespace oopse
688  
689  

Diff Legend

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