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 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 2069 by gezelter, Thu Mar 5 16:30:23 2015 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 +
50 + #ifdef IS_MPI
51 + #include "mpi.h"
52 + #include "math/ParallelRandNumGen.hpp"
53 + #endif
54 +
55   #include <exception>
56   #include <iostream>
57   #include <sstream>
# Line 55 | Line 61
61   #include "brains/SimCreator.hpp"
62   #include "brains/SimSnapshotManager.hpp"
63   #include "io/DumpReader.hpp"
64 < #include "UseTheForce/ForceFieldFactory.hpp"
64 > #include "brains/ForceField.hpp"
65   #include "utils/simError.h"
66   #include "utils/StringUtils.hpp"
67 + #include "utils/Revision.hpp"
68   #include "math/SeqRandNumGen.hpp"
69   #include "mdParser/MDLexer.hpp"
70   #include "mdParser/MDParser.hpp"
# Line 75 | Line 82
82   #include "antlr/NoViableAltForCharException.hpp"
83   #include "antlr/NoViableAltException.hpp"
84  
85 < #ifdef IS_MPI
86 < #include "math/ParallelRandNumGen.hpp"
87 < #endif
85 > #include "types/DirectionalAdapter.hpp"
86 > #include "types/MultipoleAdapter.hpp"
87 > #include "types/EAMAdapter.hpp"
88 > #include "types/SuttonChenAdapter.hpp"
89 > #include "types/PolarizableAdapter.hpp"
90 > #include "types/FixedChargeAdapter.hpp"
91 > #include "types/FluctuatingChargeAdapter.hpp"
92  
93 < namespace oopse {
93 >
94 > namespace OpenMD {
95    
96 < Globals* SimCreator::parseFile(const std::string mdFileName){
97 <        Globals* simParams = NULL;
98 <        try {
96 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
97 >    Globals* simParams = NULL;
98 >    try {
99  
100 <            // Create a preprocessor that preprocesses md file into an ostringstream
101 <            std::stringstream ppStream;
100 >      // Create a preprocessor that preprocesses md file into an ostringstream
101 >      std::stringstream ppStream;
102   #ifdef IS_MPI            
103 <            int streamSize;
104 <            const int masterNode = 0;
105 <            int commStatus;
106 <            if (worldRank == masterNode) {
107 < #endif
103 >      int streamSize;
104 >      const int masterNode = 0;
105 >
106 >      if (worldRank == masterNode) {
107 >        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
108 > #endif                
109 >        SimplePreprocessor preprocessor;
110 >        preprocessor.preprocess(rawMetaDataStream, filename,
111 >                                startOfMetaDataBlock, ppStream);
112                  
97                SimplePreprocessor preprocessor;
98                preprocessor.preprocess(mdFileName, ppStream);
99                
113   #ifdef IS_MPI            
114 <                //brocasting the stream size
115 <                streamSize = ppStream.str().size() +1;
116 <                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
114 >        //broadcasting the stream size
115 >        streamSize = ppStream.str().size() +1;
116 >        MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
117 >        MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
118 >                  streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
119 >      } else {
120  
121 <                commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
122 <            
121 >        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
122 >
123 >        //get stream size
124 >        MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
125 >        char* buf = new char[streamSize];
126 >        assert(buf);
127                  
128 <            } else {
129 <                //get stream size
110 <                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
111 <                
112 <                  char* buf = new char[streamSize];
113 <                  assert(buf);
114 <                
115 <                  //receive file content
116 <                  commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
117 <                
118 <                  ppStream.str(buf);
119 <                  delete buf;
128 >        //receive file content
129 >        MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
130  
131 <            }
131 >        ppStream.str(buf);
132 >        delete [] buf;
133 >      }
134   #endif            
135 <            // Create a scanner that reads from the input stream
136 <            MDLexer lexer(ppStream);
137 <            lexer.setFilename(mdFileName);
138 <            lexer.initDeferredLineCount();
135 >      // Create a scanner that reads from the input stream
136 >      MDLexer lexer(ppStream);
137 >      lexer.setFilename(filename);
138 >      lexer.initDeferredLineCount();
139      
140 <            // Create a parser that reads from the scanner
141 <            MDParser parser(lexer);
142 <            parser.setFilename(mdFileName);
140 >      // Create a parser that reads from the scanner
141 >      MDParser parser(lexer);
142 >      parser.setFilename(filename);
143  
144 <            // Create an observer that synchorizes file name change
145 <            FilenameObserver observer;
146 <            observer.setLexer(&lexer);
147 <            observer.setParser(&parser);
148 <            lexer.setObserver(&observer);
144 >      // Create an observer that synchorizes file name change
145 >      FilenameObserver observer;
146 >      observer.setLexer(&lexer);
147 >      observer.setParser(&parser);
148 >      lexer.setObserver(&observer);
149      
150 <            antlr::ASTFactory factory;
151 <            parser.initializeASTFactory(factory);
152 <            parser.setASTFactory(&factory);
153 <            parser.mdfile();
150 >      antlr::ASTFactory factory;
151 >      parser.initializeASTFactory(factory);
152 >      parser.setASTFactory(&factory);
153 >      parser.mdfile();
154 >      // Create a tree parser that reads information into Globals
155 >      MDTreeParser treeParser;
156 >      treeParser.initializeASTFactory(factory);
157 >      treeParser.setASTFactory(&factory);
158 >      simParams = treeParser.walkTree(parser.getAST());
159 >    }
160  
143            // Create a tree parser that reads information into Globals
144            MDTreeParser treeParser;
145            treeParser.initializeASTFactory(factory);
146            treeParser.setASTFactory(&factory);
147             simParams = treeParser.walkTree(parser.getAST());
148
149        }
150
161        
162 <      catch(antlr::MismatchedCharException& e) {
163 <          sprintf(painCave.errMsg,
164 <                  "parser exception: %s %s:%d:%d\n",
165 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
166 <          painCave.isFatal = 1;
167 <          simError();          
168 <      }
169 <      catch(antlr::MismatchedTokenException &e) {
170 <          sprintf(painCave.errMsg,
171 <                  "parser exception: %s %s:%d:%d\n",
172 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
173 <          painCave.isFatal = 1;
174 <          simError();  
175 <      }
176 <      catch(antlr::NoViableAltForCharException &e) {
177 <          sprintf(painCave.errMsg,
178 <                  "parser exception: %s %s:%d:%d\n",
179 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
180 <          painCave.isFatal = 1;
181 <          simError();  
182 <      }
183 <      catch(antlr::NoViableAltException &e) {
184 <          sprintf(painCave.errMsg,
185 <                  "parser exception: %s %s:%d:%d\n",
186 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
187 <          painCave.isFatal = 1;
188 <          simError();  
189 <      }
162 >    catch(antlr::MismatchedCharException& e) {
163 >      sprintf(painCave.errMsg,
164 >              "parser exception: %s %s:%d:%d\n",
165 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
166 >      painCave.isFatal = 1;
167 >      simError();          
168 >    }
169 >    catch(antlr::MismatchedTokenException &e) {
170 >      sprintf(painCave.errMsg,
171 >              "parser exception: %s %s:%d:%d\n",
172 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
173 >      painCave.isFatal = 1;
174 >      simError();  
175 >    }
176 >    catch(antlr::NoViableAltForCharException &e) {
177 >      sprintf(painCave.errMsg,
178 >              "parser exception: %s %s:%d:%d\n",
179 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
180 >      painCave.isFatal = 1;
181 >      simError();  
182 >    }
183 >    catch(antlr::NoViableAltException &e) {
184 >      sprintf(painCave.errMsg,
185 >              "parser exception: %s %s:%d:%d\n",
186 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
187 >      painCave.isFatal = 1;
188 >      simError();  
189 >    }
190        
191 <        catch(antlr::TokenStreamRecognitionException& e) {
192 <          sprintf(painCave.errMsg,
193 <                  "parser exception: %s %s:%d:%d\n",
194 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
195 <          painCave.isFatal = 1;
196 <          simError();  
197 <        }
191 >    catch(antlr::TokenStreamRecognitionException& e) {
192 >      sprintf(painCave.errMsg,
193 >              "parser exception: %s %s:%d:%d\n",
194 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
195 >      painCave.isFatal = 1;
196 >      simError();  
197 >    }
198          
199 <        catch(antlr::TokenStreamIOException& e) {
200 <          sprintf(painCave.errMsg,
201 <                  "parser exception: %s\n",
202 <                  e.getMessage().c_str());
203 <          painCave.isFatal = 1;
204 <          simError();
205 <        }
199 >    catch(antlr::TokenStreamIOException& e) {
200 >      sprintf(painCave.errMsg,
201 >              "parser exception: %s\n",
202 >              e.getMessage().c_str());
203 >      painCave.isFatal = 1;
204 >      simError();
205 >    }
206          
207 <        catch(antlr::TokenStreamException& e) {
208 <          sprintf(painCave.errMsg,
209 <                  "parser exception: %s\n",
210 <                  e.getMessage().c_str());
211 <          painCave.isFatal = 1;
212 <          simError();
213 <        }        
214 <       catch (antlr::RecognitionException& e) {
215 <          sprintf(painCave.errMsg,
216 <                  "parser exception: %s %s:%d:%d\n",
217 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
218 <          painCave.isFatal = 1;
219 <          simError();          
220 <       }
221 <       catch (antlr::CharStreamException& e) {
222 <          sprintf(painCave.errMsg,
223 <                  "parser exception: %s\n",
224 <                  e.getMessage().c_str());
225 <          painCave.isFatal = 1;
226 <          simError();        
227 <       }
228 <       catch (OOPSEException& e) {
229 <          sprintf(painCave.errMsg,
230 <                  "%s\n",
231 <                  e.getMessage().c_str());
232 <          painCave.isFatal = 1;
233 <          simError();
234 <       }
235 <       catch (std::exception& e) {
236 <          sprintf(painCave.errMsg,
237 <                  "parser exception: %s\n",
238 <                  e.what());
239 <          painCave.isFatal = 1;
240 <          simError();
241 <       }
207 >    catch(antlr::TokenStreamException& e) {
208 >      sprintf(painCave.errMsg,
209 >              "parser exception: %s\n",
210 >              e.getMessage().c_str());
211 >      painCave.isFatal = 1;
212 >      simError();
213 >    }        
214 >    catch (antlr::RecognitionException& e) {
215 >      sprintf(painCave.errMsg,
216 >              "parser exception: %s %s:%d:%d\n",
217 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
218 >      painCave.isFatal = 1;
219 >      simError();          
220 >    }
221 >    catch (antlr::CharStreamException& e) {
222 >      sprintf(painCave.errMsg,
223 >              "parser exception: %s\n",
224 >              e.getMessage().c_str());
225 >      painCave.isFatal = 1;
226 >      simError();        
227 >    }
228 >    catch (OpenMDException& e) {
229 >      sprintf(painCave.errMsg,
230 >              "%s\n",
231 >              e.what());
232 >      painCave.isFatal = 1;
233 >      simError();
234 >    }
235 >    catch (std::exception& e) {
236 >      sprintf(painCave.errMsg,
237 >              "parser exception: %s\n",
238 >              e.what());
239 >      painCave.isFatal = 1;
240 >      simError();
241 >    }
242  
243 <        return simParams;
243 >    simParams->setMDfileVersion(mdFileVersion);
244 >    return simParams;
245    }
246    
247    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
248                                    bool loadInitCoords) {
249 +    
250 +    const int bufferSize = 65535;
251 +    char buffer[bufferSize];
252 +    int lineNo = 0;
253 +    std::string mdRawData;
254 +    int metaDataBlockStart = -1;
255 +    int metaDataBlockEnd = -1;
256 +    streamoff mdOffset;
257 +    int mdFileVersion;
258  
259 +    // Create a string for embedding the version information in the MetaData
260 +    std::string version;
261 +    version.assign("## Last run using OpenMD Version: ");
262 +    version.append(OPENMD_VERSION_MAJOR);
263 +    version.append(".");
264 +    version.append(OPENMD_VERSION_MINOR);
265 +
266 +    std::string svnrev(g_REVISION, strnlen(g_REVISION, 20));
267 +    //convert a macro from compiler to a string in c++
268 +    // STR_DEFINE(svnrev, SVN_REV );
269 +    version.append(" Revision: ");
270 +    // If there's no SVN revision, just call this the RELEASE revision.
271 +    if (!svnrev.empty()) {
272 +      version.append(svnrev);
273 +    } else {
274 +      version.append("RELEASE");
275 +    }
276 +  
277 + #ifdef IS_MPI            
278 +    const int masterNode = 0;
279 +    if (worldRank == masterNode) {
280 + #endif
281 +
282 +      std::ifstream mdFile_;
283 +      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
284 +      
285 +      if (mdFile_.fail()) {
286 +        sprintf(painCave.errMsg,
287 +                "SimCreator: Cannot open file: %s\n",
288 +                mdFileName.c_str());
289 +        painCave.isFatal = 1;
290 +        simError();
291 +      }
292 +
293 +      mdFile_.getline(buffer, bufferSize);
294 +      ++lineNo;
295 +      std::string line = trimLeftCopy(buffer);
296 +      std::size_t i = CaseInsensitiveFind(line, "<OpenMD");
297 +      if (i == string::npos) {
298 +        // try the older file strings to see if that works:
299 +        i = CaseInsensitiveFind(line, "<OOPSE");
300 +      }
301 +      
302 +      if (i == string::npos) {
303 +        // still no luck!
304 +        sprintf(painCave.errMsg,
305 +                "SimCreator: File: %s is not a valid OpenMD file!\n",
306 +                mdFileName.c_str());
307 +        painCave.isFatal = 1;
308 +        simError();
309 +      }
310 +      
311 +      // found the correct opening string, now try to get the file
312 +      // format version number.
313 +
314 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
315 +      std::string fileType = tokenizer.nextToken();
316 +      toUpper(fileType);
317 +
318 +      mdFileVersion = 0;
319 +
320 +      if (fileType == "OPENMD") {
321 +        while (tokenizer.hasMoreTokens()) {
322 +          std::string token(tokenizer.nextToken());
323 +          toUpper(token);
324 +          if (token == "VERSION") {
325 +            mdFileVersion = tokenizer.nextTokenAsInt();
326 +            break;
327 +          }
328 +        }
329 +      }
330 +            
331 +      //scan through the input stream and find MetaData tag        
332 +      while(mdFile_.getline(buffer, bufferSize)) {
333 +        ++lineNo;
334 +        
335 +        std::string line = trimLeftCopy(buffer);
336 +        if (metaDataBlockStart == -1) {
337 +          std::size_t i = CaseInsensitiveFind(line, "<MetaData>");
338 +          if (i != string::npos) {
339 +            metaDataBlockStart = lineNo;
340 +            mdOffset = mdFile_.tellg();
341 +          }
342 +        } else {
343 +          std::size_t i = CaseInsensitiveFind(line, "</MetaData>");
344 +          if (i != string::npos) {
345 +            metaDataBlockEnd = lineNo;
346 +          }
347 +        }
348 +      }
349 +
350 +      if (metaDataBlockStart == -1) {
351 +        sprintf(painCave.errMsg,
352 +                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
353 +                mdFileName.c_str());
354 +        painCave.isFatal = 1;
355 +        simError();
356 +      }
357 +      if (metaDataBlockEnd == -1) {
358 +        sprintf(painCave.errMsg,
359 +                "SimCreator: File: %s did not contain a closed MetaData block!\n",
360 +                mdFileName.c_str());
361 +        painCave.isFatal = 1;
362 +        simError();
363 +      }
364 +        
365 +      mdFile_.clear();
366 +      mdFile_.seekg(0);
367 +      mdFile_.seekg(mdOffset);
368 +
369 +      mdRawData.clear();
370 +
371 +      bool foundVersion = false;
372 +
373 +      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
374 +        mdFile_.getline(buffer, bufferSize);
375 +        std::string line = trimLeftCopy(buffer);
376 +        std::size_t j = CaseInsensitiveFind(line,
377 +                                            "## Last run using OpenMD Version");
378 +        if (j != string::npos) {
379 +          foundVersion = true;
380 +          mdRawData += version;
381 +        } else {
382 +          mdRawData += buffer;
383 +        }
384 +        mdRawData += "\n";
385 +      }
386 +      
387 +      if (!foundVersion) mdRawData += version + "\n";
388 +      
389 +      mdFile_.close();
390 +
391 + #ifdef IS_MPI
392 +    }
393 + #endif
394 +
395 +    std::stringstream rawMetaDataStream(mdRawData);
396 +
397      //parse meta-data file
398 <    Globals* simParams = parseFile(mdFileName);
398 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
399 >                                   metaDataBlockStart + 1);
400      
401      //create the force field
402 <    ForceField * ff = ForceFieldFactory::getInstance()
403 <      ->createForceField(simParams->getForceField());
245 <    
402 >    ForceField * ff = new ForceField(simParams->getForceField());
403 >
404      if (ff == NULL) {
405        sprintf(painCave.errMsg,
406                "ForceField Factory can not create %s force field\n",
# Line 275 | Line 433 | Globals* SimCreator::parseFile(const std::string mdFil
433      }
434      
435      ff->parse(forcefieldFileName);
278    ff->setFortranForceOptions();
436      //create SimInfo
437      SimInfo * info = new SimInfo(ff, simParams);
438 +
439 +    info->setRawMetaData(mdRawData);
440      
441      //gather parameters (SimCreator only retrieves part of the
442      //parameters)
# Line 291 | Line 450 | Globals* SimCreator::parseFile(const std::string mdFil
450      //create the molecules
451      createMolecules(info);
452      
453 <    
453 >    //find the storage layout
454 >
455 >    int storageLayout = computeStorageLayout(info);
456 >
457      //allocate memory for DataStorage(circular reference, need to
458      //break it)
459 <    info->setSnapshotManager(new SimSnapshotManager(info));
459 >    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
460      
461      //set the global index of atoms, rigidbodies and cutoffgroups
462      //(only need to be set once, the global index will never change
# Line 303 | Line 465 | Globals* SimCreator::parseFile(const std::string mdFil
465      //responsibility to LocalIndexManager.
466      setGlobalIndex(info);
467      
468 <    //Although addExcludePairs is called inside SimInfo's addMolecule
468 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
469      //method, at that point atoms don't have the global index yet
470      //(their global index are all initialized to -1).  Therefore we
471 <    //have to call addExcludePairs explicitly here. A way to work
471 >    //have to call addInteractionPairs explicitly here. A way to work
472      //around is that we can determine the beginning global indices of
473      //atoms before they get created.
474      SimInfo::MoleculeIterator mi;
475      Molecule* mol;
476      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
477 <      info->addExcludePairs(mol);
477 >      info->addInteractionPairs(mol);
478      }
479      
480      if (loadInitCoords)
481 <      loadCoordinates(info);    
320 <    
481 >      loadCoordinates(info, mdFileName);    
482      return info;
483    }
484    
# Line 352 | Line 513 | Globals* SimCreator::parseFile(const std::string mdFil
513    
514   #ifdef IS_MPI
515    void SimCreator::divideMolecules(SimInfo *info) {
355    RealType numerator;
356    RealType denominator;
357    RealType precast;
358    RealType x;
359    RealType y;
516      RealType a;
361    int old_atoms;
362    int add_atoms;
363    int new_atoms;
364    int nTarget;
365    int done;
366    int i;
367    int j;
368    int loops;
369    int which_proc;
517      int nProcessors;
518      std::vector<int> atomsPerProc;
519      int nGlobalMols = info->getNGlobalMolecules();
520 <    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
520 >    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an
521 >                                                    // error
522 >                                                    // condition:
523      
524 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
524 >    MPI_Comm_size( MPI_COMM_WORLD, &nProcessors);    
525      
526      if (nProcessors > nGlobalMols) {
527        sprintf(painCave.errMsg,
# Line 381 | Line 530 | Globals* SimCreator::parseFile(const std::string mdFil
530                "\tthe number of molecules.  This will not result in a \n"
531                "\tusable division of atoms for force decomposition.\n"
532                "\tEither try a smaller number of processors, or run the\n"
533 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
533 >              "\tsingle-processor version of OpenMD.\n", nProcessors,
534 >              nGlobalMols);
535        
536        painCave.isFatal = 1;
537        simError();
538      }
539      
390    int seedValue;
540      Globals * simParams = info->getSimParams();
541 <    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
541 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel
542 >                             //random number generator
543      if (simParams->haveSeed()) {
544 <      seedValue = simParams->getSeed();
544 >      int seedValue = simParams->getSeed();
545        myRandom = new SeqRandNumGen(seedValue);
546      }else {
547        myRandom = new SeqRandNumGen();
# Line 404 | Line 554 | Globals* SimCreator::parseFile(const std::string mdFil
554      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
555      
556      if (worldRank == 0) {
557 <      numerator = info->getNGlobalAtoms();
558 <      denominator = nProcessors;
559 <      precast = numerator / denominator;
560 <      nTarget = (int)(precast + 0.5);
557 >      RealType numerator = info->getNGlobalAtoms();
558 >      RealType denominator = nProcessors;
559 >      RealType precast = numerator / denominator;
560 >      int nTarget = (int)(precast + 0.5);
561        
562 <      for(i = 0; i < nGlobalMols; i++) {
563 <        done = 0;
564 <        loops = 0;
562 >      for(int i = 0; i < nGlobalMols; i++) {
563 >
564 >        int done = 0;
565 >        int loops = 0;
566          
567          while (!done) {
568            loops++;
569            
570            // Pick a processor at random
571            
572 <          which_proc = (int) (myRandom->rand() * nProcessors);
572 >          int which_proc = (int) (myRandom->rand() * nProcessors);
573            
574            //get the molecule stamp first
575            int stampId = info->getMoleculeStampId(i);
576            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
577            
578            // How many atoms does this processor have so far?
579 <          old_atoms = atomsPerProc[which_proc];
580 <          add_atoms = moleculeStamp->getNAtoms();
581 <          new_atoms = old_atoms + add_atoms;
579 >          int old_atoms = atomsPerProc[which_proc];
580 >          int add_atoms = moleculeStamp->getNAtoms();
581 >          int new_atoms = old_atoms + add_atoms;
582            
583            // If we've been through this loop too many times, we need
584            // to just give up and assign the molecule to this processor
585            // and be done with it.
586            
587            if (loops > 100) {
588 +
589              sprintf(painCave.errMsg,
590 <                    "I've tried 100 times to assign molecule %d to a "
591 <                    " processor, but can't find a good spot.\n"
592 <                    "I'm assigning it at random to processor %d.\n",
590 >                    "There have been 100 attempts to assign molecule %d to an\n"
591 >                    "\tunderworked processor, but there's no good place to\n"
592 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
593                      i, which_proc);
594 <            
594 >          
595              painCave.isFatal = 0;
596 +            painCave.severity = OPENMD_INFO;
597              simError();
598              
599              molToProcMap[i] = which_proc;
# Line 469 | Line 622 | Globals* SimCreator::parseFile(const std::string mdFil
622            //           Pacc(x) = exp(- a * x)
623            // where a = penalty / (average atoms per molecule)
624            
625 <          x = (RealType)(new_atoms - nTarget);
626 <          y = myRandom->rand();
625 >          RealType x = (RealType)(new_atoms - nTarget);
626 >          RealType y = myRandom->rand();
627            
628            if (y < exp(- a * x)) {
629              molToProcMap[i] = which_proc;
# Line 485 | Line 638 | Globals* SimCreator::parseFile(const std::string mdFil
638        }
639        
640        delete myRandom;
641 <      
641 >
642        // Spray out this nonsense to all other processors:
490      
643        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
644 +
645      } else {
646        
647        // Listen to your marching orders from processor 0:
495      
648        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
649 +
650      }
651      
652      info->setMolToProcMap(molToProcMap);
653      sprintf(checkPointMsg,
654              "Successfully divided the molecules among the processors.\n");
655 <    MPIcheckPoint();
655 >    errorCheckPoint();
656    }
657    
658   #endif
# Line 516 | Line 669 | Globals* SimCreator::parseFile(const std::string mdFil
669   #endif
670          
671          stampId = info->getMoleculeStampId(i);
672 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
673 <                                                   stampId, i, info->getLocalIndexManager());
672 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
673 >                                                   info->getMoleculeStamp(stampId),
674 >                                                   stampId, i,
675 >                                                   info->getLocalIndexManager());
676          
677          info->addMolecule(mol);
678          
# Line 529 | Line 684 | Globals* SimCreator::parseFile(const std::string mdFil
684        
685      } //end for(int i=0)  
686    }
687 +    
688 +  int SimCreator::computeStorageLayout(SimInfo* info) {
689 +
690 +    Globals* simParams = info->getSimParams();
691 +    int nRigidBodies = info->getNGlobalRigidBodies();
692 +    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
693 +    set<AtomType*>::iterator i;
694 +    bool hasDirectionalAtoms = false;
695 +    bool hasFixedCharge = false;
696 +    bool hasDipoles = false;    
697 +    bool hasQuadrupoles = false;    
698 +    bool hasPolarizable = false;    
699 +    bool hasFluctuatingCharge = false;    
700 +    bool hasMetallic = false;
701 +    int storageLayout = 0;
702 +    storageLayout |= DataStorage::dslPosition;
703 +    storageLayout |= DataStorage::dslVelocity;
704 +    storageLayout |= DataStorage::dslForce;
705 +
706 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
707 +
708 +      DirectionalAdapter da = DirectionalAdapter( (*i) );
709 +      MultipoleAdapter ma = MultipoleAdapter( (*i) );
710 +      EAMAdapter ea = EAMAdapter( (*i) );
711 +      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
712 +      PolarizableAdapter pa = PolarizableAdapter( (*i) );
713 +      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
714 +      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
715 +
716 +      if (da.isDirectional()){
717 +        hasDirectionalAtoms = true;
718 +      }
719 +      if (ma.isDipole()){
720 +        hasDipoles = true;
721 +      }
722 +      if (ma.isQuadrupole()){
723 +        hasQuadrupoles = true;
724 +      }
725 +      if (ea.isEAM() || sca.isSuttonChen()){
726 +        hasMetallic = true;
727 +      }
728 +      if ( fca.isFixedCharge() ){
729 +        hasFixedCharge = true;
730 +      }
731 +      if ( fqa.isFluctuatingCharge() ){
732 +        hasFluctuatingCharge = true;
733 +      }
734 +      if ( pa.isPolarizable() ){
735 +        hasPolarizable = true;
736 +      }
737 +    }
738      
739 +    if (nRigidBodies > 0 || hasDirectionalAtoms) {
740 +      storageLayout |= DataStorage::dslAmat;
741 +      if(storageLayout & DataStorage::dslVelocity) {
742 +        storageLayout |= DataStorage::dslAngularMomentum;
743 +      }
744 +      if (storageLayout & DataStorage::dslForce) {
745 +        storageLayout |= DataStorage::dslTorque;
746 +      }
747 +    }
748 +    if (hasDipoles) {
749 +      storageLayout |= DataStorage::dslDipole;
750 +    }
751 +    if (hasQuadrupoles) {
752 +      storageLayout |= DataStorage::dslQuadrupole;
753 +    }
754 +    if (hasFixedCharge || hasFluctuatingCharge) {
755 +      storageLayout |= DataStorage::dslSkippedCharge;
756 +    }
757 +    if (hasMetallic) {
758 +      storageLayout |= DataStorage::dslDensity;
759 +      storageLayout |= DataStorage::dslFunctional;
760 +      storageLayout |= DataStorage::dslFunctionalDerivative;
761 +    }
762 +    if (hasPolarizable) {
763 +      storageLayout |= DataStorage::dslElectricField;
764 +    }
765 +    if (hasFluctuatingCharge){
766 +      storageLayout |= DataStorage::dslFlucQPosition;
767 +      if(storageLayout & DataStorage::dslVelocity) {
768 +        storageLayout |= DataStorage::dslFlucQVelocity;
769 +      }
770 +      if (storageLayout & DataStorage::dslForce) {
771 +        storageLayout |= DataStorage::dslFlucQForce;
772 +      }
773 +    }
774 +    
775 +    // if the user has asked for them, make sure we've got the memory for the
776 +    // objects defined.
777 +
778 +    if (simParams->getOutputParticlePotential()) {
779 +      storageLayout |= DataStorage::dslParticlePot;
780 +    }
781 +
782 +    if (simParams->havePrintHeatFlux()) {
783 +      if (simParams->getPrintHeatFlux()) {
784 +        storageLayout |= DataStorage::dslParticlePot;
785 +      }
786 +    }
787 +
788 +    if (simParams->getOutputElectricField() |
789 +        simParams->haveElectricField() | simParams->haveUniformField() |
790 +        simParams->haveUniformGradientStrength() |
791 +        simParams->haveUniformGradientDirection1() |
792 +        simParams->haveUniformGradientDirection2() ) {
793 +      storageLayout |= DataStorage::dslElectricField;
794 +    }
795 +
796 +    if (simParams->getOutputSitePotential() ) {
797 +      storageLayout |= DataStorage::dslSitePotential;
798 +    }
799 +
800 +    if (simParams->getOutputFluctuatingCharges()) {
801 +      storageLayout |= DataStorage::dslFlucQPosition;
802 +      storageLayout |= DataStorage::dslFlucQVelocity;
803 +      storageLayout |= DataStorage::dslFlucQForce;
804 +    }
805 +
806 +    info->setStorageLayout(storageLayout);
807 +
808 +    return storageLayout;
809 +  }
810 +
811    void SimCreator::setGlobalIndex(SimInfo *info) {
812      SimInfo::MoleculeIterator mi;
813      Molecule::AtomIterator ai;
814      Molecule::RigidBodyIterator ri;
815      Molecule::CutoffGroupIterator ci;
816 <    Molecule * mol;
817 <    Atom * atom;
818 <    RigidBody * rb;
819 <    CutoffGroup * cg;
816 >    Molecule::BondIterator boi;
817 >    Molecule::BendIterator bei;
818 >    Molecule::TorsionIterator ti;
819 >    Molecule::InversionIterator ii;
820 >    Molecule::IntegrableObjectIterator  ioi;
821 >    Molecule* mol;
822 >    Atom* atom;
823 >    RigidBody* rb;
824 >    CutoffGroup* cg;
825 >    Bond* bond;
826 >    Bend* bend;
827 >    Torsion* torsion;
828 >    Inversion* inversion;
829      int beginAtomIndex;
830      int beginRigidBodyIndex;
831      int beginCutoffGroupIndex;
832 +    int beginBondIndex;
833 +    int beginBendIndex;
834 +    int beginTorsionIndex;
835 +    int beginInversionIndex;
836 + #ifdef IS_MPI
837      int nGlobalAtoms = info->getNGlobalAtoms();
838 <    
839 < #ifndef IS_MPI
838 >    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
839 > #endif
840      
841      beginAtomIndex = 0;
842 <    beginRigidBodyIndex = 0;
842 >    // The rigid body indices begin immediately after the atom indices:
843 >    beginRigidBodyIndex = info->getNGlobalAtoms();
844      beginCutoffGroupIndex = 0;
845 <    
846 < #else
847 <    
848 <    int nproc;
849 <    int myNode;
850 <    
558 <    myNode = worldRank;
559 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
560 <    
561 <    std::vector < int > tmpAtomsInProc(nproc, 0);
562 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
563 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
564 <    std::vector < int > NumAtomsInProc(nproc, 0);
565 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
566 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
567 <    
568 <    tmpAtomsInProc[myNode] = info->getNAtoms();
569 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
570 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
571 <    
572 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
573 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
574 <                  MPI_SUM, MPI_COMM_WORLD);
575 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
576 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
577 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
578 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
579 <    
580 <    beginAtomIndex = 0;
581 <    beginRigidBodyIndex = 0;
582 <    beginCutoffGroupIndex = 0;
583 <    
584 <    for(int i = 0; i < myNode; i++) {
585 <      beginAtomIndex += NumAtomsInProc[i];
586 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
587 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
588 <    }
589 <    
590 < #endif
591 <    
592 <    //rigidbody's index begins right after atom's
593 <    beginRigidBodyIndex += info->getNGlobalAtoms();
594 <    
595 <    for(mol = info->beginMolecule(mi); mol != NULL;
596 <        mol = info->nextMolecule(mi)) {
845 >    beginBondIndex = 0;
846 >    beginBendIndex = 0;
847 >    beginTorsionIndex = 0;
848 >    beginInversionIndex = 0;
849 >  
850 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
851        
852 <      //local index(index in DataStorge) of atom is important
853 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
854 <        atom->setGlobalIndex(beginAtomIndex++);
852 > #ifdef IS_MPI      
853 >      if (info->getMolToProc(i) == worldRank) {
854 > #endif        
855 >        // stuff to do if I own this molecule
856 >        mol = info->getMoleculeByGlobalIndex(i);
857 >
858 >        // The local index(index in DataStorge) of the atom is important:
859 >        for(atom = mol->beginAtom(ai); atom != NULL;
860 >            atom = mol->nextAtom(ai)) {
861 >          atom->setGlobalIndex(beginAtomIndex++);
862 >        }
863 >        
864 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
865 >            rb = mol->nextRigidBody(ri)) {
866 >          rb->setGlobalIndex(beginRigidBodyIndex++);
867 >        }
868 >        
869 >        // The local index of other objects only depends on the order
870 >        // of traversal:
871 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
872 >            cg = mol->nextCutoffGroup(ci)) {
873 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
874 >        }        
875 >        for(bond = mol->beginBond(boi); bond != NULL;
876 >            bond = mol->nextBond(boi)) {
877 >          bond->setGlobalIndex(beginBondIndex++);
878 >        }        
879 >        for(bend = mol->beginBend(bei); bend != NULL;
880 >            bend = mol->nextBend(bei)) {
881 >          bend->setGlobalIndex(beginBendIndex++);
882 >        }        
883 >        for(torsion = mol->beginTorsion(ti); torsion != NULL;
884 >            torsion = mol->nextTorsion(ti)) {
885 >          torsion->setGlobalIndex(beginTorsionIndex++);
886 >        }        
887 >        for(inversion = mol->beginInversion(ii); inversion != NULL;
888 >            inversion = mol->nextInversion(ii)) {
889 >          inversion->setGlobalIndex(beginInversionIndex++);
890 >        }        
891 >        
892 > #ifdef IS_MPI        
893 >      }  else {
894 >
895 >        // stuff to do if I don't own this molecule
896 >        
897 >        int stampId = info->getMoleculeStampId(i);
898 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
899 >
900 >        beginAtomIndex += stamp->getNAtoms();
901 >        beginRigidBodyIndex += stamp->getNRigidBodies();
902 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
903 >        beginBondIndex += stamp->getNBonds();
904 >        beginBendIndex += stamp->getNBends();
905 >        beginTorsionIndex += stamp->getNTorsions();
906 >        beginInversionIndex += stamp->getNInversions();
907        }
908 <      
909 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
910 <          rb = mol->nextRigidBody(ri)) {
911 <        rb->setGlobalIndex(beginRigidBodyIndex++);
606 <      }
607 <      
608 <      //local index of cutoff group is trivial, it only depends on the order of travesing
609 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
610 <          cg = mol->nextCutoffGroup(ci)) {
611 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
612 <      }
613 <    }
614 <    
908 > #endif          
909 >
910 >    } //end for(int i=0)  
911 >
912      //fill globalGroupMembership
913      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
914 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
915 <      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
916 <        
914 >    for(mol = info->beginMolecule(mi); mol != NULL;
915 >        mol = info->nextMolecule(mi)) {        
916 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
917 >           cg = mol->nextCutoffGroup(ci)) {        
918          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
919            globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
920          }
921          
922        }      
923      }
924 <    
924 >  
925   #ifdef IS_MPI    
926      // Since the globalGroupMembership has been zero filled and we've only
927      // poked values into the atoms we know, we can do an Allreduce
928      // to get the full globalGroupMembership array (We think).
929      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
930      // docs said we could.
931 <    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
932 <    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
931 >    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
932 >    MPI_Allreduce(&globalGroupMembership[0],
933 >                  &tmpGroupMembership[0], nGlobalAtoms,
934                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
935 +
936      info->setGlobalGroupMembership(tmpGroupMembership);
937   #else
938      info->setGlobalGroupMembership(globalGroupMembership);
939   #endif
940      
941      //fill molMembership
942 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
942 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
943 >                                         info->getNGlobalRigidBodies(), 0);
944      
945 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
946 <      
945 >    for(mol = info->beginMolecule(mi); mol != NULL;
946 >        mol = info->nextMolecule(mi)) {
947        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
948          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
949        }
950 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
951 +           rb = mol->nextRigidBody(ri)) {
952 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
953 +      }
954      }
955      
956   #ifdef IS_MPI
957 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
958 <    
959 <    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
957 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
958 >                                      info->getNGlobalRigidBodies(), 0);
959 >    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
960 >                  nGlobalAtoms + nGlobalRigidBodies,
961                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
962      
963      info->setGlobalMolMembership(tmpMolMembership);
964   #else
965      info->setGlobalMolMembership(globalMolMembership);
966   #endif
967 +
968 +    // nIOPerMol holds the number of integrable objects per molecule
969 +    // here the molecules are listed by their global indices.
970 +
971 +    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
972 +    for (mol = info->beginMolecule(mi); mol != NULL;
973 +         mol = info->nextMolecule(mi)) {
974 +      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
975 +    }
976      
977 <  }
978 <  
979 <  void SimCreator::loadCoordinates(SimInfo* info) {
980 <    Globals* simParams;
981 <    simParams = info->getSimParams();
977 > #ifdef IS_MPI
978 >    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
979 >    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
980 >      info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
981 > #else
982 >    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
983 > #endif    
984 >
985 >    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
986      
987 <    if (!simParams->haveInitialConfig()) {
988 <      sprintf(painCave.errMsg,
989 <              "Cannot intialize a simulation without an initial configuration file.\n");
990 <      painCave.isFatal = 1;;
672 <      simError();
987 >    int startingIndex = 0;
988 >    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
989 >      startingIOIndexForMol[i] = startingIndex;
990 >      startingIndex += numIntegrableObjectsPerMol[i];
991      }
992      
993 <    DumpReader reader(info, simParams->getInitialConfig());
993 >    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
994 >    for (mol = info->beginMolecule(mi); mol != NULL;
995 >         mol = info->nextMolecule(mi)) {
996 >      int myGlobalIndex = mol->getGlobalIndex();
997 >      int globalIO = startingIOIndexForMol[myGlobalIndex];
998 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
999 >           sd = mol->nextIntegrableObject(ioi)) {
1000 >        sd->setGlobalIntegrableObjectIndex(globalIO);
1001 >        IOIndexToIntegrableObject[globalIO] = sd;
1002 >        globalIO++;
1003 >      }
1004 >    }
1005 >      
1006 >    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
1007 >    
1008 >  }
1009 >  
1010 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
1011 >    
1012 >    DumpReader reader(info, mdFileName);
1013      int nframes = reader.getNFrames();
1014      
1015      if (nframes > 0) {
# Line 681 | Line 1018 | Globals* SimCreator::parseFile(const std::string mdFil
1018        //invalid initial coordinate file
1019        sprintf(painCave.errMsg,
1020                "Initial configuration file %s should at least contain one frame\n",
1021 <              simParams->getInitialConfig().c_str());
1021 >              mdFileName.c_str());
1022        painCave.isFatal = 1;
1023        simError();
1024      }
688    
1025      //copy the current snapshot to previous snapshot
1026      info->getSnapshotManager()->advance();
1027    }
1028    
1029 < } //end namespace oopse
1029 > } //end namespace OpenMD
1030  
1031  

Comparing trunk/src/brains/SimCreator.cpp (property svn:keywords):
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 2069 by gezelter, Thu Mar 5 16:30:23 2015 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines