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

Comparing trunk/src/brains/SimCreator.cpp (file contents):
Revision 297 by tim, Mon Feb 7 19:14:26 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 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 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  
50 – #include <sprng.h>
51 –
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 + #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
79 < #include "io/mpiBASS.h"
62 < #include "math/randomSPRNG.hpp"
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
151 >      
152 >      catch(antlr::MismatchedCharException& e) {
153 >          sprintf(painCave.errMsg,
154 >                  "parser exception: %s %s:%d:%d\n",
155 >                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
156 >          painCave.isFatal = 1;
157 >          simError();          
158 >      }
159 >      catch(antlr::MismatchedTokenException &e) {
160 >          sprintf(painCave.errMsg,
161 >                  "parser exception: %s %s:%d:%d\n",
162 >                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
163 >          painCave.isFatal = 1;
164 >          simError();  
165 >      }
166 >      catch(antlr::NoViableAltForCharException &e) {
167 >          sprintf(painCave.errMsg,
168 >                  "parser exception: %s %s:%d:%d\n",
169 >                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
170 >          painCave.isFatal = 1;
171 >          simError();  
172 >      }
173 >      catch(antlr::NoViableAltException &e) {
174 >          sprintf(painCave.errMsg,
175 >                  "parser exception: %s %s:%d:%d\n",
176 >                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
177 >          painCave.isFatal = 1;
178 >          simError();  
179 >      }
180 >      
181 >        catch(antlr::TokenStreamRecognitionException& e) {
182 >          sprintf(painCave.errMsg,
183 >                  "parser exception: %s %s:%d:%d\n",
184 >                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
185 >          painCave.isFatal = 1;
186 >          simError();  
187 >        }
188 >        
189 >        catch(antlr::TokenStreamIOException& e) {
190 >          sprintf(painCave.errMsg,
191 >                  "parser exception: %s\n",
192 >                  e.getMessage().c_str());
193 >          painCave.isFatal = 1;
194 >          simError();
195 >        }
196 >        
197 >        catch(antlr::TokenStreamException& e) {
198 >          sprintf(painCave.errMsg,
199 >                  "parser exception: %s\n",
200 >                  e.getMessage().c_str());
201 >          painCave.isFatal = 1;
202 >          simError();
203 >        }        
204 >       catch (antlr::RecognitionException& e) {
205 >          sprintf(painCave.errMsg,
206 >                  "parser exception: %s %s:%d:%d\n",
207 >                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
208 >          painCave.isFatal = 1;
209 >          simError();          
210 >       }
211 >       catch (antlr::CharStreamException& e) {
212 >          sprintf(painCave.errMsg,
213 >                  "parser exception: %s\n",
214 >                  e.getMessage().c_str());
215 >          painCave.isFatal = 1;
216 >          simError();        
217 >       }
218 >       catch (OOPSEException& e) {
219 >          sprintf(painCave.errMsg,
220 >                  "%s\n",
221 >                  e.getMessage().c_str());
222 >          painCave.isFatal = 1;
223 >          simError();
224 >       }
225 >       catch (std::exception& e) {
226 >          sprintf(painCave.errMsg,
227 >                  "parser exception: %s\n",
228 >                  e.what());
229 >          painCave.isFatal = 1;
230 >          simError();
231 >       }
232  
233 <        yacc_BASS(mdFileName.c_str());
233 >        return simParams;
234 >  }
235 >  
236 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
237 >                                  bool loadInitCoords) {
238  
85 – #ifdef IS_MPI
86 –
87 –        throwMPIEvent(NULL);
88 –    } else {
89 –        set_interface_stamps(stamps, simParams);
90 –        mpiEventInit();
91 –        MPIcheckPoint();
92 –        mpiEventLoop();
93 –    }
94 –
95 – #endif
96 –
97 – }
98 –
99 – SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
100 –    
101 –    MakeStamps * stamps = new MakeStamps();
102 –
103 –    Globals * simParams = new Globals();
104 –
239      //parse meta-data file
240 <    parseFile(mdFileName, stamps, simParams);
241 <
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());
249 <        painCave.isFatal = 1;
250 <        simError();
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());
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 <        
265 <        std::string variant = simParams->getForceFieldVariant();
266 <
267 <        std::string::size_type pos = forcefieldFileName.rfind(".frc");
268 <        variant = "." + variant;
269 <        if (pos != std::string::npos) {
270 <            forcefieldFileName.insert(pos, variant);
271 <        } else {
272 <            //If the default force field file name does not containt .frc suffix, just append the .variant
273 <            forcefieldFileName.append(variant);
274 <        }
262 >      //If the force field has variant, the variant force field name will be
263 >      //Base.variant.frc. For exampel EAM.u6.frc
264 >      
265 >      std::string variant = simParams->getForceFieldVariant();
266 >      
267 >      std::string::size_type pos = forcefieldFileName.rfind(".frc");
268 >      variant = "." + variant;
269 >      if (pos != std::string::npos) {
270 >        forcefieldFileName.insert(pos, variant);
271 >      } else {
272 >        //If the default force field file name does not containt .frc suffix, just append the .variant
273 >        forcefieldFileName.append(variant);
274 >      }
275      }
276      
277      ff->parse(forcefieldFileName);
278 <    
144 <    //extract the molecule stamps
145 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
146 <    compList(stamps, simParams, moleculeStampPairs);
147 <
278 >    ff->setFortranForceOptions();
279      //create SimInfo
280 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
281 <
282 <    //gather parameters (SimCreator only retrieves part of the parameters)
280 >    SimInfo * info = new SimInfo(ff, simParams);
281 >    
282 >    //gather parameters (SimCreator only retrieves part of the
283 >    //parameters)
284      gatherParameters(info, mdFileName);
285 <
285 >    
286      //divide the molecules and determine the global index of molecules
287   #ifdef IS_MPI
288      divideMolecules(info);
289   #endif
290 <
290 >    
291      //create the molecules
292      createMolecules(info);
293 <
294 <
295 <    //allocate memory for DataStorage(circular reference, need to break it)
293 >    
294 >    
295 >    //allocate memory for DataStorage(circular reference, need to
296 >    //break it)
297      info->setSnapshotManager(new SimSnapshotManager(info));
298      
299 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
300 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
301 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
299 >    //set the global index of atoms, rigidbodies and cutoffgroups
300 >    //(only need to be set once, the global index will never change
301 >    //again). Local indices of atoms and rigidbodies are already set
302 >    //by MoleculeCreator class which actually delegates the
303 >    //responsibility to LocalIndexManager.
304      setGlobalIndex(info);
305 <
306 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
307 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
308 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
309 <    //we can determine the beginning global indices of atoms before they get created.
305 >    
306 >    //Although addExcludePairs is called inside SimInfo's addMolecule
307 >    //method, at that point atoms don't have the global index yet
308 >    //(their global index are all initialized to -1).  Therefore we
309 >    //have to call addExcludePairs explicitly here. A way to work
310 >    //around is that we can determine the beginning global indices of
311 >    //atoms before they get created.
312      SimInfo::MoleculeIterator mi;
313      Molecule* mol;
314      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
315 <        info->addExcludePairs(mol);
315 >      info->addExcludePairs(mol);
316      }
317      
181 –
182 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
183 –    //eta, chi for NPT integrator)
318      if (loadInitCoords)
319 <        loadCoordinates(info);    
319 >      loadCoordinates(info);    
320      
321      return info;
322 < }
323 <
324 < void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
325 <
326 <    //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
322 >  }
323 >  
324 >  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
325 >    
326 >    //figure out the output file names
327      std::string prefix;
328 <
328 >    
329   #ifdef IS_MPI
330 <
330 >    
331      if (worldRank == 0) {
332   #endif // is_mpi
333 <
334 <        if (simParams->haveFinalConfig()) {
335 <            prefix = getPrefix(simParams->getFinalConfig());
336 <        } else {
337 <            prefix = getPrefix(mdfile);
338 <        }
339 <
340 <        info->setFinalConfigFileName(prefix + ".eor");
341 <        info->setDumpFileName(prefix + ".dump");
342 <        info->setStatFileName(prefix + ".stat");
343 <
333 >      Globals * simParams = info->getSimParams();
334 >      if (simParams->haveFinalConfig()) {
335 >        prefix = getPrefix(simParams->getFinalConfig());
336 >      } else {
337 >        prefix = getPrefix(mdfile);
338 >      }
339 >      
340 >      info->setFinalConfigFileName(prefix + ".eor");
341 >      info->setDumpFileName(prefix + ".dump");
342 >      info->setStatFileName(prefix + ".stat");
343 >      info->setRestFileName(prefix + ".zang");
344 >      
345   #ifdef IS_MPI
346 <
346 >      
347      }
348 <
348 >    
349   #endif
350 <
351 < }
352 <
350 >    
351 >  }
352 >  
353   #ifdef IS_MPI
354 < void SimCreator::divideMolecules(SimInfo *info) {
355 <    double numerator;
356 <    double denominator;
357 <    double precast;
358 <    double x;
359 <    double y;
360 <    double a;
354 >  void SimCreator::divideMolecules(SimInfo *info) {
355 >    RealType numerator;
356 >    RealType denominator;
357 >    RealType precast;
358 >    RealType x;
359 >    RealType y;
360 >    RealType a;
361      int old_atoms;
362      int add_atoms;
363      int new_atoms;
# Line 287 | Line 369 | void SimCreator::divideMolecules(SimInfo *info) {
369      int which_proc;
370      int nProcessors;
371      std::vector<int> atomsPerProc;
290 –    randomSPRNG myRandom(info->getSeed());
372      int nGlobalMols = info->getNGlobalMolecules();
373      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
374      
375      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
376 <
376 >    
377      if (nProcessors > nGlobalMols) {
378 <        sprintf(painCave.errMsg,
379 <                "nProcessors (%d) > nMol (%d)\n"
380 <                    "\tThe number of processors is larger than\n"
381 <                    "\tthe number of molecules.  This will not result in a \n"
382 <                    "\tusable division of atoms for force decomposition.\n"
383 <                    "\tEither try a smaller number of processors, or run the\n"
384 <                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
385 <
386 <        painCave.isFatal = 1;
387 <        simError();
378 >      sprintf(painCave.errMsg,
379 >              "nProcessors (%d) > nMol (%d)\n"
380 >              "\tThe number of processors is larger than\n"
381 >              "\tthe number of molecules.  This will not result in a \n"
382 >              "\tusable division of atoms for force decomposition.\n"
383 >              "\tEither try a smaller number of processors, or run the\n"
384 >              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
385 >      
386 >      painCave.isFatal = 1;
387 >      simError();
388      }
389 <
389 >    
390 >    int seedValue;
391 >    Globals * simParams = info->getSimParams();
392 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
393 >    if (simParams->haveSeed()) {
394 >      seedValue = simParams->getSeed();
395 >      myRandom = new SeqRandNumGen(seedValue);
396 >    }else {
397 >      myRandom = new SeqRandNumGen();
398 >    }  
399 >    
400 >    
401      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
402 <
402 >    
403      //initialize atomsPerProc
404      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
405 <
405 >    
406      if (worldRank == 0) {
407 <        numerator = info->getNGlobalAtoms();
408 <        denominator = nProcessors;
409 <        precast = numerator / denominator;
410 <        nTarget = (int)(precast + 0.5);
411 <
412 <        for(i = 0; i < nGlobalMols; i++) {
413 <            done = 0;
414 <            loops = 0;
415 <
416 <            while (!done) {
417 <                loops++;
418 <
419 <                // Pick a processor at random
420 <
421 <                which_proc = (int) (myRandom.getRandom() * nProcessors);
422 <
423 <                //get the molecule stamp first
424 <                int stampId = info->getMoleculeStampId(i);
425 <                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
426 <
427 <                // How many atoms does this processor have so far?
428 <                old_atoms = atomsPerProc[which_proc];
429 <                add_atoms = moleculeStamp->getNAtoms();
430 <                new_atoms = old_atoms + add_atoms;
431 <
432 <                // If we've been through this loop too many times, we need
433 <                // to just give up and assign the molecule to this processor
434 <                // and be done with it.
435 <
436 <                if (loops > 100) {
437 <                    sprintf(painCave.errMsg,
438 <                            "I've tried 100 times to assign molecule %d to a "
439 <                                " processor, but can't find a good spot.\n"
440 <                                "I'm assigning it at random to processor %d.\n",
441 <                            i, which_proc);
442 <
443 <                    painCave.isFatal = 0;
444 <                    simError();
445 <
446 <                    molToProcMap[i] = which_proc;
447 <                    atomsPerProc[which_proc] += add_atoms;
448 <
449 <                    done = 1;
450 <                    continue;
451 <                }
452 <
453 <                // If we can add this molecule to this processor without sending
454 <                // it above nTarget, then go ahead and do it:
455 <
456 <                if (new_atoms <= nTarget) {
457 <                    molToProcMap[i] = which_proc;
458 <                    atomsPerProc[which_proc] += add_atoms;
459 <
460 <                    done = 1;
461 <                    continue;
462 <                }
463 <
464 <                // The only situation left is when new_atoms > nTarget.  We
465 <                // want to accept this with some probability that dies off the
466 <                // farther we are from nTarget
467 <
468 <                // roughly:  x = new_atoms - nTarget
469 <                //           Pacc(x) = exp(- a * x)
470 <                // where a = penalty / (average atoms per molecule)
471 <
472 <                x = (double)(new_atoms - nTarget);
473 <                y = myRandom.getRandom();
474 <
475 <                if (y < exp(- a * x)) {
476 <                    molToProcMap[i] = which_proc;
477 <                    atomsPerProc[which_proc] += add_atoms;
478 <
479 <                    done = 1;
480 <                    continue;
481 <                } else {
482 <                    continue;
483 <                }
392 <            }
407 >      numerator = info->getNGlobalAtoms();
408 >      denominator = nProcessors;
409 >      precast = numerator / denominator;
410 >      nTarget = (int)(precast + 0.5);
411 >      
412 >      for(i = 0; i < nGlobalMols; i++) {
413 >        done = 0;
414 >        loops = 0;
415 >        
416 >        while (!done) {
417 >          loops++;
418 >          
419 >          // Pick a processor at random
420 >          
421 >          which_proc = (int) (myRandom->rand() * nProcessors);
422 >          
423 >          //get the molecule stamp first
424 >          int stampId = info->getMoleculeStampId(i);
425 >          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
426 >          
427 >          // How many atoms does this processor have so far?
428 >          old_atoms = atomsPerProc[which_proc];
429 >          add_atoms = moleculeStamp->getNAtoms();
430 >          new_atoms = old_atoms + add_atoms;
431 >          
432 >          // If we've been through this loop too many times, we need
433 >          // to just give up and assign the molecule to this processor
434 >          // and be done with it.
435 >          
436 >          if (loops > 100) {
437 >            sprintf(painCave.errMsg,
438 >                    "I've tried 100 times to assign molecule %d to a "
439 >                    " processor, but can't find a good spot.\n"
440 >                    "I'm assigning it at random to processor %d.\n",
441 >                    i, which_proc);
442 >            
443 >            painCave.isFatal = 0;
444 >            simError();
445 >            
446 >            molToProcMap[i] = which_proc;
447 >            atomsPerProc[which_proc] += add_atoms;
448 >            
449 >            done = 1;
450 >            continue;
451 >          }
452 >          
453 >          // If we can add this molecule to this processor without sending
454 >          // it above nTarget, then go ahead and do it:
455 >          
456 >          if (new_atoms <= nTarget) {
457 >            molToProcMap[i] = which_proc;
458 >            atomsPerProc[which_proc] += add_atoms;
459 >            
460 >            done = 1;
461 >            continue;
462 >          }
463 >          
464 >          // The only situation left is when new_atoms > nTarget.  We
465 >          // want to accept this with some probability that dies off the
466 >          // farther we are from nTarget
467 >          
468 >          // roughly:  x = new_atoms - nTarget
469 >          //           Pacc(x) = exp(- a * x)
470 >          // where a = penalty / (average atoms per molecule)
471 >          
472 >          x = (RealType)(new_atoms - nTarget);
473 >          y = myRandom->rand();
474 >          
475 >          if (y < exp(- a * x)) {
476 >            molToProcMap[i] = which_proc;
477 >            atomsPerProc[which_proc] += add_atoms;
478 >            
479 >            done = 1;
480 >            continue;
481 >          } else {
482 >            continue;
483 >          }
484          }
485 <
486 <        // Spray out this nonsense to all other processors:
487 <
488 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
485 >      }
486 >      
487 >      delete myRandom;
488 >      
489 >      // Spray out this nonsense to all other processors:
490 >      
491 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
492      } else {
493 <
494 <        // Listen to your marching orders from processor 0:
495 <
496 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
493 >      
494 >      // Listen to your marching orders from processor 0:
495 >      
496 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
497      }
498 <
498 >    
499      info->setMolToProcMap(molToProcMap);
500      sprintf(checkPointMsg,
501              "Successfully divided the molecules among the processors.\n");
502      MPIcheckPoint();
503 < }
504 <
503 >  }
504 >  
505   #endif
506 <
507 < void SimCreator::createMolecules(SimInfo *info) {
506 >  
507 >  void SimCreator::createMolecules(SimInfo *info) {
508      MoleculeCreator molCreator;
509      int stampId;
510 <
510 >    
511      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
512 <
512 >      
513   #ifdef IS_MPI
514 <
515 <        if (info->getMolToProc(i) == worldRank) {
514 >      
515 >      if (info->getMolToProc(i) == worldRank) {
516   #endif
517 <
518 <            stampId = info->getMoleculeStampId(i);
519 <            Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
520 <                                                                                    stampId, i, info->getLocalIndexManager());
521 <
522 <            info->addMolecule(mol);
523 <
517 >        
518 >        stampId = info->getMoleculeStampId(i);
519 >        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
520 >                                                   stampId, i, info->getLocalIndexManager());
521 >        
522 >        info->addMolecule(mol);
523 >        
524   #ifdef IS_MPI
525 <
526 <        }
527 <
525 >        
526 >      }
527 >      
528   #endif
529 <
529 >      
530      } //end for(int i=0)  
531 < }
532 <
533 < 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) {
531 >  }
532 >    
533 >  void SimCreator::setGlobalIndex(SimInfo *info) {
534      SimInfo::MoleculeIterator mi;
535      Molecule::AtomIterator ai;
536      Molecule::RigidBodyIterator ri;
# Line 509 | Line 545 | void SimCreator::setGlobalIndex(SimInfo *info) {
545      int nGlobalAtoms = info->getNGlobalAtoms();
546      
547   #ifndef IS_MPI
548 <
548 >    
549      beginAtomIndex = 0;
550      beginRigidBodyIndex = 0;
551      beginCutoffGroupIndex = 0;
552 <
552 >    
553   #else
554 <
554 >    
555      int nproc;
556      int myNode;
557 <
557 >    
558      myNode = worldRank;
559      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
560 <
560 >    
561      std::vector < int > tmpAtomsInProc(nproc, 0);
562      std::vector < int > tmpRigidBodiesInProc(nproc, 0);
563      std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
564      std::vector < int > NumAtomsInProc(nproc, 0);
565      std::vector < int > NumRigidBodiesInProc(nproc, 0);
566      std::vector < int > NumCutoffGroupsInProc(nproc, 0);
567 <
567 >    
568      tmpAtomsInProc[myNode] = info->getNAtoms();
569      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
570      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
571 <
571 >    
572      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
573      MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
574                    MPI_SUM, MPI_COMM_WORLD);
# Line 540 | Line 576 | void SimCreator::setGlobalIndex(SimInfo *info) {
576                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
577      MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
578                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
579 <
579 >    
580      beginAtomIndex = 0;
581      beginRigidBodyIndex = 0;
582      beginCutoffGroupIndex = 0;
583 <
583 >    
584      for(int i = 0; i < myNode; i++) {
585 <        beginAtomIndex += NumAtomsInProc[i];
586 <        beginRigidBodyIndex += NumRigidBodiesInProc[i];
587 <        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
585 >      beginAtomIndex += NumAtomsInProc[i];
586 >      beginRigidBodyIndex += NumRigidBodiesInProc[i];
587 >      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
588      }
589 <
589 >    
590   #endif
591 <
591 >    
592      //rigidbody's index begins right after atom's
593      beginRigidBodyIndex += info->getNGlobalAtoms();
594 <
594 >    
595      for(mol = info->beginMolecule(mi); mol != NULL;
596          mol = info->nextMolecule(mi)) {
597 <
598 <        //local index(index in DataStorge) of atom is important
599 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
600 <            atom->setGlobalIndex(beginAtomIndex++);
601 <        }
602 <
603 <        for(rb = mol->beginRigidBody(ri); rb != NULL;
604 <            rb = mol->nextRigidBody(ri)) {
605 <            rb->setGlobalIndex(beginRigidBodyIndex++);
606 <        }
607 <
608 <        //local index of cutoff group is trivial, it only depends on the order of travesing
609 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
610 <            cg = mol->nextCutoffGroup(ci)) {
611 <            cg->setGlobalIndex(beginCutoffGroupIndex++);
612 <        }
597 >      
598 >      //local index(index in DataStorge) of atom is important
599 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
600 >        atom->setGlobalIndex(beginAtomIndex++);
601 >      }
602 >      
603 >      for(rb = mol->beginRigidBody(ri); rb != NULL;
604 >          rb = mol->nextRigidBody(ri)) {
605 >        rb->setGlobalIndex(beginRigidBodyIndex++);
606 >      }
607 >      
608 >      //local index of cutoff group is trivial, it only depends on the order of travesing
609 >      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
610 >          cg = mol->nextCutoffGroup(ci)) {
611 >        cg->setGlobalIndex(beginCutoffGroupIndex++);
612 >      }
613      }
614 <
614 >    
615      //fill globalGroupMembership
616      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
617      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
618 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
619 <
620 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
621 <                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
622 <            }
623 <
624 <        }      
618 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
619 >        
620 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
621 >          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
622 >        }
623 >        
624 >      }      
625      }
626 <
626 >    
627   #ifdef IS_MPI    
628      // Since the globalGroupMembership has been zero filled and we've only
629      // poked values into the atoms we know, we can do an Allreduce
# Line 597 | Line 633 | void SimCreator::setGlobalIndex(SimInfo *info) {
633      std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
634      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
635                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
636 <     info->setGlobalGroupMembership(tmpGroupMembership);
636 >    info->setGlobalGroupMembership(tmpGroupMembership);
637   #else
638      info->setGlobalGroupMembership(globalGroupMembership);
639   #endif
640 <
640 >    
641      //fill molMembership
642      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
643      
644      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
645 <
646 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
647 <            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
648 <        }
645 >      
646 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
647 >        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
648 >      }
649      }
650 <
650 >    
651   #ifdef IS_MPI
652      std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
653 <
653 >    
654      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
655                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
656      
# Line 622 | Line 658 | void SimCreator::setGlobalIndex(SimInfo *info) {
658   #else
659      info->setGlobalMolMembership(globalMolMembership);
660   #endif
661 <
662 < }
663 <
664 < void SimCreator::loadCoordinates(SimInfo* info) {
661 >    
662 >  }
663 >  
664 >  void SimCreator::loadCoordinates(SimInfo* info) {
665      Globals* simParams;
666      simParams = info->getSimParams();
667      
668      if (!simParams->haveInitialConfig()) {
669 <        sprintf(painCave.errMsg,
670 <                "Cannot intialize a simulation without an initial configuration file.\n");
671 <        painCave.isFatal = 1;;
672 <        simError();
669 >      sprintf(painCave.errMsg,
670 >              "Cannot intialize a simulation without an initial configuration file.\n");
671 >      painCave.isFatal = 1;;
672 >      simError();
673      }
674 <        
674 >    
675      DumpReader reader(info, simParams->getInitialConfig());
676      int nframes = reader.getNFrames();
677 <
677 >    
678      if (nframes > 0) {
679 <        reader.readFrame(nframes - 1);
679 >      reader.readFrame(nframes - 1);
680      } else {
681 <        //invalid initial coordinate file
682 <        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
683 <                simParams->getInitialConfig());
684 <        painCave.isFatal = 1;
685 <        simError();
681 >      //invalid initial coordinate file
682 >      sprintf(painCave.errMsg,
683 >              "Initial configuration file %s should at least contain one frame\n",
684 >              simParams->getInitialConfig().c_str());
685 >      painCave.isFatal = 1;
686 >      simError();
687      }
688 <
688 >    
689      //copy the current snapshot to previous snapshot
690      info->getSnapshotManager()->advance();
691 < }
692 <
691 >  }
692 >  
693   } //end namespace oopse
694  
695  

Diff Legend

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