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 403 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
Revision 1810 by gezelter, Thu Nov 8 14:23:43 2012 UTC

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

Comparing trunk/src/brains/SimCreator.cpp (property svn:keywords):
Revision 403 by gezelter, Tue Mar 8 21:06:49 2005 UTC vs.
Revision 1810 by gezelter, Thu Nov 8 14:23:43 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines