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 1874 by gezelter, Wed May 15 15:09: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, 234107 (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;
225 <    int old_atoms;
226 <    int add_atoms;
227 <    int new_atoms;
228 <    int nTarget;
229 <    int done;
230 <    int i;
231 <    int j;
232 <    int loops;
233 <    int which_proc;
513 >    RealType numerator;
514 >    RealType denominator;
515 >    RealType precast;
516 >    RealType x;
517 >    RealType y;
518 >    RealType a;
519      int nProcessors;
520      std::vector<int> atomsPerProc;
521      int nGlobalMols = info->getNGlobalMolecules();
522 <    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
522 >    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an
523 >                                                    // error
524 >                                                    // condition:
525      
526 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
527 <
526 >    nProcessors = MPI::COMM_WORLD.Get_size();
527 >    
528      if (nProcessors > nGlobalMols) {
529        sprintf(painCave.errMsg,
530                "nProcessors (%d) > nMol (%d)\n"
# Line 245 | Line 532 | namespace oopse {
532                "\tthe number of molecules.  This will not result in a \n"
533                "\tusable division of atoms for force decomposition.\n"
534                "\tEither try a smaller number of processors, or run the\n"
535 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
536 <
535 >              "\tsingle-processor version of OpenMD.\n", nProcessors,
536 >              nGlobalMols);
537 >      
538        painCave.isFatal = 1;
539        simError();
540      }
541 <
254 <    int seedValue;
541 >    
542      Globals * simParams = info->getSimParams();
543 <    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
543 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel
544 >                             //random number generator
545      if (simParams->haveSeed()) {
546 <      seedValue = simParams->getSeed();
546 >      int seedValue = simParams->getSeed();
547        myRandom = new SeqRandNumGen(seedValue);
548      }else {
549        myRandom = new SeqRandNumGen();
550      }  
551 <
552 <
551 >    
552 >    
553      a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
554 <
554 >    
555      //initialize atomsPerProc
556      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
557 <
557 >    
558      if (worldRank == 0) {
559        numerator = info->getNGlobalAtoms();
560        denominator = nProcessors;
561        precast = numerator / denominator;
562 <      nTarget = (int)(precast + 0.5);
562 >      int nTarget = (int)(precast + 0.5);
563 >      
564 >      for(int i = 0; i < nGlobalMols; i++) {
565  
566 <      for(i = 0; i < nGlobalMols; i++) {
567 <        done = 0;
568 <        loops = 0;
279 <
566 >        int done = 0;
567 >        int loops = 0;
568 >        
569          while (!done) {
570            loops++;
571 <
571 >          
572            // Pick a processor at random
573 <
574 <          which_proc = (int) (myRandom->rand() * nProcessors);
575 <
573 >          
574 >          int which_proc = (int) (myRandom->rand() * nProcessors);
575 >          
576            //get the molecule stamp first
577            int stampId = info->getMoleculeStampId(i);
578            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
579 <
579 >          
580            // How many atoms does this processor have so far?
581 <          old_atoms = atomsPerProc[which_proc];
582 <          add_atoms = moleculeStamp->getNAtoms();
583 <          new_atoms = old_atoms + add_atoms;
584 <
581 >          int old_atoms = atomsPerProc[which_proc];
582 >          int add_atoms = moleculeStamp->getNAtoms();
583 >          int new_atoms = old_atoms + add_atoms;
584 >          
585            // If we've been through this loop too many times, we need
586            // to just give up and assign the molecule to this processor
587            // and be done with it.
588 <
588 >          
589            if (loops > 100) {
590 +
591              sprintf(painCave.errMsg,
592 <                    "I've tried 100 times to assign molecule %d to a "
593 <                    " processor, but can't find a good spot.\n"
594 <                    "I'm assigning it at random to processor %d.\n",
592 >                    "There have been 100 attempts to assign molecule %d to an\n"
593 >                    "\tunderworked processor, but there's no good place to\n"
594 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
595                      i, which_proc);
596 <
596 >          
597              painCave.isFatal = 0;
598 +            painCave.severity = OPENMD_INFO;
599              simError();
600 <
600 >            
601              molToProcMap[i] = which_proc;
602              atomsPerProc[which_proc] += add_atoms;
603 <
603 >            
604              done = 1;
605              continue;
606            }
607 <
607 >          
608            // If we can add this molecule to this processor without sending
609            // it above nTarget, then go ahead and do it:
610 <
610 >          
611            if (new_atoms <= nTarget) {
612              molToProcMap[i] = which_proc;
613              atomsPerProc[which_proc] += add_atoms;
614 <
614 >            
615              done = 1;
616              continue;
617            }
618 <
618 >          
619            // The only situation left is when new_atoms > nTarget.  We
620            // want to accept this with some probability that dies off the
621            // farther we are from nTarget
622 <
622 >          
623            // roughly:  x = new_atoms - nTarget
624            //           Pacc(x) = exp(- a * x)
625            // where a = penalty / (average atoms per molecule)
626 <
627 <          x = (double)(new_atoms - nTarget);
626 >          
627 >          x = (RealType)(new_atoms - nTarget);
628            y = myRandom->rand();
629 <
629 >          
630            if (y < exp(- a * x)) {
631              molToProcMap[i] = which_proc;
632              atomsPerProc[which_proc] += add_atoms;
633 <
633 >            
634              done = 1;
635              continue;
636            } else {
# Line 347 | Line 638 | namespace oopse {
638            }
639          }
640        }
641 <
641 >      
642        delete myRandom;
352        
353      // Spray out this nonsense to all other processors:
643  
644 <      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
644 >      // Spray out this nonsense to all other processors:
645 >      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
646      } else {
647 <
647 >      
648        // Listen to your marching orders from processor 0:
649 +      MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
650  
360      MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
651      }
652 <
652 >    
653      info->setMolToProcMap(molToProcMap);
654      sprintf(checkPointMsg,
655              "Successfully divided the molecules among the processors.\n");
656 <    MPIcheckPoint();
656 >    errorCheckPoint();
657    }
658 <
658 >  
659   #endif
660 <
660 >  
661    void SimCreator::createMolecules(SimInfo *info) {
662      MoleculeCreator molCreator;
663      int stampId;
664 <
664 >    
665      for(int i = 0; i < info->getNGlobalMolecules(); i++) {
666 <
666 >      
667   #ifdef IS_MPI
668 <
668 >      
669        if (info->getMolToProc(i) == worldRank) {
670   #endif
671 <
671 >        
672          stampId = info->getMoleculeStampId(i);
673 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
674 <                                                   stampId, i, info->getLocalIndexManager());
675 <
673 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
674 >                                                   info->getMoleculeStamp(stampId),
675 >                                                   stampId, i,
676 >                                                   info->getLocalIndexManager());
677 >        
678          info->addMolecule(mol);
679 <
679 >        
680   #ifdef IS_MPI
681 <
681 >        
682        }
683 <
683 >      
684   #endif
685 <
685 >      
686      } //end for(int i=0)  
687    }
688 +    
689 +  int SimCreator::computeStorageLayout(SimInfo* info) {
690  
691 <  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
692 <                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
693 <    int i;
694 <    char * id;
695 <    LinkedMolStamp* extractedStamp = NULL;
696 <    MoleculeStamp * currentStamp;
697 <    Component** the_components = simParams->getComponents();
698 <    int n_components = simParams->getNComponents();
691 >    Globals* simParams = info->getSimParams();
692 >    int nRigidBodies = info->getNGlobalRigidBodies();
693 >    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
694 >    set<AtomType*>::iterator i;
695 >    bool hasDirectionalAtoms = false;
696 >    bool hasFixedCharge = false;
697 >    bool hasDipoles = false;    
698 >    bool hasQuadrupoles = false;    
699 >    bool hasPolarizable = false;    
700 >    bool hasFluctuatingCharge = false;    
701 >    bool hasMetallic = false;
702 >    int storageLayout = 0;
703 >    storageLayout |= DataStorage::dslPosition;
704 >    storageLayout |= DataStorage::dslVelocity;
705 >    storageLayout |= DataStorage::dslForce;
706  
707 <    if (!simParams->haveNMol()) {
407 <      // we don't have the total number of molecules, so we assume it is
408 <      // given in each component
707 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
708  
709 <      for(i = 0; i < n_components; i++) {
710 <        if (!the_components[i]->haveNMol()) {
711 <          // we have a problem
712 <          sprintf(painCave.errMsg,
713 <                  "SimCreator Error. No global NMol or component NMol given.\n"
714 <                  "\tCannot calculate the number of atoms.\n");
709 >      DirectionalAdapter da = DirectionalAdapter( (*i) );
710 >      MultipoleAdapter ma = MultipoleAdapter( (*i) );
711 >      EAMAdapter ea = EAMAdapter( (*i) );
712 >      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
713 >      PolarizableAdapter pa = PolarizableAdapter( (*i) );
714 >      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
715 >      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
716  
717 <          painCave.isFatal = 1;
718 <          simError();
719 <        }
720 <      
721 <        id = the_components[i]->getType();
717 >      if (da.isDirectional()){
718 >        hasDirectionalAtoms = true;
719 >      }
720 >      if (ma.isDipole()){
721 >        hasDipoles = true;
722 >      }
723 >      if (ma.isQuadrupole()){
724 >        hasQuadrupoles = true;
725 >      }
726 >      if (ea.isEAM() || sca.isSuttonChen()){
727 >        hasMetallic = true;
728 >      }
729 >      if ( fca.isFixedCharge() ){
730 >        hasFixedCharge = true;
731 >      }
732 >      if ( fqa.isFluctuatingCharge() ){
733 >        hasFluctuatingCharge = true;
734 >      }
735 >      if ( pa.isPolarizable() ){
736 >        hasPolarizable = true;
737 >      }
738 >    }
739 >    
740 >    if (nRigidBodies > 0 || hasDirectionalAtoms) {
741 >      storageLayout |= DataStorage::dslAmat;
742 >      if(storageLayout & DataStorage::dslVelocity) {
743 >        storageLayout |= DataStorage::dslAngularMomentum;
744 >      }
745 >      if (storageLayout & DataStorage::dslForce) {
746 >        storageLayout |= DataStorage::dslTorque;
747 >      }
748 >    }
749 >    if (hasDipoles) {
750 >      storageLayout |= DataStorage::dslDipole;
751 >    }
752 >    if (hasQuadrupoles) {
753 >      storageLayout |= DataStorage::dslQuadrupole;
754 >    }
755 >    if (hasFixedCharge || hasFluctuatingCharge) {
756 >      storageLayout |= DataStorage::dslSkippedCharge;
757 >    }
758 >    if (hasMetallic) {
759 >      storageLayout |= DataStorage::dslDensity;
760 >      storageLayout |= DataStorage::dslFunctional;
761 >      storageLayout |= DataStorage::dslFunctionalDerivative;
762 >    }
763 >    if (hasPolarizable) {
764 >      storageLayout |= DataStorage::dslElectricField;
765 >    }
766 >    if (hasFluctuatingCharge){
767 >      storageLayout |= DataStorage::dslFlucQPosition;
768 >      if(storageLayout & DataStorage::dslVelocity) {
769 >        storageLayout |= DataStorage::dslFlucQVelocity;
770 >      }
771 >      if (storageLayout & DataStorage::dslForce) {
772 >        storageLayout |= DataStorage::dslFlucQForce;
773 >      }
774 >    }
775 >    
776 >    // if the user has asked for them, make sure we've got the memory for the
777 >    // objects defined.
778  
779 <        extractedStamp = stamps->extractMolStamp(id);
780 <        if (extractedStamp == NULL) {
781 <          sprintf(painCave.errMsg,
426 <                  "SimCreator error: Component \"%s\" was not found in the "
427 <                  "list of declared molecules\n", id);
779 >    if (simParams->getOutputParticlePotential()) {
780 >      storageLayout |= DataStorage::dslParticlePot;
781 >    }
782  
783 <          painCave.isFatal = 1;
784 <          simError();
785 <        }
783 >    if (simParams->havePrintHeatFlux()) {
784 >      if (simParams->getPrintHeatFlux()) {
785 >        storageLayout |= DataStorage::dslParticlePot;
786 >      }
787 >    }
788  
789 <        currentStamp = extractedStamp->getStamp();
790 <
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();
789 >    if (simParams->getOutputElectricField() | simParams->haveElectricField()) {
790 >      storageLayout |= DataStorage::dslElectricField;
791      }
792  
793 < #ifdef IS_MPI
793 >    if (simParams->getOutputFluctuatingCharges()) {
794 >      storageLayout |= DataStorage::dslFlucQPosition;
795 >      storageLayout |= DataStorage::dslFlucQVelocity;
796 >      storageLayout |= DataStorage::dslFlucQForce;
797 >    }
798  
799 <    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
453 <    MPIcheckPoint();
799 >    info->setStorageLayout(storageLayout);
800  
801 < #endif // is_mpi
456 <
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 <
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.
branches/development/src/brains/SimCreator.cpp (property svn:keywords), Revision 1874 by gezelter, Wed May 15 15:09:35 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines