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 749 by tim, Wed Nov 16 23:10:02 2005 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 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>
58 + #include <string>
59 +
60   #include "brains/MoleculeCreator.hpp"
61   #include "brains/SimCreator.hpp"
62   #include "brains/SimSnapshotManager.hpp"
63   #include "io/DumpReader.hpp"
64 < #include "io/parse_me.h"
55 < #include "UseTheForce/ForceFieldFactory.hpp"
64 > #include "brains/ForceField.hpp"
65   #include "utils/simError.h"
66   #include "utils/StringUtils.hpp"
67   #include "math/SeqRandNumGen.hpp"
68 < #ifdef IS_MPI
69 < #include "io/mpiBASS.h"
70 < #include "math/ParallelRandNumGen.hpp"
71 < #endif
68 > #include "mdParser/MDLexer.hpp"
69 > #include "mdParser/MDParser.hpp"
70 > #include "mdParser/MDTreeParser.hpp"
71 > #include "mdParser/SimplePreprocessor.hpp"
72 > #include "antlr/ANTLRException.hpp"
73 > #include "antlr/TokenStreamRecognitionException.hpp"
74 > #include "antlr/TokenStreamIOException.hpp"
75 > #include "antlr/TokenStreamException.hpp"
76 > #include "antlr/RecognitionException.hpp"
77 > #include "antlr/CharStreamException.hpp"
78  
79 < namespace oopse {
79 > #include "antlr/MismatchedCharException.hpp"
80 > #include "antlr/MismatchedTokenException.hpp"
81 > #include "antlr/NoViableAltForCharException.hpp"
82 > #include "antlr/NoViableAltException.hpp"
83 >
84 > #include "types/DirectionalAdapter.hpp"
85 > #include "types/MultipoleAdapter.hpp"
86 > #include "types/EAMAdapter.hpp"
87 > #include "types/SuttonChenAdapter.hpp"
88 > #include "types/PolarizableAdapter.hpp"
89 > #include "types/FixedChargeAdapter.hpp"
90 > #include "types/FluctuatingChargeAdapter.hpp"
91 >
92 >
93 > namespace OpenMD {
94    
95 <  void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps,
96 <                             Globals* simParams){
95 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
96 >    Globals* simParams = NULL;
97 >    try {
98 >
99 >      // Create a preprocessor that preprocesses md file into an ostringstream
100 >      std::stringstream ppStream;
101 > #ifdef IS_MPI            
102 >      int streamSize;
103 >      const int masterNode = 0;
104 >
105 >      if (worldRank == masterNode) {
106 >        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
107 >        // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
108 > #endif                
109 >        SimplePreprocessor preprocessor;
110 >        preprocessor.preprocess(rawMetaDataStream, filename,
111 >                                startOfMetaDataBlock, ppStream);
112 >                
113 > #ifdef IS_MPI            
114 >        //broadcasting the stream size
115 >        streamSize = ppStream.str().size() +1;
116 >        MPI_Bcast(&streamSize, 1, MPI_LONG, 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 >
120 >        // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
121 >        // MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())),
122 >        //                       streamSize, MPI::CHAR, masterNode);
123 >                          
124 >      } else {
125 >
126 >        MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
127 >        // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode);
128 >
129 >        //get stream size
130 >        MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);
131 >        // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode);
132 >        char* buf = new char[streamSize];
133 >        assert(buf);
134 >                
135 >        //receive file content
136 >        MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
137 >        // MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode);
138 >
139 >        ppStream.str(buf);
140 >        delete [] buf;
141 >      }
142 > #endif            
143 >      // Create a scanner that reads from the input stream
144 >      MDLexer lexer(ppStream);
145 >      lexer.setFilename(filename);
146 >      lexer.initDeferredLineCount();
147      
148 < #ifdef IS_MPI
148 >      // Create a parser that reads from the scanner
149 >      MDParser parser(lexer);
150 >      parser.setFilename(filename);
151 >
152 >      // Create an observer that synchorizes file name change
153 >      FilenameObserver observer;
154 >      observer.setLexer(&lexer);
155 >      observer.setParser(&parser);
156 >      lexer.setObserver(&observer);
157      
158 <    if (worldRank == 0) {
159 < #endif // is_mpi
158 >      antlr::ASTFactory factory;
159 >      parser.initializeASTFactory(factory);
160 >      parser.setASTFactory(&factory);
161 >      parser.mdfile();
162 >      // Create a tree parser that reads information into Globals
163 >      MDTreeParser treeParser;
164 >      treeParser.initializeASTFactory(factory);
165 >      treeParser.setASTFactory(&factory);
166 >      simParams = treeParser.walkTree(parser.getAST());
167 >    }
168 >
169        
170 <      set_interface_stamps(stamps, simParams);
170 >    catch(antlr::MismatchedCharException& e) {
171 >      sprintf(painCave.errMsg,
172 >              "parser exception: %s %s:%d:%d\n",
173 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
174 >      painCave.isFatal = 1;
175 >      simError();          
176 >    }
177 >    catch(antlr::MismatchedTokenException &e) {
178 >      sprintf(painCave.errMsg,
179 >              "parser exception: %s %s:%d:%d\n",
180 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
181 >      painCave.isFatal = 1;
182 >      simError();  
183 >    }
184 >    catch(antlr::NoViableAltForCharException &e) {
185 >      sprintf(painCave.errMsg,
186 >              "parser exception: %s %s:%d:%d\n",
187 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
188 >      painCave.isFatal = 1;
189 >      simError();  
190 >    }
191 >    catch(antlr::NoViableAltException &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 < #ifdef IS_MPI
200 <      
201 <      mpiEventInit();
202 <      
203 < #endif
204 <      
82 <      yacc_BASS(mdFileName.c_str());
83 <      
84 < #ifdef IS_MPI
85 <      
86 <      throwMPIEvent(NULL);
87 <    } else {
88 <      set_interface_stamps(stamps, simParams);
89 <      mpiEventInit();
90 <      MPIcheckPoint();
91 <      mpiEventLoop();
199 >    catch(antlr::TokenStreamRecognitionException& e) {
200 >      sprintf(painCave.errMsg,
201 >              "parser exception: %s %s:%d:%d\n",
202 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
203 >      painCave.isFatal = 1;
204 >      simError();  
205      }
206 <    
207 < #endif
208 <    
206 >        
207 >    catch(antlr::TokenStreamIOException& e) {
208 >      sprintf(painCave.errMsg,
209 >              "parser exception: %s\n",
210 >              e.getMessage().c_str());
211 >      painCave.isFatal = 1;
212 >      simError();
213 >    }
214 >        
215 >    catch(antlr::TokenStreamException& e) {
216 >      sprintf(painCave.errMsg,
217 >              "parser exception: %s\n",
218 >              e.getMessage().c_str());
219 >      painCave.isFatal = 1;
220 >      simError();
221 >    }        
222 >    catch (antlr::RecognitionException& e) {
223 >      sprintf(painCave.errMsg,
224 >              "parser exception: %s %s:%d:%d\n",
225 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
226 >      painCave.isFatal = 1;
227 >      simError();          
228 >    }
229 >    catch (antlr::CharStreamException& e) {
230 >      sprintf(painCave.errMsg,
231 >              "parser exception: %s\n",
232 >              e.getMessage().c_str());
233 >      painCave.isFatal = 1;
234 >      simError();        
235 >    }
236 >    catch (OpenMDException& e) {
237 >      sprintf(painCave.errMsg,
238 >              "%s\n",
239 >              e.getMessage().c_str());
240 >      painCave.isFatal = 1;
241 >      simError();
242 >    }
243 >    catch (std::exception& e) {
244 >      sprintf(painCave.errMsg,
245 >              "parser exception: %s\n",
246 >              e.what());
247 >      painCave.isFatal = 1;
248 >      simError();
249 >    }
250 >
251 >    simParams->setMDfileVersion(mdFileVersion);
252 >    return simParams;
253    }
254    
255    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
256                                    bool loadInitCoords) {
257      
258 <    MakeStamps * stamps = new MakeStamps();
259 <    
260 <    Globals * simParams = new Globals();
261 <    
258 >    const int bufferSize = 65535;
259 >    char buffer[bufferSize];
260 >    int lineNo = 0;
261 >    std::string mdRawData;
262 >    int metaDataBlockStart = -1;
263 >    int metaDataBlockEnd = -1;
264 >    int i, j;
265 >    streamoff mdOffset;
266 >    int mdFileVersion;
267 >
268 >    // Create a string for embedding the version information in the MetaData
269 >    std::string version;
270 >    version.assign("## Last run using OpenMD Version: ");
271 >    version.append(OPENMD_VERSION_MAJOR);
272 >    version.append(".");
273 >    version.append(OPENMD_VERSION_MINOR);
274 >
275 >    std::string svnrev;
276 >    //convert a macro from compiler to a string in c++
277 >    STR_DEFINE(svnrev, SVN_REV );
278 >    version.append(" Revision: ");
279 >    // If there's no SVN revision, just call this the RELEASE revision.
280 >    if (!svnrev.empty()) {
281 >      version.append(svnrev);
282 >    } else {
283 >      version.append("RELEASE");
284 >    }
285 >  
286 > #ifdef IS_MPI            
287 >    const int masterNode = 0;
288 >    if (worldRank == masterNode) {
289 > #endif
290 >
291 >      std::ifstream mdFile_;
292 >      mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary);
293 >      
294 >      if (mdFile_.fail()) {
295 >        sprintf(painCave.errMsg,
296 >                "SimCreator: Cannot open file: %s\n",
297 >                mdFileName.c_str());
298 >        painCave.isFatal = 1;
299 >        simError();
300 >      }
301 >
302 >      mdFile_.getline(buffer, bufferSize);
303 >      ++lineNo;
304 >      std::string line = trimLeftCopy(buffer);
305 >      i = CaseInsensitiveFind(line, "<OpenMD");
306 >      if (static_cast<size_t>(i) == string::npos) {
307 >        // try the older file strings to see if that works:
308 >        i = CaseInsensitiveFind(line, "<OOPSE");
309 >      }
310 >      
311 >      if (static_cast<size_t>(i) == string::npos) {
312 >        // still no luck!
313 >        sprintf(painCave.errMsg,
314 >                "SimCreator: File: %s is not a valid OpenMD file!\n",
315 >                mdFileName.c_str());
316 >        painCave.isFatal = 1;
317 >        simError();
318 >      }
319 >      
320 >      // found the correct opening string, now try to get the file
321 >      // format version number.
322 >
323 >      StringTokenizer tokenizer(line, "=<> \t\n\r");
324 >      std::string fileType = tokenizer.nextToken();
325 >      toUpper(fileType);
326 >
327 >      mdFileVersion = 0;
328 >
329 >      if (fileType == "OPENMD") {
330 >        while (tokenizer.hasMoreTokens()) {
331 >          std::string token(tokenizer.nextToken());
332 >          toUpper(token);
333 >          if (token == "VERSION") {
334 >            mdFileVersion = tokenizer.nextTokenAsInt();
335 >            break;
336 >          }
337 >        }
338 >      }
339 >            
340 >      //scan through the input stream and find MetaData tag        
341 >      while(mdFile_.getline(buffer, bufferSize)) {
342 >        ++lineNo;
343 >        
344 >        std::string line = trimLeftCopy(buffer);
345 >        if (metaDataBlockStart == -1) {
346 >          i = CaseInsensitiveFind(line, "<MetaData>");
347 >          if (i != string::npos) {
348 >            metaDataBlockStart = lineNo;
349 >            mdOffset = mdFile_.tellg();
350 >          }
351 >        } else {
352 >          i = CaseInsensitiveFind(line, "</MetaData>");
353 >          if (i != string::npos) {
354 >            metaDataBlockEnd = lineNo;
355 >          }
356 >        }
357 >      }
358 >
359 >      if (metaDataBlockStart == -1) {
360 >        sprintf(painCave.errMsg,
361 >                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
362 >                mdFileName.c_str());
363 >        painCave.isFatal = 1;
364 >        simError();
365 >      }
366 >      if (metaDataBlockEnd == -1) {
367 >        sprintf(painCave.errMsg,
368 >                "SimCreator: File: %s did not contain a closed MetaData block!\n",
369 >                mdFileName.c_str());
370 >        painCave.isFatal = 1;
371 >        simError();
372 >      }
373 >        
374 >      mdFile_.clear();
375 >      mdFile_.seekg(0);
376 >      mdFile_.seekg(mdOffset);
377 >
378 >      mdRawData.clear();
379 >
380 >      bool foundVersion = false;
381 >
382 >      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
383 >        mdFile_.getline(buffer, bufferSize);
384 >        std::string line = trimLeftCopy(buffer);
385 >        j = CaseInsensitiveFind(line, "## Last run using OpenMD Version");
386 >        if (static_cast<size_t>(j) != string::npos) {
387 >          foundVersion = true;
388 >          mdRawData += version;
389 >        } else {
390 >          mdRawData += buffer;
391 >        }
392 >        mdRawData += "\n";
393 >      }
394 >      
395 >      if (!foundVersion) mdRawData += version + "\n";
396 >      
397 >      mdFile_.close();
398 >
399 > #ifdef IS_MPI
400 >    }
401 > #endif
402 >
403 >    std::stringstream rawMetaDataStream(mdRawData);
404 >
405      //parse meta-data file
406 <    parseFile(mdFileName, stamps, simParams);
406 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
407 >                                   metaDataBlockStart + 1);
408      
409      //create the force field
410 <    ForceField * ff = ForceFieldFactory::getInstance()
411 <      ->createForceField(simParams->getForceField());
111 <    
410 >    ForceField * ff = new ForceField(simParams->getForceField());
411 >
412      if (ff == NULL) {
413        sprintf(painCave.errMsg,
414                "ForceField Factory can not create %s force field\n",
# Line 141 | Line 441 | namespace oopse {
441      }
442      
443      ff->parse(forcefieldFileName);
144    
145    //extract the molecule stamps
146    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
147    compList(stamps, simParams, moleculeStampPairs);
148    
444      //create SimInfo
445 <    SimInfo * info = new SimInfo(stamps, moleculeStampPairs, ff, simParams);
445 >    SimInfo * info = new SimInfo(ff, simParams);
446 >
447 >    info->setRawMetaData(mdRawData);
448      
449 <    //gather parameters (SimCreator only retrieves part of the parameters)
449 >    //gather parameters (SimCreator only retrieves part of the
450 >    //parameters)
451      gatherParameters(info, mdFileName);
452      
453      //divide the molecules and determine the global index of molecules
# Line 160 | Line 458 | namespace oopse {
458      //create the molecules
459      createMolecules(info);
460      
461 +    //find the storage layout
462 +
463 +    int storageLayout = computeStorageLayout(info);
464 +
465 +    //allocate memory for DataStorage(circular reference, need to
466 +    //break it)
467 +    info->setSnapshotManager(new SimSnapshotManager(info, storageLayout));
468      
469 <    //allocate memory for DataStorage(circular reference, need to break it)
470 <    info->setSnapshotManager(new SimSnapshotManager(info));
471 <    
472 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
473 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
169 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
469 >    //set the global index of atoms, rigidbodies and cutoffgroups
470 >    //(only need to be set once, the global index will never change
471 >    //again). Local indices of atoms and rigidbodies are already set
472 >    //by MoleculeCreator class which actually delegates the
473 >    //responsibility to LocalIndexManager.
474      setGlobalIndex(info);
475      
476 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
477 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
478 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
479 <    //we can determine the beginning global indices of atoms before they get created.
476 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
477 >    //method, at that point atoms don't have the global index yet
478 >    //(their global index are all initialized to -1).  Therefore we
479 >    //have to call addInteractionPairs explicitly here. A way to work
480 >    //around is that we can determine the beginning global indices of
481 >    //atoms before they get created.
482      SimInfo::MoleculeIterator mi;
483      Molecule* mol;
484      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
485 <      info->addExcludePairs(mol);
485 >      info->addInteractionPairs(mol);
486      }
487      
488      if (loadInitCoords)
489 <      loadCoordinates(info);    
184 <    
489 >      loadCoordinates(info, mdFileName);    
490      return info;
491    }
492    
# Line 216 | Line 521 | namespace oopse {
521    
522   #ifdef IS_MPI
523    void SimCreator::divideMolecules(SimInfo *info) {
524 <    double numerator;
220 <    double denominator;
221 <    double precast;
222 <    double x;
223 <    double y;
224 <    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;
524 >    RealType a;
525      int nProcessors;
526      std::vector<int> atomsPerProc;
527      int nGlobalMols = info->getNGlobalMolecules();
528 <    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
528 >    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an
529 >                                                    // error
530 >                                                    // condition:
531      
532 <    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
532 >    MPI_Comm_size( MPI_COMM_WORLD, &nProcessors);    
533 >    //nProcessors = MPI::COMM_WORLD.Get_size();
534      
535      if (nProcessors > nGlobalMols) {
536        sprintf(painCave.errMsg,
# Line 245 | Line 539 | namespace oopse {
539                "\tthe number of molecules.  This will not result in a \n"
540                "\tusable division of atoms for force decomposition.\n"
541                "\tEither try a smaller number of processors, or run the\n"
542 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
542 >              "\tsingle-processor version of OpenMD.\n", nProcessors,
543 >              nGlobalMols);
544        
545        painCave.isFatal = 1;
546        simError();
547      }
548      
254    int seedValue;
549      Globals * simParams = info->getSimParams();
550 <    SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator
550 >    SeqRandNumGen* myRandom; //divide labor does not need Parallel
551 >                             //random number generator
552      if (simParams->haveSeed()) {
553 <      seedValue = simParams->getSeed();
553 >      int seedValue = simParams->getSeed();
554        myRandom = new SeqRandNumGen(seedValue);
555      }else {
556        myRandom = new SeqRandNumGen();
# Line 268 | Line 563 | namespace oopse {
563      atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
564      
565      if (worldRank == 0) {
566 <      numerator = info->getNGlobalAtoms();
567 <      denominator = nProcessors;
568 <      precast = numerator / denominator;
569 <      nTarget = (int)(precast + 0.5);
566 >      RealType numerator = info->getNGlobalAtoms();
567 >      RealType denominator = nProcessors;
568 >      RealType precast = numerator / denominator;
569 >      int nTarget = (int)(precast + 0.5);
570        
571 <      for(i = 0; i < nGlobalMols; i++) {
572 <        done = 0;
573 <        loops = 0;
571 >      for(int i = 0; i < nGlobalMols; i++) {
572 >
573 >        int done = 0;
574 >        int loops = 0;
575          
576          while (!done) {
577            loops++;
578            
579            // Pick a processor at random
580            
581 <          which_proc = (int) (myRandom->rand() * nProcessors);
581 >          int which_proc = (int) (myRandom->rand() * nProcessors);
582            
583            //get the molecule stamp first
584            int stampId = info->getMoleculeStampId(i);
585            MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
586            
587            // How many atoms does this processor have so far?
588 <          old_atoms = atomsPerProc[which_proc];
589 <          add_atoms = moleculeStamp->getNAtoms();
590 <          new_atoms = old_atoms + add_atoms;
588 >          int old_atoms = atomsPerProc[which_proc];
589 >          int add_atoms = moleculeStamp->getNAtoms();
590 >          int new_atoms = old_atoms + add_atoms;
591            
592            // If we've been through this loop too many times, we need
593            // to just give up and assign the molecule to this processor
594            // and be done with it.
595            
596            if (loops > 100) {
597 +
598              sprintf(painCave.errMsg,
599 <                    "I've tried 100 times to assign molecule %d to a "
600 <                    " processor, but can't find a good spot.\n"
601 <                    "I'm assigning it at random to processor %d.\n",
599 >                    "There have been 100 attempts to assign molecule %d to an\n"
600 >                    "\tunderworked processor, but there's no good place to\n"
601 >                    "\tleave it.  OpenMD is assigning it at random to processor %d.\n",
602                      i, which_proc);
603 <            
603 >          
604              painCave.isFatal = 0;
605 +            painCave.severity = OPENMD_INFO;
606              simError();
607              
608              molToProcMap[i] = which_proc;
# Line 333 | Line 631 | namespace oopse {
631            //           Pacc(x) = exp(- a * x)
632            // where a = penalty / (average atoms per molecule)
633            
634 <          x = (double)(new_atoms - nTarget);
635 <          y = myRandom->rand();
634 >          RealType x = (RealType)(new_atoms - nTarget);
635 >          RealType y = myRandom->rand();
636            
637            if (y < exp(- a * x)) {
638              molToProcMap[i] = which_proc;
# Line 349 | Line 647 | namespace oopse {
647        }
648        
649        delete myRandom;
650 <      
650 >
651        // Spray out this nonsense to all other processors:
354      
652        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
653 +      // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
654      } else {
655        
656        // Listen to your marching orders from processor 0:
359      
657        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
658 +      // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0);
659 +
660      }
661      
662      info->setMolToProcMap(molToProcMap);
663      sprintf(checkPointMsg,
664              "Successfully divided the molecules among the processors.\n");
665 <    MPIcheckPoint();
665 >    errorCheckPoint();
666    }
667    
668   #endif
# Line 380 | Line 679 | namespace oopse {
679   #endif
680          
681          stampId = info->getMoleculeStampId(i);
682 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
683 <                                                   stampId, i, info->getLocalIndexManager());
682 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
683 >                                                   info->getMoleculeStamp(stampId),
684 >                                                   stampId, i,
685 >                                                   info->getLocalIndexManager());
686          
687          info->addMolecule(mol);
688          
# Line 393 | Line 694 | namespace oopse {
694        
695      } //end for(int i=0)  
696    }
396  
397  void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
398                            std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
399    int i;
400    char * id;
401    MoleculeStamp * currentStamp;
402    Component** the_components = simParams->getComponents();
403    int n_components = simParams->getNComponents();
697      
698 <    if (!simParams->haveNMol()) {
406 <      // we don't have the total number of molecules, so we assume it is
407 <      // given in each component
408 <      
409 <      for(i = 0; i < n_components; i++) {
410 <        if (!the_components[i]->haveNMol()) {
411 <          // we have a problem
412 <          sprintf(painCave.errMsg,
413 <                  "SimCreator Error. No global NMol or component NMol given.\n"
414 <                  "\tCannot calculate the number of atoms.\n");
415 <          
416 <          painCave.isFatal = 1;
417 <          simError();
418 <        }
419 <        
420 <        id = the_components[i]->getType();
698 >  int SimCreator::computeStorageLayout(SimInfo* info) {
699  
700 <        currentStamp = stamps->getMolStamp(id);
701 <        if (currentStamp == NULL) {
702 <          sprintf(painCave.errMsg,
703 <                  "SimCreator error: Component \"%s\" was not found in the "
704 <                  "list of declared molecules\n", id);
705 <          
706 <          painCave.isFatal = 1;
707 <          simError();
708 <        }
709 <        
710 <        moleculeStampPairs.push_back(
711 <                                     std::make_pair(currentStamp, the_components[i]->getNMol()));
712 <      } //end for (i = 0; i < n_components; i++)
713 <    } else {
714 <      sprintf(painCave.errMsg, "SimSetup error.\n"
715 <              "\tSorry, the ability to specify total"
716 <              " nMols and then give molfractions in the components\n"
717 <              "\tis not currently supported."
718 <              " Please give nMol in the components.\n");
719 <      
720 <      painCave.isFatal = 1;
721 <      simError();
700 >    Globals* simParams = info->getSimParams();
701 >    int nRigidBodies = info->getNGlobalRigidBodies();
702 >    set<AtomType*> atomTypes = info->getSimulatedAtomTypes();
703 >    set<AtomType*>::iterator i;
704 >    bool hasDirectionalAtoms = false;
705 >    bool hasFixedCharge = false;
706 >    bool hasDipoles = false;    
707 >    bool hasQuadrupoles = false;    
708 >    bool hasPolarizable = false;    
709 >    bool hasFluctuatingCharge = false;    
710 >    bool hasMetallic = false;
711 >    int storageLayout = 0;
712 >    storageLayout |= DataStorage::dslPosition;
713 >    storageLayout |= DataStorage::dslVelocity;
714 >    storageLayout |= DataStorage::dslForce;
715 >
716 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
717 >
718 >      DirectionalAdapter da = DirectionalAdapter( (*i) );
719 >      MultipoleAdapter ma = MultipoleAdapter( (*i) );
720 >      EAMAdapter ea = EAMAdapter( (*i) );
721 >      SuttonChenAdapter sca = SuttonChenAdapter( (*i) );
722 >      PolarizableAdapter pa = PolarizableAdapter( (*i) );
723 >      FixedChargeAdapter fca = FixedChargeAdapter( (*i) );
724 >      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) );
725 >
726 >      if (da.isDirectional()){
727 >        hasDirectionalAtoms = true;
728 >      }
729 >      if (ma.isDipole()){
730 >        hasDipoles = true;
731 >      }
732 >      if (ma.isQuadrupole()){
733 >        hasQuadrupoles = true;
734 >      }
735 >      if (ea.isEAM() || sca.isSuttonChen()){
736 >        hasMetallic = true;
737 >      }
738 >      if ( fca.isFixedCharge() ){
739 >        hasFixedCharge = true;
740 >      }
741 >      if ( fqa.isFluctuatingCharge() ){
742 >        hasFluctuatingCharge = true;
743 >      }
744 >      if ( pa.isPolarizable() ){
745 >        hasPolarizable = true;
746 >      }
747      }
748      
749 < #ifdef IS_MPI
749 >    if (nRigidBodies > 0 || hasDirectionalAtoms) {
750 >      storageLayout |= DataStorage::dslAmat;
751 >      if(storageLayout & DataStorage::dslVelocity) {
752 >        storageLayout |= DataStorage::dslAngularMomentum;
753 >      }
754 >      if (storageLayout & DataStorage::dslForce) {
755 >        storageLayout |= DataStorage::dslTorque;
756 >      }
757 >    }
758 >    if (hasDipoles) {
759 >      storageLayout |= DataStorage::dslDipole;
760 >    }
761 >    if (hasQuadrupoles) {
762 >      storageLayout |= DataStorage::dslQuadrupole;
763 >    }
764 >    if (hasFixedCharge || hasFluctuatingCharge) {
765 >      storageLayout |= DataStorage::dslSkippedCharge;
766 >    }
767 >    if (hasMetallic) {
768 >      storageLayout |= DataStorage::dslDensity;
769 >      storageLayout |= DataStorage::dslFunctional;
770 >      storageLayout |= DataStorage::dslFunctionalDerivative;
771 >    }
772 >    if (hasPolarizable) {
773 >      storageLayout |= DataStorage::dslElectricField;
774 >    }
775 >    if (hasFluctuatingCharge){
776 >      storageLayout |= DataStorage::dslFlucQPosition;
777 >      if(storageLayout & DataStorage::dslVelocity) {
778 >        storageLayout |= DataStorage::dslFlucQVelocity;
779 >      }
780 >      if (storageLayout & DataStorage::dslForce) {
781 >        storageLayout |= DataStorage::dslFlucQForce;
782 >      }
783 >    }
784      
785 <    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
786 <    MPIcheckPoint();
787 <    
788 < #endif // is_mpi
789 <    
785 >    // if the user has asked for them, make sure we've got the memory for the
786 >    // objects defined.
787 >
788 >    if (simParams->getOutputParticlePotential()) {
789 >      storageLayout |= DataStorage::dslParticlePot;
790 >    }
791 >
792 >    if (simParams->havePrintHeatFlux()) {
793 >      if (simParams->getPrintHeatFlux()) {
794 >        storageLayout |= DataStorage::dslParticlePot;
795 >      }
796 >    }
797 >
798 >    if (simParams->getOutputElectricField() | simParams->haveElectricField()) {
799 >      storageLayout |= DataStorage::dslElectricField;
800 >    }
801 >
802 >    if (simParams->getOutputFluctuatingCharges()) {
803 >      storageLayout |= DataStorage::dslFlucQPosition;
804 >      storageLayout |= DataStorage::dslFlucQVelocity;
805 >      storageLayout |= DataStorage::dslFlucQForce;
806 >    }
807 >
808 >    info->setStorageLayout(storageLayout);
809 >
810 >    return storageLayout;
811    }
812 <  
812 >
813    void SimCreator::setGlobalIndex(SimInfo *info) {
814      SimInfo::MoleculeIterator mi;
815      Molecule::AtomIterator ai;
816      Molecule::RigidBodyIterator ri;
817      Molecule::CutoffGroupIterator ci;
818 <    Molecule * mol;
819 <    Atom * atom;
820 <    RigidBody * rb;
821 <    CutoffGroup * cg;
818 >    Molecule::BondIterator boi;
819 >    Molecule::BendIterator bei;
820 >    Molecule::TorsionIterator ti;
821 >    Molecule::InversionIterator ii;
822 >    Molecule::IntegrableObjectIterator  ioi;
823 >    Molecule* mol;
824 >    Atom* atom;
825 >    RigidBody* rb;
826 >    CutoffGroup* cg;
827 >    Bond* bond;
828 >    Bend* bend;
829 >    Torsion* torsion;
830 >    Inversion* inversion;
831      int beginAtomIndex;
832      int beginRigidBodyIndex;
833      int beginCutoffGroupIndex;
834 +    int beginBondIndex;
835 +    int beginBendIndex;
836 +    int beginTorsionIndex;
837 +    int beginInversionIndex;
838      int nGlobalAtoms = info->getNGlobalAtoms();
839 +    int nGlobalRigidBodies = info->getNGlobalRigidBodies();
840      
469 #ifndef IS_MPI
470    
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 <    
480 <    myNode = worldRank;
481 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
482 <    
483 <    std::vector < int > tmpAtomsInProc(nproc, 0);
484 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
485 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
486 <    std::vector < int > NumAtomsInProc(nproc, 0);
487 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
488 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
489 <    
490 <    tmpAtomsInProc[myNode] = info->getNAtoms();
491 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
492 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
493 <    
494 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
495 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
496 <                  MPI_SUM, MPI_COMM_WORLD);
497 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
498 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
499 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
500 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
501 <    
502 <    beginAtomIndex = 0;
503 <    beginRigidBodyIndex = 0;
504 <    beginCutoffGroupIndex = 0;
505 <    
506 <    for(int i = 0; i < myNode; i++) {
507 <      beginAtomIndex += NumAtomsInProc[i];
508 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
509 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
510 <    }
511 <    
512 < #endif
513 <    
514 <    //rigidbody's index begins right after atom's
515 <    beginRigidBodyIndex += info->getNGlobalAtoms();
516 <    
517 <    for(mol = info->beginMolecule(mi); mol != NULL;
518 <        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++);
528 <      }
529 <      
530 <      //local index of cutoff group is trivial, it only depends on the order of travesing
531 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
532 <          cg = mol->nextCutoffGroup(ci)) {
533 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
534 <      }
535 <    }
536 <    
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 +    // MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0],
936 +    //                           &tmpGroupMembership[0], nGlobalAtoms,
937 +    //                           MPI::INT, MPI::SUM);
938      info->setGlobalGroupMembership(tmpGroupMembership);
939   #else
940      info->setGlobalGroupMembership(globalGroupMembership);
941   #endif
942      
943      //fill molMembership
944 <    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
944 >    std::vector<int> globalMolMembership(info->getNGlobalAtoms() +
945 >                                         info->getNGlobalRigidBodies(), 0);
946      
947 <    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
948 <      
947 >    for(mol = info->beginMolecule(mi); mol != NULL;
948 >        mol = info->nextMolecule(mi)) {
949        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
950          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
951        }
952 +      for (rb = mol->beginRigidBody(ri); rb != NULL;
953 +           rb = mol->nextRigidBody(ri)) {
954 +        globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex();
955 +      }
956      }
957      
958   #ifdef IS_MPI
959 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
960 <    
961 <    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
959 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms() +
960 >                                      info->getNGlobalRigidBodies(), 0);
961 >    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
962 >                  nGlobalAtoms + nGlobalRigidBodies,
963                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
964 +    // MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0],
965 +    //                           nGlobalAtoms + nGlobalRigidBodies,
966 +    //                           MPI::INT, MPI::SUM);
967      
968      info->setGlobalMolMembership(tmpMolMembership);
969   #else
970      info->setGlobalMolMembership(globalMolMembership);
971   #endif
972 +
973 +    // nIOPerMol holds the number of integrable objects per molecule
974 +    // here the molecules are listed by their global indices.
975 +
976 +    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
977 +    for (mol = info->beginMolecule(mi); mol != NULL;
978 +         mol = info->nextMolecule(mi)) {
979 +      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
980 +    }
981      
982 <  }
983 <  
984 <  void SimCreator::loadCoordinates(SimInfo* info) {
985 <    Globals* simParams;
986 <    simParams = info->getSimParams();
982 > #ifdef IS_MPI
983 >    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
984 >    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
985 >      info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
986 >    // MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
987 >    //                           info->getNGlobalMolecules(), MPI::INT, MPI::SUM);
988 > #else
989 >    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
990 > #endif    
991 >
992 >    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
993      
994 <    if (!simParams->haveInitialConfig()) {
995 <      sprintf(painCave.errMsg,
996 <              "Cannot intialize a simulation without an initial configuration file.\n");
997 <      painCave.isFatal = 1;;
594 <      simError();
994 >    int startingIndex = 0;
995 >    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
996 >      startingIOIndexForMol[i] = startingIndex;
997 >      startingIndex += numIntegrableObjectsPerMol[i];
998      }
999      
1000 <    DumpReader reader(info, simParams->getInitialConfig());
1000 >    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
1001 >    for (mol = info->beginMolecule(mi); mol != NULL;
1002 >         mol = info->nextMolecule(mi)) {
1003 >      int myGlobalIndex = mol->getGlobalIndex();
1004 >      int globalIO = startingIOIndexForMol[myGlobalIndex];
1005 >      for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL;
1006 >           sd = mol->nextIntegrableObject(ioi)) {
1007 >        sd->setGlobalIntegrableObjectIndex(globalIO);
1008 >        IOIndexToIntegrableObject[globalIO] = sd;
1009 >        globalIO++;
1010 >      }
1011 >    }
1012 >      
1013 >    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
1014 >    
1015 >  }
1016 >  
1017 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
1018 >    
1019 >    DumpReader reader(info, mdFileName);
1020      int nframes = reader.getNFrames();
1021      
1022      if (nframes > 0) {
# Line 603 | Line 1025 | namespace oopse {
1025        //invalid initial coordinate file
1026        sprintf(painCave.errMsg,
1027                "Initial configuration file %s should at least contain one frame\n",
1028 <              simParams->getInitialConfig().c_str());
1028 >              mdFileName.c_str());
1029        painCave.isFatal = 1;
1030        simError();
1031      }
610    
1032      //copy the current snapshot to previous snapshot
1033      info->getSnapshotManager()->advance();
1034    }
1035    
1036 < } //end namespace oopse
1036 > } //end namespace OpenMD
1037  
1038  

Comparing trunk/src/brains/SimCreator.cpp (property svn:keywords):
Revision 749 by tim, Wed Nov 16 23:10:02 2005 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines