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 1981 by tim, Mon Feb 7 19:14:26 2005 UTC vs.
Revision 2469 by tim, Fri Dec 2 15:38:03 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 47 | Line 47
47   * @version 1.0
48   */
49  
50 < #include <sprng.h>
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"
56 – #include "io/parse_me.h"
58   #include "UseTheForce/ForceFieldFactory.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61 < #ifdef IS_MPI
62 < #include "io/mpiBASS.h"
63 < #include "math/randomSPRNG.hpp"
64 < #endif
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  
65 – namespace oopse {
67  
67 – void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){
68 –
68   #ifdef IS_MPI
69 <
71 <    if (worldRank == 0) {
72 < #endif // is_mpi
73 <
74 <        simParams->initalize();
75 <        set_interface_stamps(stamps, simParams);
76 <
77 < #ifdef IS_MPI
78 <
79 <        mpiEventInit();
80 <
69 > #include "math/ParallelRandNumGen.hpp"
70   #endif
71  
72 <        yacc_BASS(mdFileName.c_str());
72 > namespace oopse {
73 >  
74 > Globals* SimCreator::parseFile(const std::string mdFileName){
75 >        Globals* simParams = NULL;
76 >        try {
77  
78 < #ifdef IS_MPI
78 >            // Create a preprocessor that preprocesses md file into an ostringstream
79 >            std::stringstream ppStream;
80 > #ifdef IS_MPI            
81 >            int streamSize;
82 >            const int masterNode = 0;
83 >            int commStatus;
84 >            if (worldRank == masterNode) {
85 > #endif
86 >                
87 >                SimplePreprocessor preprocessor;
88 >                preprocessor.preprocess(mdFileName, ppStream);
89 >                
90 > #ifdef IS_MPI            
91 >                //brocasting the stream size
92 >                streamSize = ppStream.str().size() +1;
93 >                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
94  
95 <        throwMPIEvent(NULL);
96 <    } else {
97 <        set_interface_stamps(stamps, simParams);
98 <        mpiEventInit();
99 <        MPIcheckPoint();
100 <        mpiEventLoop();
101 <    }
95 >                commStatus = MPI_Bcast(ppStream.str().c_str(), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
96 >            
97 >                
98 >            } else {
99 >                //get stream size
100 >                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
101 >                
102 >                  char* buf = new char[streamSize];
103 >                  assert(buf);
104 >                
105 >                  //receive file content
106 >                  commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
107 >                
108 >                  ppStream.str(buf);
109 >                  delete buf;
110  
111 < #endif
111 >            }
112 > #endif            
113 >            // Create a scanner that reads from the input stream
114 >            MDLexer lexer(ppStream);
115 >            lexer.setFilename(mdFileName);
116 >            lexer.initDeferredLineCount();
117 >    
118 >            // Create a parser that reads from the scanner
119 >            MDParser parser(lexer);
120 >            parser.setFilename(mdFileName);
121  
122 < }
123 <
124 < SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
122 >            // Create an observer that synchorizes file name change
123 >            FilenameObserver observer;
124 >            observer.setLexer(&lexer);
125 >            observer.setParser(&parser);
126 >            lexer.setObserver(&observer);
127      
128 <    MakeStamps * stamps = new MakeStamps();
128 >            antlr::ASTFactory factory;
129 >            parser.initializeASTFactory(factory);
130 >            parser.setASTFactory(&factory);
131 >            parser.mdfile();
132  
133 <    Globals * simParams = new Globals();
133 >            // Create a tree parser that reads information into Globals
134 >            MDTreeParser treeParser;
135 >            treeParser.initializeASTFactory(factory);
136 >            treeParser.setASTFactory(&factory);
137 >             simParams = treeParser.walkTree(parser.getAST());
138  
139 <    //parse meta-data file
140 <    parseFile(mdFileName, stamps, simParams);
139 >        }
140 >        catch (exception& e) {
141 >            cerr << "parser exception: " << e.what() << endl;
142 >        }
143  
144 +        return simParams;
145 +  }
146 +  
147 +  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
148 +                                  bool loadInitCoords) {
149 +
150 +    //parse meta-data file
151 +    Globals* simParams = parseFile(mdFileName);
152 +    
153      //create the force field
154 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
155 <                          simParams->getForceField());
154 >    ForceField * ff = ForceFieldFactory::getInstance()
155 >      ->createForceField(simParams->getForceField());
156      
157      if (ff == NULL) {
158 <        sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
159 <                simParams->getForceField());
160 <        painCave.isFatal = 1;
161 <        simError();
158 >      sprintf(painCave.errMsg,
159 >              "ForceField Factory can not create %s force field\n",
160 >              simParams->getForceField().c_str());
161 >      painCave.isFatal = 1;
162 >      simError();
163      }
164 <
164 >    
165      if (simParams->haveForceFieldFileName()) {
166 <        ff->setForceFieldFileName(simParams->getForceFieldFileName());
166 >      ff->setForceFieldFileName(simParams->getForceFieldFileName());
167      }
168      
169      std::string forcefieldFileName;
170      forcefieldFileName = ff->getForceFieldFileName();
171 <
171 >    
172      if (simParams->haveForceFieldVariant()) {
173 <        //If the force field has variant, the variant force field name will be
174 <        //Base.variant.frc. For exampel EAM.u6.frc
175 <        
176 <        std::string variant = simParams->getForceFieldVariant();
177 <
178 <        std::string::size_type pos = forcefieldFileName.rfind(".frc");
179 <        variant = "." + variant;
180 <        if (pos != std::string::npos) {
181 <            forcefieldFileName.insert(pos, variant);
182 <        } else {
183 <            //If the default force field file name does not containt .frc suffix, just append the .variant
184 <            forcefieldFileName.append(variant);
185 <        }
173 >      //If the force field has variant, the variant force field name will be
174 >      //Base.variant.frc. For exampel EAM.u6.frc
175 >      
176 >      std::string variant = simParams->getForceFieldVariant();
177 >      
178 >      std::string::size_type pos = forcefieldFileName.rfind(".frc");
179 >      variant = "." + variant;
180 >      if (pos != std::string::npos) {
181 >        forcefieldFileName.insert(pos, variant);
182 >      } else {
183 >        //If the default force field file name does not containt .frc suffix, just append the .variant
184 >        forcefieldFileName.append(variant);
185 >      }
186      }
187      
188      ff->parse(forcefieldFileName);
189 <    
144 <    //extract the molecule stamps
145 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
146 <    compList(stamps, simParams, moleculeStampPairs);
147 <
189 >        
190      //create SimInfo
191 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
192 <
191 >    SimInfo * info = new SimInfo(ff, simParams);
192 >    
193      //gather parameters (SimCreator only retrieves part of the parameters)
194      gatherParameters(info, mdFileName);
195 <
195 >    
196      //divide the molecules and determine the global index of molecules
197   #ifdef IS_MPI
198      divideMolecules(info);
199   #endif
200 <
200 >    
201      //create the molecules
202      createMolecules(info);
203 <
204 <
203 >    
204 >    
205      //allocate memory for DataStorage(circular reference, need to break it)
206      info->setSnapshotManager(new SimSnapshotManager(info));
207      
# Line 167 | Line 209 | SimInfo*  SimCreator::createSim(const std::string & md
209      //global index will never change again). Local indices of atoms and rigidbodies are already set by
210      //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
211      setGlobalIndex(info);
212 <
212 >    
213      //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
214      //atoms don't have the global index yet  (their global index are all initialized to -1).
215      //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
# Line 175 | Line 217 | SimInfo*  SimCreator::createSim(const std::string & md
217      SimInfo::MoleculeIterator mi;
218      Molecule* mol;
219      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
220 <        info->addExcludePairs(mol);
220 >      info->addExcludePairs(mol);
221      }
222      
181 –
182 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
183 –    //eta, chi for NPT integrator)
223      if (loadInitCoords)
224 <        loadCoordinates(info);    
224 >      loadCoordinates(info);    
225      
226      return info;
227 < }
228 <
229 < void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
230 <
231 <    //setup seed for random number generator
193 <    int seedValue;
194 <    Globals * simParams = info->getSimParams();
195 <
196 <    if (simParams->haveSeed()) {
197 <        seedValue = simParams->getSeed();
198 <
199 <        if (seedValue < 100000000 ) {
200 <            sprintf(painCave.errMsg,
201 <                    "Seed for sprng library should contain at least 9 digits\n"
202 <                        "OOPSE will generate a seed for user\n");
203 <
204 <            painCave.isFatal = 0;
205 <            simError();
206 <
207 <            //using seed generated by system instead of invalid seed set by user
208 <
209 < #ifndef IS_MPI
210 <
211 <            seedValue = make_sprng_seed();
212 <
213 < #else
214 <
215 <            if (worldRank == 0) {
216 <                seedValue = make_sprng_seed();
217 <            }
218 <
219 <            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
220 <
221 < #endif
222 <
223 <        } //end if (seedValue /1000000000 == 0)
224 <    } else {
225 <
226 < #ifndef IS_MPI
227 <
228 <        seedValue = make_sprng_seed();
229 <
230 < #else
231 <
232 <        if (worldRank == 0) {
233 <            seedValue = make_sprng_seed();
234 <        }
235 <
236 <        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
237 <
238 < #endif
239 <
240 <    } //end of simParams->haveSeed()
241 <
242 <    info->setSeed(seedValue);
243 <
244 <
245 <    //figure out the ouput file names
227 >  }
228 >  
229 >  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
230 >    
231 >    //figure out the output file names
232      std::string prefix;
233 <
233 >    
234   #ifdef IS_MPI
235 <
235 >    
236      if (worldRank == 0) {
237   #endif // is_mpi
238 <
239 <        if (simParams->haveFinalConfig()) {
240 <            prefix = getPrefix(simParams->getFinalConfig());
241 <        } else {
242 <            prefix = getPrefix(mdfile);
243 <        }
244 <
245 <        info->setFinalConfigFileName(prefix + ".eor");
246 <        info->setDumpFileName(prefix + ".dump");
247 <        info->setStatFileName(prefix + ".stat");
248 <
238 >      Globals * simParams = info->getSimParams();
239 >      if (simParams->haveFinalConfig()) {
240 >        prefix = getPrefix(simParams->getFinalConfig());
241 >      } else {
242 >        prefix = getPrefix(mdfile);
243 >      }
244 >      
245 >      info->setFinalConfigFileName(prefix + ".eor");
246 >      info->setDumpFileName(prefix + ".dump");
247 >      info->setStatFileName(prefix + ".stat");
248 >      info->setRestFileName(prefix + ".zang");
249 >      
250   #ifdef IS_MPI
251 <
251 >      
252      }
253 <
253 >    
254   #endif
255 <
256 < }
257 <
255 >    
256 >  }
257 >  
258   #ifdef IS_MPI
259 < void SimCreator::divideMolecules(SimInfo *info) {
259 >  void SimCreator::divideMolecules(SimInfo *info) {
260      double numerator;
261      double denominator;
262      double precast;
# Line 287 | Line 274 | void SimCreator::divideMolecules(SimInfo *info) {
274      int which_proc;
275      int nProcessors;
276      std::vector<int> atomsPerProc;
290 –    randomSPRNG myRandom(info->getSeed());
277      int nGlobalMols = info->getNGlobalMolecules();
278      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
279      
280      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
281 <
281 >    
282      if (nProcessors > nGlobalMols) {
283 <        sprintf(painCave.errMsg,
284 <                "nProcessors (%d) > nMol (%d)\n"
285 <                    "\tThe number of processors is larger than\n"
286 <                    "\tthe number of molecules.  This will not result in a \n"
287 <                    "\tusable division of atoms for force decomposition.\n"
288 <                    "\tEither try a smaller number of processors, or run the\n"
289 <                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
290 <
291 <        painCave.isFatal = 1;
292 <        simError();
283 >      sprintf(painCave.errMsg,
284 >              "nProcessors (%d) > nMol (%d)\n"
285 >              "\tThe number of processors is larger than\n"
286 >              "\tthe number of molecules.  This will not result in a \n"
287 >              "\tusable division of atoms for force decomposition.\n"
288 >              "\tEither try a smaller number of processors, or run the\n"
289 >              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
290 >      
291 >      painCave.isFatal = 1;
292 >      simError();
293      }
294 <
294 >    
295 >    int seedValue;
296 >    Globals * simParams = info->getSimParams();
297 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
298 >    if (simParams->haveSeed()) {
299 >      seedValue = simParams->getSeed();
300 >      myRandom = new SeqRandNumGen(seedValue);
301 >    }else {
302 >      myRandom = new SeqRandNumGen();
303 >    }  
304 >    
305 >    
306      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
307 <
307 >    
308      //initialize atomsPerProc
309      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
310 <
310 >    
311      if (worldRank == 0) {
312 <        numerator = info->getNGlobalAtoms();
313 <        denominator = nProcessors;
314 <        precast = numerator / denominator;
315 <        nTarget = (int)(precast + 0.5);
316 <
317 <        for(i = 0; i < nGlobalMols; i++) {
318 <            done = 0;
319 <            loops = 0;
320 <
321 <            while (!done) {
322 <                loops++;
323 <
324 <                // Pick a processor at random
325 <
326 <                which_proc = (int) (myRandom.getRandom() * nProcessors);
327 <
328 <                //get the molecule stamp first
329 <                int stampId = info->getMoleculeStampId(i);
330 <                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
331 <
332 <                // How many atoms does this processor have so far?
333 <                old_atoms = atomsPerProc[which_proc];
334 <                add_atoms = moleculeStamp->getNAtoms();
335 <                new_atoms = old_atoms + add_atoms;
336 <
337 <                // If we've been through this loop too many times, we need
338 <                // to just give up and assign the molecule to this processor
339 <                // and be done with it.
340 <
341 <                if (loops > 100) {
342 <                    sprintf(painCave.errMsg,
343 <                            "I've tried 100 times to assign molecule %d to a "
344 <                                " processor, but can't find a good spot.\n"
345 <                                "I'm assigning it at random to processor %d.\n",
346 <                            i, which_proc);
347 <
348 <                    painCave.isFatal = 0;
349 <                    simError();
350 <
351 <                    molToProcMap[i] = which_proc;
352 <                    atomsPerProc[which_proc] += add_atoms;
353 <
354 <                    done = 1;
355 <                    continue;
356 <                }
357 <
358 <                // If we can add this molecule to this processor without sending
359 <                // it above nTarget, then go ahead and do it:
360 <
361 <                if (new_atoms <= nTarget) {
362 <                    molToProcMap[i] = which_proc;
363 <                    atomsPerProc[which_proc] += add_atoms;
364 <
365 <                    done = 1;
366 <                    continue;
367 <                }
368 <
369 <                // The only situation left is when new_atoms > nTarget.  We
370 <                // want to accept this with some probability that dies off the
371 <                // farther we are from nTarget
372 <
373 <                // roughly:  x = new_atoms - nTarget
374 <                //           Pacc(x) = exp(- a * x)
375 <                // where a = penalty / (average atoms per molecule)
376 <
377 <                x = (double)(new_atoms - nTarget);
378 <                y = myRandom.getRandom();
379 <
380 <                if (y < exp(- a * x)) {
381 <                    molToProcMap[i] = which_proc;
382 <                    atomsPerProc[which_proc] += add_atoms;
383 <
384 <                    done = 1;
385 <                    continue;
386 <                } else {
387 <                    continue;
388 <                }
392 <            }
312 >      numerator = info->getNGlobalAtoms();
313 >      denominator = nProcessors;
314 >      precast = numerator / denominator;
315 >      nTarget = (int)(precast + 0.5);
316 >      
317 >      for(i = 0; i < nGlobalMols; i++) {
318 >        done = 0;
319 >        loops = 0;
320 >        
321 >        while (!done) {
322 >          loops++;
323 >          
324 >          // Pick a processor at random
325 >          
326 >          which_proc = (int) (myRandom->rand() * nProcessors);
327 >          
328 >          //get the molecule stamp first
329 >          int stampId = info->getMoleculeStampId(i);
330 >          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
331 >          
332 >          // How many atoms does this processor have so far?
333 >          old_atoms = atomsPerProc[which_proc];
334 >          add_atoms = moleculeStamp->getNAtoms();
335 >          new_atoms = old_atoms + add_atoms;
336 >          
337 >          // If we've been through this loop too many times, we need
338 >          // to just give up and assign the molecule to this processor
339 >          // and be done with it.
340 >          
341 >          if (loops > 100) {
342 >            sprintf(painCave.errMsg,
343 >                    "I've tried 100 times to assign molecule %d to a "
344 >                    " processor, but can't find a good spot.\n"
345 >                    "I'm assigning it at random to processor %d.\n",
346 >                    i, which_proc);
347 >            
348 >            painCave.isFatal = 0;
349 >            simError();
350 >            
351 >            molToProcMap[i] = which_proc;
352 >            atomsPerProc[which_proc] += add_atoms;
353 >            
354 >            done = 1;
355 >            continue;
356 >          }
357 >          
358 >          // If we can add this molecule to this processor without sending
359 >          // it above nTarget, then go ahead and do it:
360 >          
361 >          if (new_atoms <= nTarget) {
362 >            molToProcMap[i] = which_proc;
363 >            atomsPerProc[which_proc] += add_atoms;
364 >            
365 >            done = 1;
366 >            continue;
367 >          }
368 >          
369 >          // The only situation left is when new_atoms > nTarget.  We
370 >          // want to accept this with some probability that dies off the
371 >          // farther we are from nTarget
372 >          
373 >          // roughly:  x = new_atoms - nTarget
374 >          //           Pacc(x) = exp(- a * x)
375 >          // where a = penalty / (average atoms per molecule)
376 >          
377 >          x = (double)(new_atoms - nTarget);
378 >          y = myRandom->rand();
379 >          
380 >          if (y < exp(- a * x)) {
381 >            molToProcMap[i] = which_proc;
382 >            atomsPerProc[which_proc] += add_atoms;
383 >            
384 >            done = 1;
385 >            continue;
386 >          } else {
387 >            continue;
388 >          }
389          }
390 <
391 <        // Spray out this nonsense to all other processors:
392 <
393 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
390 >      }
391 >      
392 >      delete myRandom;
393 >      
394 >      // Spray out this nonsense to all other processors:
395 >      
396 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
397      } else {
398 <
399 <        // Listen to your marching orders from processor 0:
400 <
401 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
398 >      
399 >      // Listen to your marching orders from processor 0:
400 >      
401 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
402      }
403 <
403 >    
404      info->setMolToProcMap(molToProcMap);
405      sprintf(checkPointMsg,
406              "Successfully divided the molecules among the processors.\n");
407      MPIcheckPoint();
408 < }
409 <
408 >  }
409 >  
410   #endif
411 <
412 < void SimCreator::createMolecules(SimInfo *info) {
411 >  
412 >  void SimCreator::createMolecules(SimInfo *info) {
413      MoleculeCreator molCreator;
414      int stampId;
415 <
415 >    
416      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
417 <
417 >      
418   #ifdef IS_MPI
419 <
420 <        if (info->getMolToProc(i) == worldRank) {
419 >      
420 >      if (info->getMolToProc(i) == worldRank) {
421   #endif
422 <
423 <            stampId = info->getMoleculeStampId(i);
424 <            Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
425 <                                                                                    stampId, i, info->getLocalIndexManager());
426 <
427 <            info->addMolecule(mol);
428 <
422 >        
423 >        stampId = info->getMoleculeStampId(i);
424 >        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
425 >                                                   stampId, i, info->getLocalIndexManager());
426 >        
427 >        info->addMolecule(mol);
428 >        
429   #ifdef IS_MPI
430 <
431 <        }
432 <
430 >        
431 >      }
432 >      
433   #endif
434 <
434 >      
435      } //end for(int i=0)  
436 < }
437 <
438 < void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
440 <                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
441 <    int i;
442 <    char * id;
443 <    MoleculeStamp * currentStamp;
444 <    Component** the_components = simParams->getComponents();
445 <    int n_components = simParams->getNComponents();
446 <
447 <    if (!simParams->haveNMol()) {
448 <        // we don't have the total number of molecules, so we assume it is
449 <        // given in each component
450 <
451 <        for(i = 0; i < n_components; i++) {
452 <            if (!the_components[i]->haveNMol()) {
453 <                // we have a problem
454 <                sprintf(painCave.errMsg,
455 <                        "SimCreator Error. No global NMol or component NMol given.\n"
456 <                            "\tCannot calculate the number of atoms.\n");
457 <
458 <                painCave.isFatal = 1;
459 <                simError();
460 <            }
461 <
462 <            id = the_components[i]->getType();
463 <            currentStamp = (stamps->extractMolStamp(id))->getStamp();
464 <
465 <            if (currentStamp == NULL) {
466 <                sprintf(painCave.errMsg,
467 <                        "SimCreator error: Component \"%s\" was not found in the "
468 <                            "list of declared molecules\n", id);
469 <
470 <                painCave.isFatal = 1;
471 <                simError();
472 <            }
473 <
474 <            moleculeStampPairs.push_back(
475 <                std::make_pair(currentStamp, the_components[i]->getNMol()));
476 <        } //end for (i = 0; i < n_components; i++)
477 <    } else {
478 <        sprintf(painCave.errMsg, "SimSetup error.\n"
479 <                                     "\tSorry, the ability to specify total"
480 <                                     " nMols and then give molfractions in the components\n"
481 <                                     "\tis not currently supported."
482 <                                     " Please give nMol in the components.\n");
483 <
484 <        painCave.isFatal = 1;
485 <        simError();
486 <    }
487 <
488 < #ifdef IS_MPI
489 <
490 <    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
491 <    MPIcheckPoint();
492 <
493 < #endif // is_mpi
494 <
495 < }
496 <
497 < void SimCreator::setGlobalIndex(SimInfo *info) {
436 >  }
437 >    
438 >  void SimCreator::setGlobalIndex(SimInfo *info) {
439      SimInfo::MoleculeIterator mi;
440      Molecule::AtomIterator ai;
441      Molecule::RigidBodyIterator ri;
# Line 509 | Line 450 | void SimCreator::setGlobalIndex(SimInfo *info) {
450      int nGlobalAtoms = info->getNGlobalAtoms();
451      
452   #ifndef IS_MPI
453 <
453 >    
454      beginAtomIndex = 0;
455      beginRigidBodyIndex = 0;
456      beginCutoffGroupIndex = 0;
457 <
457 >    
458   #else
459 <
459 >    
460      int nproc;
461      int myNode;
462 <
462 >    
463      myNode = worldRank;
464      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
465 <
465 >    
466      std::vector < int > tmpAtomsInProc(nproc, 0);
467      std::vector < int > tmpRigidBodiesInProc(nproc, 0);
468      std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
469      std::vector < int > NumAtomsInProc(nproc, 0);
470      std::vector < int > NumRigidBodiesInProc(nproc, 0);
471      std::vector < int > NumCutoffGroupsInProc(nproc, 0);
472 <
472 >    
473      tmpAtomsInProc[myNode] = info->getNAtoms();
474      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
475      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
476 <
476 >    
477      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
478      MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
479                    MPI_SUM, MPI_COMM_WORLD);
# Line 540 | Line 481 | void SimCreator::setGlobalIndex(SimInfo *info) {
481                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
482      MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
483                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
484 <
484 >    
485      beginAtomIndex = 0;
486      beginRigidBodyIndex = 0;
487      beginCutoffGroupIndex = 0;
488 <
488 >    
489      for(int i = 0; i < myNode; i++) {
490 <        beginAtomIndex += NumAtomsInProc[i];
491 <        beginRigidBodyIndex += NumRigidBodiesInProc[i];
492 <        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
490 >      beginAtomIndex += NumAtomsInProc[i];
491 >      beginRigidBodyIndex += NumRigidBodiesInProc[i];
492 >      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
493      }
494 <
494 >    
495   #endif
496 <
496 >    
497      //rigidbody's index begins right after atom's
498      beginRigidBodyIndex += info->getNGlobalAtoms();
499 <
499 >    
500      for(mol = info->beginMolecule(mi); mol != NULL;
501          mol = info->nextMolecule(mi)) {
502 <
503 <        //local index(index in DataStorge) of atom is important
504 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
505 <            atom->setGlobalIndex(beginAtomIndex++);
506 <        }
507 <
508 <        for(rb = mol->beginRigidBody(ri); rb != NULL;
509 <            rb = mol->nextRigidBody(ri)) {
510 <            rb->setGlobalIndex(beginRigidBodyIndex++);
511 <        }
512 <
513 <        //local index of cutoff group is trivial, it only depends on the order of travesing
514 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
515 <            cg = mol->nextCutoffGroup(ci)) {
516 <            cg->setGlobalIndex(beginCutoffGroupIndex++);
517 <        }
502 >      
503 >      //local index(index in DataStorge) of atom is important
504 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
505 >        atom->setGlobalIndex(beginAtomIndex++);
506 >      }
507 >      
508 >      for(rb = mol->beginRigidBody(ri); rb != NULL;
509 >          rb = mol->nextRigidBody(ri)) {
510 >        rb->setGlobalIndex(beginRigidBodyIndex++);
511 >      }
512 >      
513 >      //local index of cutoff group is trivial, it only depends on the order of travesing
514 >      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
515 >          cg = mol->nextCutoffGroup(ci)) {
516 >        cg->setGlobalIndex(beginCutoffGroupIndex++);
517 >      }
518      }
519 <
519 >    
520      //fill globalGroupMembership
521      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
522      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
523 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
524 <
525 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
526 <                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
527 <            }
528 <
529 <        }      
523 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
524 >        
525 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
526 >          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
527 >        }
528 >        
529 >      }      
530      }
531 <
531 >    
532   #ifdef IS_MPI    
533      // Since the globalGroupMembership has been zero filled and we've only
534      // poked values into the atoms we know, we can do an Allreduce
# Line 597 | Line 538 | void SimCreator::setGlobalIndex(SimInfo *info) {
538      std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
539      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
540                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
541 <     info->setGlobalGroupMembership(tmpGroupMembership);
541 >    info->setGlobalGroupMembership(tmpGroupMembership);
542   #else
543      info->setGlobalGroupMembership(globalGroupMembership);
544   #endif
545 <
545 >    
546      //fill molMembership
547      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
548      
549      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
550 <
551 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
552 <            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
553 <        }
550 >      
551 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
552 >        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
553 >      }
554      }
555 <
555 >    
556   #ifdef IS_MPI
557      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
558 <
558 >    
559      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
560                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
561      
# Line 622 | Line 563 | void SimCreator::setGlobalIndex(SimInfo *info) {
563   #else
564      info->setGlobalMolMembership(globalMolMembership);
565   #endif
566 <
567 < }
568 <
569 < void SimCreator::loadCoordinates(SimInfo* info) {
566 >    
567 >  }
568 >  
569 >  void SimCreator::loadCoordinates(SimInfo* info) {
570      Globals* simParams;
571      simParams = info->getSimParams();
572      
573      if (!simParams->haveInitialConfig()) {
574 <        sprintf(painCave.errMsg,
575 <                "Cannot intialize a simulation without an initial configuration file.\n");
576 <        painCave.isFatal = 1;;
577 <        simError();
574 >      sprintf(painCave.errMsg,
575 >              "Cannot intialize a simulation without an initial configuration file.\n");
576 >      painCave.isFatal = 1;;
577 >      simError();
578      }
579 <        
579 >    
580      DumpReader reader(info, simParams->getInitialConfig());
581      int nframes = reader.getNFrames();
582 <
582 >    
583      if (nframes > 0) {
584 <        reader.readFrame(nframes - 1);
584 >      reader.readFrame(nframes - 1);
585      } else {
586 <        //invalid initial coordinate file
587 <        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
588 <                simParams->getInitialConfig());
589 <        painCave.isFatal = 1;
590 <        simError();
586 >      //invalid initial coordinate file
587 >      sprintf(painCave.errMsg,
588 >              "Initial configuration file %s should at least contain one frame\n",
589 >              simParams->getInitialConfig().c_str());
590 >      painCave.isFatal = 1;
591 >      simError();
592      }
593 <
593 >    
594      //copy the current snapshot to previous snapshot
595      info->getSnapshotManager()->advance();
596 < }
597 <
596 >  }
597 >  
598   } //end namespace oopse
599  
600  

Diff Legend

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