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

Diff Legend

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