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 1390 by gezelter, Wed Nov 25 20:02:06 2009 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 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   /**
# 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 < #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 + #include "math/ParallelRandNumGen.hpp"
80 + #endif
81  
82 <    if (worldRank == 0) {
83 < #endif // is_mpi
82 > namespace OpenMD {
83 >  
84 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){
85 >    Globals* simParams = NULL;
86 >    try {
87  
88 <        simParams->initalize();
89 <        set_interface_stamps(stamps, simParams);
88 >      // Create a preprocessor that preprocesses md file into an ostringstream
89 >      std::stringstream ppStream;
90 > #ifdef IS_MPI            
91 >      int streamSize;
92 >      const int masterNode = 0;
93 >      int commStatus;
94 >      if (worldRank == masterNode) {
95 > #endif
96 >                
97 >        SimplePreprocessor preprocessor;
98 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
99 >                
100 > #ifdef IS_MPI            
101 >        //brocasting the stream size
102 >        streamSize = ppStream.str().size() +1;
103 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
104  
105 < #ifdef IS_MPI
105 >        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
106 >            
107 >                
108 >      } else {
109 >        //get stream size
110 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
111  
112 <        mpiEventInit();
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(filename);
126 >      lexer.initDeferredLineCount();
127 >    
128 >      // Create a parser that reads from the scanner
129 >      MDParser parser(lexer);
130 >      parser.setFilename(filename);
131  
132 <        yacc_BASS(mdFileName.c_str());
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  
150 <        throwMPIEvent(NULL);
151 <    } else {
152 <        set_interface_stamps(stamps, simParams);
153 <        mpiEventInit();
154 <        MPIcheckPoint();
155 <        mpiEventLoop();
150 >      
151 >    catch(antlr::MismatchedCharException& e) {
152 >      sprintf(painCave.errMsg,
153 >              "parser exception: %s %s:%d:%d\n",
154 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
155 >      painCave.isFatal = 1;
156 >      simError();          
157      }
158 +    catch(antlr::MismatchedTokenException &e) {
159 +      sprintf(painCave.errMsg,
160 +              "parser exception: %s %s:%d:%d\n",
161 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
162 +      painCave.isFatal = 1;
163 +      simError();  
164 +    }
165 +    catch(antlr::NoViableAltForCharException &e) {
166 +      sprintf(painCave.errMsg,
167 +              "parser exception: %s %s:%d:%d\n",
168 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
169 +      painCave.isFatal = 1;
170 +      simError();  
171 +    }
172 +    catch(antlr::NoViableAltException &e) {
173 +      sprintf(painCave.errMsg,
174 +              "parser exception: %s %s:%d:%d\n",
175 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
176 +      painCave.isFatal = 1;
177 +      simError();  
178 +    }
179 +      
180 +    catch(antlr::TokenStreamRecognitionException& e) {
181 +      sprintf(painCave.errMsg,
182 +              "parser exception: %s %s:%d:%d\n",
183 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
184 +      painCave.isFatal = 1;
185 +      simError();  
186 +    }
187 +        
188 +    catch(antlr::TokenStreamIOException& e) {
189 +      sprintf(painCave.errMsg,
190 +              "parser exception: %s\n",
191 +              e.getMessage().c_str());
192 +      painCave.isFatal = 1;
193 +      simError();
194 +    }
195 +        
196 +    catch(antlr::TokenStreamException& e) {
197 +      sprintf(painCave.errMsg,
198 +              "parser exception: %s\n",
199 +              e.getMessage().c_str());
200 +      painCave.isFatal = 1;
201 +      simError();
202 +    }        
203 +    catch (antlr::RecognitionException& e) {
204 +      sprintf(painCave.errMsg,
205 +              "parser exception: %s %s:%d:%d\n",
206 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
207 +      painCave.isFatal = 1;
208 +      simError();          
209 +    }
210 +    catch (antlr::CharStreamException& e) {
211 +      sprintf(painCave.errMsg,
212 +              "parser exception: %s\n",
213 +              e.getMessage().c_str());
214 +      painCave.isFatal = 1;
215 +      simError();        
216 +    }
217 +    catch (OpenMDException& e) {
218 +      sprintf(painCave.errMsg,
219 +              "%s\n",
220 +              e.getMessage().c_str());
221 +      painCave.isFatal = 1;
222 +      simError();
223 +    }
224 +    catch (std::exception& e) {
225 +      sprintf(painCave.errMsg,
226 +              "parser exception: %s\n",
227 +              e.what());
228 +      painCave.isFatal = 1;
229 +      simError();
230 +    }
231  
232 < #endif
232 >    return simParams;
233 >  }
234 >  
235 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
236 >                                  bool loadInitCoords) {
237  
238 < }
238 >    const int bufferSize = 65535;
239 >    char buffer[bufferSize];
240 >    int lineNo = 0;
241 >    std::string mdRawData;
242 >    int metaDataBlockStart = -1;
243 >    int metaDataBlockEnd = -1;
244 >    int i;
245 >    int mdOffset;
246  
247 < SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
248 <    
249 <    MakeStamps * stamps = new MakeStamps();
247 > #ifdef IS_MPI            
248 >    const int masterNode = 0;
249 >    if (worldRank == masterNode) {
250 > #endif
251  
252 <    Globals * simParams = new Globals();
252 >      std::ifstream mdFile_(mdFileName.c_str());
253 >      
254 >      if (mdFile_.fail()) {
255 >        sprintf(painCave.errMsg,
256 >                "SimCreator: Cannot open file: %s\n",
257 >                mdFileName.c_str());
258 >        painCave.isFatal = 1;
259 >        simError();
260 >      }
261  
262 <    //parse meta-data file
263 <    parseFile(mdFileName, stamps, simParams);
262 >      mdFile_.getline(buffer, bufferSize);
263 >      ++lineNo;
264 >      std::string line = trimLeftCopy(buffer);
265 >      i = CaseInsensitiveFind(line, "<OpenMD");
266 >      if (static_cast<size_t>(i) == string::npos) {
267 >        // try the older file strings to see if that works:
268 >        i = CaseInsensitiveFind(line, "<OOPSE");
269 >      }
270 >      
271 >      if (static_cast<size_t>(i) == string::npos) {
272 >        // still no luck!
273 >        sprintf(painCave.errMsg,
274 >                "SimCreator: File: %s is not a valid OpenMD file!\n",
275 >                mdFileName.c_str());
276 >        painCave.isFatal = 1;
277 >        simError();
278 >      }
279  
280 <    //create the force field
281 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
282 <                          simParams->getForceField());
283 <    
284 <    if (ff == NULL) {
285 <        sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
286 <                simParams->getForceField());
287 <        painCave.isFatal = 1;
288 <        simError();
280 >      //scan through the input stream and find MetaData tag        
281 >      while(mdFile_.getline(buffer, bufferSize)) {
282 >        ++lineNo;
283 >        
284 >        std::string line = trimLeftCopy(buffer);
285 >        if (metaDataBlockStart == -1) {
286 >          i = CaseInsensitiveFind(line, "<MetaData>");
287 >          if (i != string::npos) {
288 >            metaDataBlockStart = lineNo;
289 >            mdOffset = mdFile_.tellg();
290 >          }
291 >        } else {
292 >          i = CaseInsensitiveFind(line, "</MetaData>");
293 >          if (i != string::npos) {
294 >            metaDataBlockEnd = lineNo;
295 >          }
296 >        }
297 >      }
298 >
299 >      if (metaDataBlockStart == -1) {
300 >        sprintf(painCave.errMsg,
301 >                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
302 >                mdFileName.c_str());
303 >        painCave.isFatal = 1;
304 >        simError();
305 >      }
306 >      if (metaDataBlockEnd == -1) {
307 >        sprintf(painCave.errMsg,
308 >                "SimCreator: File: %s did not contain a closed MetaData block!\n",
309 >                mdFileName.c_str());
310 >        painCave.isFatal = 1;
311 >        simError();
312 >      }
313 >        
314 >      mdFile_.clear();
315 >      mdFile_.seekg(0);
316 >      mdFile_.seekg(mdOffset);
317 >
318 >      mdRawData.clear();
319 >
320 >      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
321 >        mdFile_.getline(buffer, bufferSize);
322 >        mdRawData += buffer;
323 >        mdRawData += "\n";
324 >      }
325 >
326 >      mdFile_.close();
327 >
328 > #ifdef IS_MPI
329      }
330 + #endif
331  
332 +    std::stringstream rawMetaDataStream(mdRawData);
333 +
334 +    //parse meta-data file
335 +    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1);
336 +    
337 +    //create the force field
338 +    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
339 +
340 +    if (ff == NULL) {
341 +      sprintf(painCave.errMsg,
342 +              "ForceField Factory can not create %s force field\n",
343 +              simParams->getForceField().c_str());
344 +      painCave.isFatal = 1;
345 +      simError();
346 +    }
347 +    
348      if (simParams->haveForceFieldFileName()) {
349 <        ff->setForceFieldFileName(simParams->getForceFieldFileName());
349 >      ff->setForceFieldFileName(simParams->getForceFieldFileName());
350      }
351      
352      std::string forcefieldFileName;
353      forcefieldFileName = ff->getForceFieldFileName();
354 <
354 >    
355      if (simParams->haveForceFieldVariant()) {
356 <        //If the force field has variant, the variant force field name will be
357 <        //Base.variant.frc. For exampel EAM.u6.frc
358 <        
359 <        std::string variant = simParams->getForceFieldVariant();
360 <
361 <        std::string::size_type pos = forcefieldFileName.rfind(".frc");
362 <        variant = "." + variant;
363 <        if (pos != std::string::npos) {
364 <            forcefieldFileName.insert(pos, variant);
365 <        } else {
366 <            //If the default force field file name does not containt .frc suffix, just append the .variant
367 <            forcefieldFileName.append(variant);
368 <        }
356 >      //If the force field has variant, the variant force field name will be
357 >      //Base.variant.frc. For exampel EAM.u6.frc
358 >      
359 >      std::string variant = simParams->getForceFieldVariant();
360 >      
361 >      std::string::size_type pos = forcefieldFileName.rfind(".frc");
362 >      variant = "." + variant;
363 >      if (pos != std::string::npos) {
364 >        forcefieldFileName.insert(pos, variant);
365 >      } else {
366 >        //If the default force field file name does not containt .frc suffix, just append the .variant
367 >        forcefieldFileName.append(variant);
368 >      }
369      }
370      
371      ff->parse(forcefieldFileName);
372 <    
144 <    //extract the molecule stamps
145 <    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
146 <    compList(stamps, simParams, moleculeStampPairs);
147 <
372 >    ff->setFortranForceOptions();
373      //create SimInfo
374 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
374 >    SimInfo * info = new SimInfo(ff, simParams);
375  
376 <    //gather parameters (SimCreator only retrieves part of the parameters)
376 >    info->setRawMetaData(mdRawData);
377 >    
378 >    //gather parameters (SimCreator only retrieves part of the
379 >    //parameters)
380      gatherParameters(info, mdFileName);
381 <
381 >    
382      //divide the molecules and determine the global index of molecules
383   #ifdef IS_MPI
384      divideMolecules(info);
385   #endif
386 <
386 >    
387      //create the molecules
388      createMolecules(info);
389 <
390 <
391 <    //allocate memory for DataStorage(circular reference, need to break it)
389 >    
390 >    
391 >    //allocate memory for DataStorage(circular reference, need to
392 >    //break it)
393      info->setSnapshotManager(new SimSnapshotManager(info));
394      
395 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
396 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
397 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
395 >    //set the global index of atoms, rigidbodies and cutoffgroups
396 >    //(only need to be set once, the global index will never change
397 >    //again). Local indices of atoms and rigidbodies are already set
398 >    //by MoleculeCreator class which actually delegates the
399 >    //responsibility to LocalIndexManager.
400      setGlobalIndex(info);
401 <
402 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
403 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
404 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
405 <    //we can determine the beginning global indices of atoms before they get created.
401 >    
402 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
403 >    //method, at that point atoms don't have the global index yet
404 >    //(their global index are all initialized to -1).  Therefore we
405 >    //have to call addInteractionPairs explicitly here. A way to work
406 >    //around is that we can determine the beginning global indices of
407 >    //atoms before they get created.
408      SimInfo::MoleculeIterator mi;
409      Molecule* mol;
410      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
411 <        info->addExcludePairs(mol);
411 >      info->addInteractionPairs(mol);
412      }
413      
181 –
182 –    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
183 –    //eta, chi for NPT integrator)
414      if (loadInitCoords)
415 <        loadCoordinates(info);    
415 >      loadCoordinates(info, mdFileName);    
416      
417      return info;
418 < }
419 <
420 < void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
421 <
422 <    //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
418 >  }
419 >  
420 >  void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
421 >    
422 >    //figure out the output file names
423      std::string prefix;
424 <
424 >    
425   #ifdef IS_MPI
426 <
426 >    
427      if (worldRank == 0) {
428   #endif // is_mpi
429 <
430 <        if (simParams->haveFinalConfig()) {
431 <            prefix = getPrefix(simParams->getFinalConfig());
432 <        } else {
433 <            prefix = getPrefix(mdfile);
434 <        }
435 <
436 <        info->setFinalConfigFileName(prefix + ".eor");
437 <        info->setDumpFileName(prefix + ".dump");
438 <        info->setStatFileName(prefix + ".stat");
439 <
429 >      Globals * simParams = info->getSimParams();
430 >      if (simParams->haveFinalConfig()) {
431 >        prefix = getPrefix(simParams->getFinalConfig());
432 >      } else {
433 >        prefix = getPrefix(mdfile);
434 >      }
435 >      
436 >      info->setFinalConfigFileName(prefix + ".eor");
437 >      info->setDumpFileName(prefix + ".dump");
438 >      info->setStatFileName(prefix + ".stat");
439 >      info->setRestFileName(prefix + ".zang");
440 >      
441   #ifdef IS_MPI
442 <
442 >      
443      }
444 <
444 >    
445   #endif
446 <
447 < }
448 <
446 >    
447 >  }
448 >  
449   #ifdef IS_MPI
450 < void SimCreator::divideMolecules(SimInfo *info) {
451 <    double numerator;
452 <    double denominator;
453 <    double precast;
454 <    double x;
455 <    double y;
456 <    double a;
450 >  void SimCreator::divideMolecules(SimInfo *info) {
451 >    RealType numerator;
452 >    RealType denominator;
453 >    RealType precast;
454 >    RealType x;
455 >    RealType y;
456 >    RealType a;
457      int old_atoms;
458      int add_atoms;
459      int new_atoms;
# Line 287 | Line 465 | void SimCreator::divideMolecules(SimInfo *info) {
465      int which_proc;
466      int nProcessors;
467      std::vector<int> atomsPerProc;
290 –    randomSPRNG myRandom(info->getSeed());
468      int nGlobalMols = info->getNGlobalMolecules();
469      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
470      
471      MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
472 <
472 >    
473      if (nProcessors > nGlobalMols) {
474 <        sprintf(painCave.errMsg,
475 <                "nProcessors (%d) > nMol (%d)\n"
476 <                    "\tThe number of processors is larger than\n"
477 <                    "\tthe number of molecules.  This will not result in a \n"
478 <                    "\tusable division of atoms for force decomposition.\n"
479 <                    "\tEither try a smaller number of processors, or run the\n"
480 <                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
481 <
482 <        painCave.isFatal = 1;
483 <        simError();
474 >      sprintf(painCave.errMsg,
475 >              "nProcessors (%d) > nMol (%d)\n"
476 >              "\tThe number of processors is larger than\n"
477 >              "\tthe number of molecules.  This will not result in a \n"
478 >              "\tusable division of atoms for force decomposition.\n"
479 >              "\tEither try a smaller number of processors, or run the\n"
480 >              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
481 >      
482 >      painCave.isFatal = 1;
483 >      simError();
484      }
485 <
485 >    
486 >    int seedValue;
487 >    Globals * simParams = info->getSimParams();
488 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
489 >    if (simParams->haveSeed()) {
490 >      seedValue = simParams->getSeed();
491 >      myRandom = new SeqRandNumGen(seedValue);
492 >    }else {
493 >      myRandom = new SeqRandNumGen();
494 >    }  
495 >    
496 >    
497      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
498 <
498 >    
499      //initialize atomsPerProc
500      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
501 <
501 >    
502      if (worldRank == 0) {
503 <        numerator = info->getNGlobalAtoms();
504 <        denominator = nProcessors;
505 <        precast = numerator / denominator;
506 <        nTarget = (int)(precast + 0.5);
507 <
508 <        for(i = 0; i < nGlobalMols; i++) {
509 <            done = 0;
510 <            loops = 0;
511 <
512 <            while (!done) {
513 <                loops++;
514 <
515 <                // Pick a processor at random
516 <
517 <                which_proc = (int) (myRandom.getRandom() * nProcessors);
518 <
519 <                //get the molecule stamp first
520 <                int stampId = info->getMoleculeStampId(i);
521 <                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
522 <
523 <                // How many atoms does this processor have so far?
524 <                old_atoms = atomsPerProc[which_proc];
525 <                add_atoms = moleculeStamp->getNAtoms();
526 <                new_atoms = old_atoms + add_atoms;
527 <
528 <                // If we've been through this loop too many times, we need
529 <                // to just give up and assign the molecule to this processor
530 <                // and be done with it.
531 <
532 <                if (loops > 100) {
533 <                    sprintf(painCave.errMsg,
534 <                            "I've tried 100 times to assign molecule %d to a "
535 <                                " processor, but can't find a good spot.\n"
536 <                                "I'm assigning it at random to processor %d.\n",
537 <                            i, which_proc);
538 <
539 <                    painCave.isFatal = 0;
540 <                    simError();
541 <
542 <                    molToProcMap[i] = which_proc;
543 <                    atomsPerProc[which_proc] += add_atoms;
544 <
545 <                    done = 1;
546 <                    continue;
547 <                }
548 <
549 <                // If we can add this molecule to this processor without sending
550 <                // it above nTarget, then go ahead and do it:
551 <
552 <                if (new_atoms <= nTarget) {
553 <                    molToProcMap[i] = which_proc;
554 <                    atomsPerProc[which_proc] += add_atoms;
555 <
556 <                    done = 1;
557 <                    continue;
558 <                }
559 <
560 <                // The only situation left is when new_atoms > nTarget.  We
561 <                // want to accept this with some probability that dies off the
562 <                // farther we are from nTarget
563 <
564 <                // roughly:  x = new_atoms - nTarget
565 <                //           Pacc(x) = exp(- a * x)
566 <                // where a = penalty / (average atoms per molecule)
567 <
568 <                x = (double)(new_atoms - nTarget);
569 <                y = myRandom.getRandom();
570 <
571 <                if (y < exp(- a * x)) {
572 <                    molToProcMap[i] = which_proc;
573 <                    atomsPerProc[which_proc] += add_atoms;
574 <
575 <                    done = 1;
576 <                    continue;
577 <                } else {
578 <                    continue;
579 <                }
392 <            }
503 >      numerator = info->getNGlobalAtoms();
504 >      denominator = nProcessors;
505 >      precast = numerator / denominator;
506 >      nTarget = (int)(precast + 0.5);
507 >      
508 >      for(i = 0; i < nGlobalMols; i++) {
509 >        done = 0;
510 >        loops = 0;
511 >        
512 >        while (!done) {
513 >          loops++;
514 >          
515 >          // Pick a processor at random
516 >          
517 >          which_proc = (int) (myRandom->rand() * nProcessors);
518 >          
519 >          //get the molecule stamp first
520 >          int stampId = info->getMoleculeStampId(i);
521 >          MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
522 >          
523 >          // How many atoms does this processor have so far?
524 >          old_atoms = atomsPerProc[which_proc];
525 >          add_atoms = moleculeStamp->getNAtoms();
526 >          new_atoms = old_atoms + add_atoms;
527 >          
528 >          // If we've been through this loop too many times, we need
529 >          // to just give up and assign the molecule to this processor
530 >          // and be done with it.
531 >          
532 >          if (loops > 100) {
533 >            sprintf(painCave.errMsg,
534 >                    "I've tried 100 times to assign molecule %d to a "
535 >                    " processor, but can't find a good spot.\n"
536 >                    "I'm assigning it at random to processor %d.\n",
537 >                    i, which_proc);
538 >            
539 >            painCave.isFatal = 0;
540 >            simError();
541 >            
542 >            molToProcMap[i] = which_proc;
543 >            atomsPerProc[which_proc] += add_atoms;
544 >            
545 >            done = 1;
546 >            continue;
547 >          }
548 >          
549 >          // If we can add this molecule to this processor without sending
550 >          // it above nTarget, then go ahead and do it:
551 >          
552 >          if (new_atoms <= nTarget) {
553 >            molToProcMap[i] = which_proc;
554 >            atomsPerProc[which_proc] += add_atoms;
555 >            
556 >            done = 1;
557 >            continue;
558 >          }
559 >          
560 >          // The only situation left is when new_atoms > nTarget.  We
561 >          // want to accept this with some probability that dies off the
562 >          // farther we are from nTarget
563 >          
564 >          // roughly:  x = new_atoms - nTarget
565 >          //           Pacc(x) = exp(- a * x)
566 >          // where a = penalty / (average atoms per molecule)
567 >          
568 >          x = (RealType)(new_atoms - nTarget);
569 >          y = myRandom->rand();
570 >          
571 >          if (y < exp(- a * x)) {
572 >            molToProcMap[i] = which_proc;
573 >            atomsPerProc[which_proc] += add_atoms;
574 >            
575 >            done = 1;
576 >            continue;
577 >          } else {
578 >            continue;
579 >          }
580          }
581 <
582 <        // Spray out this nonsense to all other processors:
583 <
584 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
581 >      }
582 >      
583 >      delete myRandom;
584 >      
585 >      // Spray out this nonsense to all other processors:
586 >      
587 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
588      } else {
589 <
590 <        // Listen to your marching orders from processor 0:
591 <
592 <        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
589 >      
590 >      // Listen to your marching orders from processor 0:
591 >      
592 >      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
593      }
594 <
594 >    
595      info->setMolToProcMap(molToProcMap);
596      sprintf(checkPointMsg,
597              "Successfully divided the molecules among the processors.\n");
598 <    MPIcheckPoint();
599 < }
600 <
598 >    errorCheckPoint();
599 >  }
600 >  
601   #endif
602 <
603 < void SimCreator::createMolecules(SimInfo *info) {
602 >  
603 >  void SimCreator::createMolecules(SimInfo *info) {
604      MoleculeCreator molCreator;
605      int stampId;
606 <
606 >    
607      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
608 <
608 >      
609   #ifdef IS_MPI
610 <
611 <        if (info->getMolToProc(i) == worldRank) {
610 >      
611 >      if (info->getMolToProc(i) == worldRank) {
612   #endif
613 <
614 <            stampId = info->getMoleculeStampId(i);
615 <            Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
616 <                                                                                    stampId, i, info->getLocalIndexManager());
617 <
618 <            info->addMolecule(mol);
619 <
613 >        
614 >        stampId = info->getMoleculeStampId(i);
615 >        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
616 >                                                   stampId, i, info->getLocalIndexManager());
617 >        
618 >        info->addMolecule(mol);
619 >        
620   #ifdef IS_MPI
621 <
622 <        }
623 <
621 >        
622 >      }
623 >      
624   #endif
625 <
625 >      
626      } //end for(int i=0)  
627 < }
628 <
629 < 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) {
627 >  }
628 >    
629 >  void SimCreator::setGlobalIndex(SimInfo *info) {
630      SimInfo::MoleculeIterator mi;
631      Molecule::AtomIterator ai;
632      Molecule::RigidBodyIterator ri;
633      Molecule::CutoffGroupIterator ci;
634 +    Molecule::IntegrableObjectIterator  ioi;
635      Molecule * mol;
636      Atom * atom;
637      RigidBody * rb;
# Line 507 | Line 640 | void SimCreator::setGlobalIndex(SimInfo *info) {
640      int beginRigidBodyIndex;
641      int beginCutoffGroupIndex;
642      int nGlobalAtoms = info->getNGlobalAtoms();
510 –    
511 – #ifndef IS_MPI
643  
644 +    /**@todo fixme */
645 + #ifndef IS_MPI
646 +    
647      beginAtomIndex = 0;
648      beginRigidBodyIndex = 0;
649      beginCutoffGroupIndex = 0;
650 <
650 >    
651   #else
652 <
652 >    
653      int nproc;
654      int myNode;
655 <
655 >    
656      myNode = worldRank;
657      MPI_Comm_size(MPI_COMM_WORLD, &nproc);
658 <
658 >    
659      std::vector < int > tmpAtomsInProc(nproc, 0);
660      std::vector < int > tmpRigidBodiesInProc(nproc, 0);
661      std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
662      std::vector < int > NumAtomsInProc(nproc, 0);
663      std::vector < int > NumRigidBodiesInProc(nproc, 0);
664      std::vector < int > NumCutoffGroupsInProc(nproc, 0);
665 <
665 >    
666      tmpAtomsInProc[myNode] = info->getNAtoms();
667      tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
668      tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
669 <
669 >    
670      //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
671      MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
672                    MPI_SUM, MPI_COMM_WORLD);
# Line 540 | Line 674 | void SimCreator::setGlobalIndex(SimInfo *info) {
674                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
675      MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
676                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
677 <
677 >    
678      beginAtomIndex = 0;
679      beginRigidBodyIndex = 0;
680      beginCutoffGroupIndex = 0;
681 <
681 >    
682      for(int i = 0; i < myNode; i++) {
683 <        beginAtomIndex += NumAtomsInProc[i];
684 <        beginRigidBodyIndex += NumRigidBodiesInProc[i];
685 <        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
683 >      beginAtomIndex += NumAtomsInProc[i];
684 >      beginRigidBodyIndex += NumRigidBodiesInProc[i];
685 >      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
686      }
687 <
687 >    
688   #endif
689 <
689 >    
690      //rigidbody's index begins right after atom's
691      beginRigidBodyIndex += info->getNGlobalAtoms();
692 <
692 >    
693      for(mol = info->beginMolecule(mi); mol != NULL;
694          mol = info->nextMolecule(mi)) {
695 <
696 <        //local index(index in DataStorge) of atom is important
697 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
698 <            atom->setGlobalIndex(beginAtomIndex++);
699 <        }
700 <
701 <        for(rb = mol->beginRigidBody(ri); rb != NULL;
702 <            rb = mol->nextRigidBody(ri)) {
703 <            rb->setGlobalIndex(beginRigidBodyIndex++);
704 <        }
705 <
706 <        //local index of cutoff group is trivial, it only depends on the order of travesing
707 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
708 <            cg = mol->nextCutoffGroup(ci)) {
709 <            cg->setGlobalIndex(beginCutoffGroupIndex++);
710 <        }
695 >      
696 >      //local index(index in DataStorge) of atom is important
697 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
698 >        atom->setGlobalIndex(beginAtomIndex++);
699 >      }
700 >      
701 >      for(rb = mol->beginRigidBody(ri); rb != NULL;
702 >          rb = mol->nextRigidBody(ri)) {
703 >        rb->setGlobalIndex(beginRigidBodyIndex++);
704 >      }
705 >      
706 >      //local index of cutoff group is trivial, it only depends on the order of travesing
707 >      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
708 >          cg = mol->nextCutoffGroup(ci)) {
709 >        cg->setGlobalIndex(beginCutoffGroupIndex++);
710 >      }
711      }
712 <
712 >    
713      //fill globalGroupMembership
714      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
715      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
716 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
717 <
718 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
719 <                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
720 <            }
721 <
722 <        }      
716 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
717 >        
718 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
719 >          globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
720 >        }
721 >        
722 >      }      
723      }
724 <
724 >    
725   #ifdef IS_MPI    
726      // Since the globalGroupMembership has been zero filled and we've only
727      // poked values into the atoms we know, we can do an Allreduce
728      // to get the full globalGroupMembership array (We think).
729      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
730      // docs said we could.
731 <    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
731 >    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
732      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
733                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
734 <     info->setGlobalGroupMembership(tmpGroupMembership);
734 >    info->setGlobalGroupMembership(tmpGroupMembership);
735   #else
736      info->setGlobalGroupMembership(globalGroupMembership);
737   #endif
738 <
738 >    
739      //fill molMembership
740      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
741      
742      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
743 <
744 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
745 <            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
612 <        }
743 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
744 >        globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
745 >      }
746      }
747 <
747 >    
748   #ifdef IS_MPI
749 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
750 <
749 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
750 >    
751      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
752                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
753      
# Line 623 | Line 756 | void SimCreator::setGlobalIndex(SimInfo *info) {
756      info->setGlobalMolMembership(globalMolMembership);
757   #endif
758  
759 < }
759 >    // nIOPerMol holds the number of integrable objects per molecule
760 >    // here the molecules are listed by their global indices.
761  
762 < void SimCreator::loadCoordinates(SimInfo* info) {
762 >    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
763 >    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
764 >      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
765 >    }
766 >    
767 > #ifdef IS_MPI
768 >    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
769 >    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
770 >                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
771 > #else
772 >    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
773 > #endif    
774 >
775 >    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
776 >    
777 >    int startingIndex = 0;
778 >    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
779 >      startingIOIndexForMol[i] = startingIndex;
780 >      startingIndex += numIntegrableObjectsPerMol[i];
781 >    }
782 >    
783 >    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
784 >    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
785 >      int myGlobalIndex = mol->getGlobalIndex();
786 >      int globalIO = startingIOIndexForMol[myGlobalIndex];
787 >      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
788 >           integrableObject = mol->nextIntegrableObject(ioi)) {
789 >        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
790 >        IOIndexToIntegrableObject[globalIO] = integrableObject;
791 >        globalIO++;
792 >      }
793 >    }
794 >    
795 >    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
796 >    
797 >  }
798 >  
799 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
800      Globals* simParams;
801      simParams = info->getSimParams();
802      
803 <    if (!simParams->haveInitialConfig()) {
804 <        sprintf(painCave.errMsg,
634 <                "Cannot intialize a simulation without an initial configuration file.\n");
635 <        painCave.isFatal = 1;;
636 <        simError();
637 <    }
638 <        
639 <    DumpReader reader(info, simParams->getInitialConfig());
803 >    
804 >    DumpReader reader(info, mdFileName);
805      int nframes = reader.getNFrames();
806 <
806 >    
807      if (nframes > 0) {
808 <        reader.readFrame(nframes - 1);
808 >      reader.readFrame(nframes - 1);
809      } else {
810 <        //invalid initial coordinate file
811 <        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
812 <                simParams->getInitialConfig());
813 <        painCave.isFatal = 1;
814 <        simError();
810 >      //invalid initial coordinate file
811 >      sprintf(painCave.errMsg,
812 >              "Initial configuration file %s should at least contain one frame\n",
813 >              mdFileName.c_str());
814 >      painCave.isFatal = 1;
815 >      simError();
816      }
817 <
817 >    
818      //copy the current snapshot to previous snapshot
819      info->getSnapshotManager()->advance();
820 < }
820 >  }
821 >  
822 > } //end namespace OpenMD
823  
656 – } //end namespace oopse
824  
658 –

Diff Legend

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