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

Comparing:
trunk/src/brains/SimCreator.cpp (file contents), Revision 403 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
branches/development/src/brains/SimCreator.cpp (file contents), Revision 1841 by gezelter, Mon Jan 28 21:21:35 2013 UTC

# 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   /**
44   * @file SimCreator.cpp
45   * @author tlin
46   * @date 11/03/2004
46 * @time 13:51am
47   * @version 1.0
48   */
49 + #include <exception>
50 + #include <iostream>
51 + #include <sstream>
52 + #include <string>
53  
54   #include "brains/MoleculeCreator.hpp"
55   #include "brains/SimCreator.hpp"
56   #include "brains/SimSnapshotManager.hpp"
57   #include "io/DumpReader.hpp"
58 < #include "io/parse_me.h"
55 < #include "UseTheForce/ForceFieldFactory.hpp"
58 > #include "brains/ForceField.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61   #include "math/SeqRandNumGen.hpp"
62 < #ifdef IS_MPI
63 < #include "io/mpiBASS.h"
64 < #include "math/ParallelRandNumGen.hpp"
65 < #endif
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  
78 <  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){
78 > #include "types/DirectionalAdapter.hpp"
79 > #include "types/MultipoleAdapter.hpp"
80 > #include "types/EAMAdapter.hpp"
81 > #include "types/SuttonChenAdapter.hpp"
82 > #include "types/PolarizableAdapter.hpp"
83 > #include "types/FixedChargeAdapter.hpp"
84 > #include "types/FluctuatingChargeAdapter.hpp"
85  
86   #ifdef IS_MPI
87 + #include "mpi.h"
88 + #include "math/ParallelRandNumGen.hpp"
89 + #endif
90  
91 <    if (worldRank == 0) {
92 < #endif // is_mpi
91 > namespace OpenMD {
92 >  
93 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
94 >    Globals* simParams = NULL;
95 >    try {
96  
97 <      simParams->initalize();
98 <      set_interface_stamps(stamps, simParams);
97 >      // Create a preprocessor that preprocesses md file into an ostringstream
98 >      std::stringstream ppStream;
99 > #ifdef IS_MPI            
100 >      int streamSize;
101 >      const int masterNode = 0;
102  
103 < #ifdef IS_MPI
103 >      if (worldRank == masterNode) {
104 >        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
105 > #endif                
106 >        SimplePreprocessor preprocessor;
107 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
108 >                
109 > #ifdef IS_MPI            
110 >        //brocasting the stream size
111 >        streamSize = ppStream.str().size() +1;
112 >        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
113 >        MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI::CHAR, masterNode);
114 >                          
115 >      } else {
116 >        MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
117  
118 <      mpiEventInit();
118 >        //get stream size
119 >        MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
120  
121 < #endif
121 >        char* buf = new char[streamSize];
122 >        assert(buf);
123 >                
124 >        //receive file content
125 >        MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode);
126 >                
127 >        ppStream.str(buf);
128 >        delete [] buf;
129 >      }
130 > #endif            
131 >      // Create a scanner that reads from the input stream
132 >      MDLexer lexer(ppStream);
133 >      lexer.setFilename(filename);
134 >      lexer.initDeferredLineCount();
135 >    
136 >      // Create a parser that reads from the scanner
137 >      MDParser parser(lexer);
138 >      parser.setFilename(filename);
139  
140 <      yacc_BASS(mdFileName.c_str());
140 >      // Create an observer that synchorizes file name change
141 >      FilenameObserver observer;
142 >      observer.setLexer(&lexer);
143 >      observer.setParser(&parser);
144 >      lexer.setObserver(&observer);
145 >    
146 >      antlr::ASTFactory factory;
147 >      parser.initializeASTFactory(factory);
148 >      parser.setASTFactory(&factory);
149 >      parser.mdfile();
150  
151 < #ifdef IS_MPI
151 >      // Create a tree parser that reads information into Globals
152 >      MDTreeParser treeParser;
153 >      treeParser.initializeASTFactory(factory);
154 >      treeParser.setASTFactory(&factory);
155 >      simParams = treeParser.walkTree(parser.getAST());
156 >    }
157  
158 <      throwMPIEvent(NULL);
159 <    } else {
160 <      set_interface_stamps(stamps, simParams);
161 <      mpiEventInit();
162 <      MPIcheckPoint();
163 <      mpiEventLoop();
158 >      
159 >    catch(antlr::MismatchedCharException& 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::MismatchedTokenException &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::NoViableAltForCharException &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 +    catch(antlr::NoViableAltException &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::TokenStreamRecognitionException& e) {
189 +      sprintf(painCave.errMsg,
190 +              "parser exception: %s %s:%d:%d\n",
191 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
192 +      painCave.isFatal = 1;
193 +      simError();  
194 +    }
195 +        
196 +    catch(antlr::TokenStreamIOException& e) {
197 +      sprintf(painCave.errMsg,
198 +              "parser exception: %s\n",
199 +              e.getMessage().c_str());
200 +      painCave.isFatal = 1;
201 +      simError();
202 +    }
203 +        
204 +    catch(antlr::TokenStreamException& e) {
205 +      sprintf(painCave.errMsg,
206 +              "parser exception: %s\n",
207 +              e.getMessage().c_str());
208 +      painCave.isFatal = 1;
209 +      simError();
210 +    }        
211 +    catch (antlr::RecognitionException& e) {
212 +      sprintf(painCave.errMsg,
213 +              "parser exception: %s %s:%d:%d\n",
214 +              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
215 +      painCave.isFatal = 1;
216 +      simError();          
217 +    }
218 +    catch (antlr::CharStreamException& e) {
219 +      sprintf(painCave.errMsg,
220 +              "parser exception: %s\n",
221 +              e.getMessage().c_str());
222 +      painCave.isFatal = 1;
223 +      simError();        
224 +    }
225 +    catch (OpenMDException& e) {
226 +      sprintf(painCave.errMsg,
227 +              "%s\n",
228 +              e.getMessage().c_str());
229 +      painCave.isFatal = 1;
230 +      simError();
231 +    }
232 +    catch (std::exception& e) {
233 +      sprintf(painCave.errMsg,
234 +              "parser exception: %s\n",
235 +              e.what());
236 +      painCave.isFatal = 1;
237 +      simError();
238 +    }
239  
240 < #endif
241 <
240 >    simParams->setMDfileVersion(mdFileVersion);
241 >    return simParams;
242    }
243 <
244 <  SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
243 >  
244 >  SimInfo*  SimCreator::createSim(const std::string & mdFileName,
245 >                                  bool loadInitCoords) {
246      
247 <    MakeStamps * stamps = new MakeStamps();
247 >    const int bufferSize = 65535;
248 >    char buffer[bufferSize];
249 >    int lineNo = 0;
250 >    std::string mdRawData;
251 >    int metaDataBlockStart = -1;
252 >    int metaDataBlockEnd = -1;
253 >    int i, j;
254 >    streamoff mdOffset;
255 >    int mdFileVersion;
256  
257 <    Globals * simParams = new Globals();
257 >    // Create a string for embedding the version information in the MetaData
258 >    std::string version;
259 >    version.assign("## Last run using OpenMD Version: ");
260 >    version.append(OPENMD_VERSION_MAJOR);
261 >    version.append(".");
262 >    version.append(OPENMD_VERSION_MINOR);
263  
264 <    //parse meta-data file
265 <    parseFile(mdFileName, stamps, simParams);
264 >    std::string svnrev;
265 >    //convert a macro from compiler to a string in c++
266 >    STR_DEFINE(svnrev, SVN_REV );
267 >    version.append(" Revision: ");
268 >    // If there's no SVN revision, just call this the RELEASE revision.
269 >    if (!svnrev.empty()) {
270 >      version.append(svnrev);
271 >    } else {
272 >      version.append("RELEASE");
273 >    }
274 >  
275 > #ifdef IS_MPI            
276 >    const int masterNode = 0;
277 >    if (worldRank == masterNode) {
278 > #endif
279  
280 <    //create the force field
281 <    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
282 <                                                                         simParams->getForceField());
280 >      std::ifstream mdFile_;
281 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
282 >      
283 >      if (mdFile_.fail()) {
284 >        sprintf(painCave.errMsg,
285 >                "SimCreator: Cannot open file: %s\n",
286 >                mdFileName.c_str());
287 >        painCave.isFatal = 1;
288 >        simError();
289 >      }
290 >
291 >      mdFile_.getline(buffer, bufferSize);
292 >      ++lineNo;
293 >      std::string line = trimLeftCopy(buffer);
294 >      i = CaseInsensitiveFind(line, "<OpenMD");
295 >      if (static_cast<size_t>(i) == string::npos) {
296 >        // try the older file strings to see if that works:
297 >        i = CaseInsensitiveFind(line, "<OOPSE");
298 >      }
299 >      
300 >      if (static_cast<size_t>(i) == string::npos) {
301 >        // still no luck!
302 >        sprintf(painCave.errMsg,
303 >                "SimCreator: File: %s is not a valid OpenMD file!\n",
304 >                mdFileName.c_str());
305 >        painCave.isFatal = 1;
306 >        simError();
307 >      }
308 >      
309 >      // found the correct opening string, now try to get the file
310 >      // format version number.
311 >
312 >      StringTokenizer tokenizer(line, "=<> \t\n\r");
313 >      std::string fileType = tokenizer.nextToken();
314 >      toUpper(fileType);
315 >
316 >      mdFileVersion = 0;
317 >
318 >      if (fileType == "OPENMD") {
319 >        while (tokenizer.hasMoreTokens()) {
320 >          std::string token(tokenizer.nextToken());
321 >          toUpper(token);
322 >          if (token == "VERSION") {
323 >            mdFileVersion = tokenizer.nextTokenAsInt();
324 >            break;
325 >          }
326 >        }
327 >      }
328 >            
329 >      //scan through the input stream and find MetaData tag        
330 >      while(mdFile_.getline(buffer, bufferSize)) {
331 >        ++lineNo;
332 >        
333 >        std::string line = trimLeftCopy(buffer);
334 >        if (metaDataBlockStart == -1) {
335 >          i = CaseInsensitiveFind(line, "<MetaData>");
336 >          if (i != string::npos) {
337 >            metaDataBlockStart = lineNo;
338 >            mdOffset = mdFile_.tellg();
339 >          }
340 >        } else {
341 >          i = CaseInsensitiveFind(line, "</MetaData>");
342 >          if (i != string::npos) {
343 >            metaDataBlockEnd = lineNo;
344 >          }
345 >        }
346 >      }
347 >
348 >      if (metaDataBlockStart == -1) {
349 >        sprintf(painCave.errMsg,
350 >                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
351 >                mdFileName.c_str());
352 >        painCave.isFatal = 1;
353 >        simError();
354 >      }
355 >      if (metaDataBlockEnd == -1) {
356 >        sprintf(painCave.errMsg,
357 >                "SimCreator: File: %s did not contain a closed MetaData block!\n",
358 >                mdFileName.c_str());
359 >        painCave.isFatal = 1;
360 >        simError();
361 >      }
362 >        
363 >      mdFile_.clear();
364 >      mdFile_.seekg(0);
365 >      mdFile_.seekg(mdOffset);
366 >
367 >      mdRawData.clear();
368 >
369 >      bool foundVersion = false;
370 >
371 >      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
372 >        mdFile_.getline(buffer, bufferSize);
373 >        std::string line = trimLeftCopy(buffer);
374 >        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version");
375 >        if (static_cast<size_t>(j) != string::npos) {
376 >          foundVersion = true;
377 >          mdRawData += version;
378 >        } else {
379 >          mdRawData += buffer;
380 >        }
381 >        mdRawData += "\n";
382 >      }
383 >      
384 >      if (!foundVersion) mdRawData += version + "\n";
385 >      
386 >      mdFile_.close();
387 >
388 > #ifdef IS_MPI
389 >    }
390 > #endif
391 >
392 >    std::stringstream rawMetaDataStream(mdRawData);
393 >
394 >    //parse meta-data file
395 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
396 >                                   metaDataBlockStart + 1);
397      
398 +    //create the force field
399 +    ForceField * ff = new ForceField(simParams->getForceField());
400 +
401      if (ff == NULL) {
402 <      sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
403 <              simParams->getForceField());
402 >      sprintf(painCave.errMsg,
403 >              "ForceField Factory can not create %s force field\n",
404 >              simParams->getForceField().c_str());
405        painCave.isFatal = 1;
406        simError();
407      }
408 <
408 >    
409      if (simParams->haveForceFieldFileName()) {
410        ff->setForceFieldFileName(simParams->getForceFieldFileName());
411      }
412      
413      std::string forcefieldFileName;
414      forcefieldFileName = ff->getForceFieldFileName();
415 <
415 >    
416      if (simParams->haveForceFieldVariant()) {
417        //If the force field has variant, the variant force field name will be
418        //Base.variant.frc. For exampel EAM.u6.frc
419 <        
419 >      
420        std::string variant = simParams->getForceFieldVariant();
421 <
421 >      
422        std::string::size_type pos = forcefieldFileName.rfind(".frc");
423        variant = "." + variant;
424        if (pos != std::string::npos) {
# Line 139 | Line 430 | namespace oopse {
430      }
431      
432      ff->parse(forcefieldFileName);
142    
143    //extract the molecule stamps
144    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
145    compList(stamps, simParams, moleculeStampPairs);
146
433      //create SimInfo
434 <    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
434 >    SimInfo * info = new SimInfo(ff, simParams);
435  
436 <    //gather parameters (SimCreator only retrieves part of the parameters)
436 >    info->setRawMetaData(mdRawData);
437 >    
438 >    //gather parameters (SimCreator only retrieves part of the
439 >    //parameters)
440      gatherParameters(info, mdFileName);
441 <
441 >    
442      //divide the molecules and determine the global index of molecules
443   #ifdef IS_MPI
444      divideMolecules(info);
445   #endif
446 <
446 >    
447      //create the molecules
448      createMolecules(info);
449 +    
450 +    //find the storage layout
451  
452 +    int storageLayout = computeStorageLayout(info);
453  
454 <    //allocate memory for DataStorage(circular reference, need to break it)
455 <    info->setSnapshotManager(new SimSnapshotManager(info));
454 >    //allocate memory for DataStorage(circular reference, need to
455 >    //break it)
456 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
457      
458 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
459 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
460 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
458 >    //set the global index of atoms, rigidbodies and cutoffgroups
459 >    //(only need to be set once, the global index will never change
460 >    //again). Local indices of atoms and rigidbodies are already set
461 >    //by MoleculeCreator class which actually delegates the
462 >    //responsibility to LocalIndexManager.
463      setGlobalIndex(info);
464 <
465 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
466 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
467 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
468 <    //we can determine the beginning global indices of atoms before they get created.
464 >    
465 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
466 >    //method, at that point atoms don't have the global index yet
467 >    //(their global index are all initialized to -1).  Therefore we
468 >    //have to call addInteractionPairs explicitly here. A way to work
469 >    //around is that we can determine the beginning global indices of
470 >    //atoms before they get created.
471      SimInfo::MoleculeIterator mi;
472      Molecule* mol;
473      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
474 <      info->addExcludePairs(mol);
474 >      info->addInteractionPairs(mol);
475      }
476      
180
181    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
182    //eta, chi for NPT integrator)
477      if (loadInitCoords)
478 <      loadCoordinates(info);    
185 <    
478 >      loadCoordinates(info, mdFileName);    
479      return info;
480    }
481 <
481 >  
482    void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
483 <
484 <    //figure out the ouput file names
483 >    
484 >    //figure out the output file names
485      std::string prefix;
486 <
486 >    
487   #ifdef IS_MPI
488 <
488 >    
489      if (worldRank == 0) {
490   #endif // is_mpi
491        Globals * simParams = info->getSimParams();
# Line 201 | Line 494 | namespace oopse {
494        } else {
495          prefix = getPrefix(mdfile);
496        }
497 <
497 >      
498        info->setFinalConfigFileName(prefix + ".eor");
499        info->setDumpFileName(prefix + ".dump");
500        info->setStatFileName(prefix + ".stat");
501 <
501 >      info->setRestFileName(prefix + ".zang");
502 >      
503   #ifdef IS_MPI
504 <
504 >      
505      }
506 <
506 >    
507   #endif
508 <
508 >    
509    }
510 <
510 >  
511   #ifdef IS_MPI
512    void SimCreator::divideMolecules(SimInfo *info) {
513 <    double numerator;
514 <    double denominator;
515 <    double precast;
516 <    double x;
517 <    double y;
518 <    double a;
513 >    RealType numerator;
514 >    RealType denominator;
515 >    RealType precast;
516 >    RealType x;
517 >    RealType y;
518 >    RealType a;
519      int old_atoms;
520      int add_atoms;
521      int new_atoms;
522      int nTarget;
523      int done;
524      int i;
231    int j;
525      int loops;
526      int which_proc;
527      int nProcessors;
# Line 236 | Line 529 | namespace oopse {
529      int nGlobalMols = info->getNGlobalMolecules();
530      std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
531      
532 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
533 <
532 >    nProcessors = MPI::COMM_WORLD.Get_size();
533 >    
534      if (nProcessors > nGlobalMols) {
535        sprintf(painCave.errMsg,
536                "nProcessors (%d) > nMol (%d)\n"
# Line 245 | Line 538 | namespace oopse {
538                "\tthe number of molecules.  This will not result in a \n"
539                "\tusable division of atoms for force decomposition.\n"
540                "\tEither try a smaller number of processors, or run the\n"
541 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
542 <
541 >              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
542 >      
543        painCave.isFatal = 1;
544        simError();
545      }
546 <
546 >    
547      int seedValue;
548      Globals * simParams = info->getSimParams();
549      SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
# Line 260 | Line 553 | namespace oopse {
553      }else {
554        myRandom = new SeqRandNumGen();
555      }  
556 <
557 <
556 >    
557 >    
558      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
559 <
559 >    
560      //initialize atomsPerProc
561      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
562 <
562 >    
563      if (worldRank == 0) {
564        numerator = info->getNGlobalAtoms();
565        denominator = nProcessors;
566        precast = numerator / denominator;
567        nTarget = (int)(precast + 0.5);
568 <
568 >      
569        for(i = 0; i < nGlobalMols; i++) {
570 +
571          done = 0;
572          loops = 0;
573 <
573 >        
574          while (!done) {
575            loops++;
576 <
576 >          
577            // Pick a processor at random
578 <
578 >          
579            which_proc = (int) (myRandom->rand() * nProcessors);
580 <
580 >          
581            //get the molecule stamp first
582            int stampId = info->getMoleculeStampId(i);
583            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
584 <
584 >          
585            // How many atoms does this processor have so far?
586            old_atoms = atomsPerProc[which_proc];
587            add_atoms = moleculeStamp->getNAtoms();
588            new_atoms = old_atoms + add_atoms;
589 <
589 >          
590            // If we've been through this loop too many times, we need
591            // to just give up and assign the molecule to this processor
592            // and be done with it.
593 <
593 >          
594            if (loops > 100) {
595 +
596              sprintf(painCave.errMsg,
597 <                    "I've tried 100 times to assign molecule %d to a "
598 <                    " processor, but can't find a good spot.\n"
599 <                    "I'm assigning it at random to processor %d.\n",
597 >                    "There have been 100 attempts to assign molecule %d to an\n"
598 >                    "\tunderworked processor, but there's no good place to\n"
599 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
600                      i, which_proc);
601 <
601 >          
602              painCave.isFatal = 0;
603 +            painCave.severity = OPENMD_INFO;
604              simError();
605 <
605 >            
606              molToProcMap[i] = which_proc;
607              atomsPerProc[which_proc] += add_atoms;
608 <
608 >            
609              done = 1;
610              continue;
611            }
612 <
612 >          
613            // If we can add this molecule to this processor without sending
614            // it above nTarget, then go ahead and do it:
615 <
615 >          
616            if (new_atoms <= nTarget) {
617              molToProcMap[i] = which_proc;
618              atomsPerProc[which_proc] += add_atoms;
619 <
619 >            
620              done = 1;
621              continue;
622            }
623 <
623 >          
624            // The only situation left is when new_atoms > nTarget.  We
625            // want to accept this with some probability that dies off the
626            // farther we are from nTarget
627 <
627 >          
628            // roughly:  x = new_atoms - nTarget
629            //           Pacc(x) = exp(- a * x)
630            // where a = penalty / (average atoms per molecule)
631 <
632 <          x = (double)(new_atoms - nTarget);
631 >          
632 >          x = (RealType)(new_atoms - nTarget);
633            y = myRandom->rand();
634 <
634 >          
635            if (y < exp(- a * x)) {
636              molToProcMap[i] = which_proc;
637              atomsPerProc[which_proc] += add_atoms;
638 <
638 >            
639              done = 1;
640              continue;
641            } else {
# Line 347 | Line 643 | namespace oopse {
643            }
644          }
645        }
646 <
646 >      
647        delete myRandom;
352        
353      // Spray out this nonsense to all other processors:
648  
649 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
649 >      // Spray out this nonsense to all other processors:
650 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
651      } else {
652 <
652 >      
653        // Listen to your marching orders from processor 0:
654 +      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
655  
360      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
656      }
657 <
657 >    
658      info->setMolToProcMap(molToProcMap);
659      sprintf(checkPointMsg,
660              "Successfully divided the molecules among the processors.\n");
661 <    MPIcheckPoint();
661 >    errorCheckPoint();
662    }
663 <
663 >  
664   #endif
665 <
665 >  
666    void SimCreator::createMolecules(SimInfo *info) {
667      MoleculeCreator molCreator;
668      int stampId;
669 <
669 >    
670      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
671 <
671 >      
672   #ifdef IS_MPI
673 <
673 >      
674        if (info->getMolToProc(i) == worldRank) {
675   #endif
676 <
676 >        
677          stampId = info->getMoleculeStampId(i);
678 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
679 <                                                   stampId, i, info->getLocalIndexManager());
680 <
678 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
679 >                                                   info->getMoleculeStamp(stampId),
680 >                                                   stampId, i,
681 >                                                   info->getLocalIndexManager());
682 >        
683          info->addMolecule(mol);
684 <
684 >        
685   #ifdef IS_MPI
686 <
686 >        
687        }
688 <
688 >      
689   #endif
690 <
690 >      
691      } //end for(int i=0)  
692    }
693 +    
694 +  int SimCreator::computeStorageLayout(SimInfo* info) {
695  
696 <  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
697 <                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
698 <    int i;
699 <    char * id;
700 <    LinkedMolStamp* extractedStamp = NULL;
701 <    MoleculeStamp * currentStamp;
702 <    Component** the_components = simParams->getComponents();
703 <    int n_components = simParams->getNComponents();
696 >    Globals* simParams = info->getSimParams();
697 >    int nRigidBodies = info->getNGlobalRigidBodies();
698 >    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
699 >    set<AtomType*>::iterator i;
700 >    bool hasDirectionalAtoms = false;
701 >    bool hasFixedCharge = false;
702 >    bool hasDipoles = false;    
703 >    bool hasQuadrupoles = false;    
704 >    bool hasPolarizable = false;    
705 >    bool hasFluctuatingCharge = false;    
706 >    bool hasMetallic = false;
707 >    int storageLayout = 0;
708 >    storageLayout |= DataStorage::dslPosition;
709 >    storageLayout |= DataStorage::dslVelocity;
710 >    storageLayout |= DataStorage::dslForce;
711  
712 <    if (!simParams->haveNMol()) {
407 <      // we don't have the total number of molecules, so we assume it is
408 <      // given in each component
712 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
713  
714 <      for(i = 0; i < n_components; i++) {
715 <        if (!the_components[i]->haveNMol()) {
716 <          // we have a problem
717 <          sprintf(painCave.errMsg,
718 <                  "SimCreator Error. No global NMol or component NMol given.\n"
719 <                  "\tCannot calculate the number of atoms.\n");
714 >      DirectionalAdapter da = DirectionalAdapter( (*i) );
715 >      MultipoleAdapter ma = MultipoleAdapter( (*i) );
716 >      EAMAdapter ea = EAMAdapter( (*i) );
717 >      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
718 >      PolarizableAdapter pa = PolarizableAdapter( (*i) );
719 >      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
720 >      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
721  
722 <          painCave.isFatal = 1;
723 <          simError();
724 <        }
725 <      
726 <        id = the_components[i]->getType();
722 >      if (da.isDirectional()){
723 >        hasDirectionalAtoms = true;
724 >      }
725 >      if (ma.isDipole()){
726 >        hasDipoles = true;
727 >      }
728 >      if (ma.isQuadrupole()){
729 >        hasQuadrupoles = true;
730 >      }
731 >      if (ea.isEAM() || sca.isSuttonChen()){
732 >        hasMetallic = true;
733 >      }
734 >      if ( fca.isFixedCharge() ){
735 >        hasFixedCharge = true;
736 >      }
737 >      if ( fqa.isFluctuatingCharge() ){
738 >        hasFluctuatingCharge = true;
739 >      }
740 >      if ( pa.isPolarizable() ){
741 >        hasPolarizable = true;
742 >      }
743 >    }
744 >    
745 >    if (nRigidBodies > 0 || hasDirectionalAtoms) {
746 >      storageLayout |= DataStorage::dslAmat;
747 >      if(storageLayout & DataStorage::dslVelocity) {
748 >        storageLayout |= DataStorage::dslAngularMomentum;
749 >      }
750 >      if (storageLayout & DataStorage::dslForce) {
751 >        storageLayout |= DataStorage::dslTorque;
752 >      }
753 >    }
754 >    if (hasDipoles) {
755 >      storageLayout |= DataStorage::dslDipole;
756 >    }
757 >    if (hasQuadrupoles) {
758 >      storageLayout |= DataStorage::dslQuadrupole;
759 >    }
760 >    if (hasFixedCharge || hasFluctuatingCharge) {
761 >      storageLayout |= DataStorage::dslSkippedCharge;
762 >    }
763 >    if (hasMetallic) {
764 >      storageLayout |= DataStorage::dslDensity;
765 >      storageLayout |= DataStorage::dslFunctional;
766 >      storageLayout |= DataStorage::dslFunctionalDerivative;
767 >    }
768 >    if (hasPolarizable) {
769 >      storageLayout |= DataStorage::dslElectricField;
770 >    }
771 >    if (hasFluctuatingCharge){
772 >      storageLayout |= DataStorage::dslFlucQPosition;
773 >      if(storageLayout & DataStorage::dslVelocity) {
774 >        storageLayout |= DataStorage::dslFlucQVelocity;
775 >      }
776 >      if (storageLayout & DataStorage::dslForce) {
777 >        storageLayout |= DataStorage::dslFlucQForce;
778 >      }
779 >    }
780 >    
781 >    // if the user has asked for them, make sure we've got the memory for the
782 >    // objects defined.
783  
784 <        extractedStamp = stamps->extractMolStamp(id);
785 <        if (extractedStamp == NULL) {
786 <          sprintf(painCave.errMsg,
426 <                  "SimCreator error: Component \"%s\" was not found in the "
427 <                  "list of declared molecules\n", id);
784 >    if (simParams->getOutputParticlePotential()) {
785 >      storageLayout |= DataStorage::dslParticlePot;
786 >    }
787  
788 <          painCave.isFatal = 1;
789 <          simError();
790 <        }
788 >    if (simParams->havePrintHeatFlux()) {
789 >      if (simParams->getPrintHeatFlux()) {
790 >        storageLayout |= DataStorage::dslParticlePot;
791 >      }
792 >    }
793  
794 <        currentStamp = extractedStamp->getStamp();
794 >    if (simParams->getOutputElectricField() | simParams->haveElectricField()) {
795 >      storageLayout |= DataStorage::dslElectricField;
796 >    }
797  
798 <
799 <        moleculeStampPairs.push_back(
800 <                                     std::make_pair(currentStamp, the_components[i]->getNMol()));
801 <      } //end for (i = 0; i < n_components; i++)
439 <    } else {
440 <      sprintf(painCave.errMsg, "SimSetup error.\n"
441 <              "\tSorry, the ability to specify total"
442 <              " nMols and then give molfractions in the components\n"
443 <              "\tis not currently supported."
444 <              " Please give nMol in the components.\n");
445 <
446 <      painCave.isFatal = 1;
447 <      simError();
798 >    if (simParams->getOutputFluctuatingCharges()) {
799 >      storageLayout |= DataStorage::dslFlucQPosition;
800 >      storageLayout |= DataStorage::dslFlucQVelocity;
801 >      storageLayout |= DataStorage::dslFlucQForce;
802      }
803  
804 < #ifdef IS_MPI
451 <
452 <    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
453 <    MPIcheckPoint();
804 >    info->setStorageLayout(storageLayout);
805  
806 < #endif // is_mpi
456 <
806 >    return storageLayout;
807    }
808  
809    void SimCreator::setGlobalIndex(SimInfo *info) {
# Line 461 | Line 811 | namespace oopse {
811      Molecule::AtomIterator ai;
812      Molecule::RigidBodyIterator ri;
813      Molecule::CutoffGroupIterator ci;
814 +    Molecule::IntegrableObjectIterator  ioi;
815      Molecule * mol;
816      Atom * atom;
817      RigidBody * rb;
# Line 469 | Line 820 | namespace oopse {
820      int beginRigidBodyIndex;
821      int beginCutoffGroupIndex;
822      int nGlobalAtoms = info->getNGlobalAtoms();
823 +    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
824      
473 #ifndef IS_MPI
474
825      beginAtomIndex = 0;
826 <    beginRigidBodyIndex = 0;
826 >    //rigidbody's index begins right after atom's
827 >    beginRigidBodyIndex = info->getNGlobalAtoms();
828      beginCutoffGroupIndex = 0;
829  
830 < #else
830 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
831 >      
832 > #ifdef IS_MPI      
833 >      if (info->getMolToProc(i) == worldRank) {
834 > #endif        
835 >        // stuff to do if I own this molecule
836 >        mol = info->getMoleculeByGlobalIndex(i);
837  
838 <    int nproc;
839 <    int myNode;
838 >        //local index(index in DataStorge) of atom is important
839 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
840 >          atom->setGlobalIndex(beginAtomIndex++);
841 >        }
842 >        
843 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
844 >            rb = mol->nextRigidBody(ri)) {
845 >          rb->setGlobalIndex(beginRigidBodyIndex++);
846 >        }
847 >        
848 >        //local index of cutoff group is trivial, it only depends on
849 >        //the order of travesing
850 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
851 >            cg = mol->nextCutoffGroup(ci)) {
852 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
853 >        }        
854 >        
855 > #ifdef IS_MPI        
856 >      }  else {
857  
858 <    myNode = worldRank;
859 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
858 >        // stuff to do if I don't own this molecule
859 >        
860 >        int stampId = info->getMoleculeStampId(i);
861 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
862  
863 <    std::vector < int > tmpAtomsInProc(nproc, 0);
864 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
865 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
490 <    std::vector < int > NumAtomsInProc(nproc, 0);
491 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
492 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
493 <
494 <    tmpAtomsInProc[myNode] = info->getNAtoms();
495 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
496 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
497 <
498 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
499 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
500 <                  MPI_SUM, MPI_COMM_WORLD);
501 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
502 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
503 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
504 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
505 <
506 <    beginAtomIndex = 0;
507 <    beginRigidBodyIndex = 0;
508 <    beginCutoffGroupIndex = 0;
509 <
510 <    for(int i = 0; i < myNode; i++) {
511 <      beginAtomIndex += NumAtomsInProc[i];
512 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
513 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
514 <    }
515 <
516 < #endif
517 <
518 <    //rigidbody's index begins right after atom's
519 <    beginRigidBodyIndex += info->getNGlobalAtoms();
520 <
521 <    for(mol = info->beginMolecule(mi); mol != NULL;
522 <        mol = info->nextMolecule(mi)) {
523 <
524 <      //local index(index in DataStorge) of atom is important
525 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
526 <        atom->setGlobalIndex(beginAtomIndex++);
863 >        beginAtomIndex += stamp->getNAtoms();
864 >        beginRigidBodyIndex += stamp->getNRigidBodies();
865 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
866        }
867 + #endif          
868  
869 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
530 <          rb = mol->nextRigidBody(ri)) {
531 <        rb->setGlobalIndex(beginRigidBodyIndex++);
532 <      }
869 >    } //end for(int i=0)  
870  
534      //local index of cutoff group is trivial, it only depends on the order of travesing
535      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
536          cg = mol->nextCutoffGroup(ci)) {
537        cg->setGlobalIndex(beginCutoffGroupIndex++);
538      }
539    }
540
871      //fill globalGroupMembership
872      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
873      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
874        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
875 <
875 >        
876          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
877            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
878          }
879 <
879 >        
880        }      
881      }
882 <
882 >  
883   #ifdef IS_MPI    
884      // Since the globalGroupMembership has been zero filled and we've only
885      // poked values into the atoms we know, we can do an Allreduce
886      // to get the full globalGroupMembership array (We think).
887      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
888      // docs said we could.
889 <    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
890 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
891 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
889 >    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
890 >    MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
891 >                              &tmpGroupMembership[0], nGlobalAtoms,
892 >                              MPI::INT, MPI::SUM);
893      info->setGlobalGroupMembership(tmpGroupMembership);
894   #else
895      info->setGlobalGroupMembership(globalGroupMembership);
896   #endif
897 <
897 >    
898      //fill molMembership
899 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
899 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
900 >                                         info->getNGlobalRigidBodies(), 0);
901      
902 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
903 <
902 >    for(mol = info->beginMolecule(mi); mol != NULL;
903 >        mol = info->nextMolecule(mi)) {
904        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
905          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
906        }
907 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
908 +           rb = mol->nextRigidBody(ri)) {
909 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
910 +      }
911      }
912 <
912 >    
913   #ifdef IS_MPI
914 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
915 <
916 <    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
917 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
914 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
915 >                                      info->getNGlobalRigidBodies(), 0);
916 >    MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
917 >                              nGlobalAtoms + nGlobalRigidBodies,
918 >                              MPI::INT, MPI::SUM);
919      
920      info->setGlobalMolMembership(tmpMolMembership);
921   #else
922      info->setGlobalMolMembership(globalMolMembership);
923   #endif
924  
925 <  }
925 >    // nIOPerMol holds the number of integrable objects per molecule
926 >    // here the molecules are listed by their global indices.
927  
928 <  void SimCreator::loadCoordinates(SimInfo* info) {
929 <    Globals* simParams;
930 <    simParams = info->getSimParams();
928 >    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
929 >    for (mol = info->beginMolecule(mi); mol != NULL;
930 >         mol = info->nextMolecule(mi)) {
931 >      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
932 >    }
933      
934 <    if (!simParams->haveInitialConfig()) {
935 <      sprintf(painCave.errMsg,
936 <              "Cannot intialize a simulation without an initial configuration file.\n");
937 <      painCave.isFatal = 1;;
938 <      simError();
934 > #ifdef IS_MPI
935 >    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
936 >    MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
937 >                              info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
938 > #else
939 >    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
940 > #endif    
941 >
942 >    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
943 >    
944 >    int startingIndex = 0;
945 >    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
946 >      startingIOIndexForMol[i] = startingIndex;
947 >      startingIndex += numIntegrableObjectsPerMol[i];
948      }
949 <        
950 <    DumpReader reader(info, simParams->getInitialConfig());
949 >    
950 >    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
951 >    for (mol = info->beginMolecule(mi); mol != NULL;
952 >         mol = info->nextMolecule(mi)) {
953 >      int myGlobalIndex = mol->getGlobalIndex();
954 >      int globalIO = startingIOIndexForMol[myGlobalIndex];
955 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
956 >           sd = mol->nextIntegrableObject(ioi)) {
957 >        sd->setGlobalIntegrableObjectIndex(globalIO);
958 >        IOIndexToIntegrableObject[globalIO] = sd;
959 >        globalIO++;
960 >      }
961 >    }
962 >      
963 >    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
964 >    
965 >  }
966 >  
967 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
968 >    
969 >    DumpReader reader(info, mdFileName);
970      int nframes = reader.getNFrames();
971 <
971 >    
972      if (nframes > 0) {
973        reader.readFrame(nframes - 1);
974      } else {
975        //invalid initial coordinate file
976 <      sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
977 <              simParams->getInitialConfig());
976 >      sprintf(painCave.errMsg,
977 >              "Initial configuration file %s should at least contain one frame\n",
978 >              mdFileName.c_str());
979        painCave.isFatal = 1;
980        simError();
981      }
613
982      //copy the current snapshot to previous snapshot
983      info->getSnapshotManager()->advance();
984    }
985 +  
986 + } //end namespace OpenMD
987  
618 } //end namespace oopse
988  
620

Comparing:
trunk/src/brains/SimCreator.cpp (property svn:keywords), Revision 403 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
branches/development/src/brains/SimCreator.cpp (property svn:keywords), Revision 1841 by gezelter, Mon Jan 28 21:21:35 2013 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines