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

Comparing:
trunk/src/brains/SimCreator.cpp (file contents), Revision 816 by tim, Fri Dec 16 18:26:41 2005 UTC vs.
branches/development/src/brains/SimCreator.cpp (file contents), Revision 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   /**
# Line 46 | Line 46
46   * @time 13:51am
47   * @version 1.0
48   */
49 <
49 > #include <exception>
50   #include <iostream>
51   #include <sstream>
52   #include <string>
# Line 79 | Line 79
79   #include "math/ParallelRandNumGen.hpp"
80   #endif
81  
82 < namespace oopse {
82 > namespace OpenMD {
83    
84 < Globals* SimCreator::parseFile(const std::string mdFileName){
85 <        Globals* simParams = NULL;
86 <        try {
84 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
85 >    Globals* simParams = NULL;
86 >    try {
87  
88 <            // Create a preprocessor that preprocesses md file into an ostringstream
89 <            std::stringstream ppStream;
88 >      // Create a preprocessor that preprocesses md file into an ostringstream
89 >      std::stringstream ppStream;
90   #ifdef IS_MPI            
91 <            int streamSize;
92 <            const int masterNode = 0;
93 <            int commStatus;
94 <            if (worldRank == masterNode) {
95 < #endif
91 >      int streamSize;
92 >      const int masterNode = 0;
93 >      int commStatus;
94 >      if (worldRank == masterNode) {
95 >        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
96 > #endif                
97 >        SimplePreprocessor preprocessor;
98 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
99                  
97                SimplePreprocessor preprocessor;
98                preprocessor.preprocess(mdFileName, ppStream);
99                
100   #ifdef IS_MPI            
101 <                //brocasting the stream size
102 <                streamSize = ppStream.str().size() +1;
103 <                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
101 >        //brocasting the stream size
102 >        streamSize = ppStream.str().size() +1;
103 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
104  
105 <                commStatus = MPI_Bcast(ppStream.str().c_str(), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
105 >        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
106              
107                  
108 <            } else {
109 <                //get stream size
110 <                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
108 >      } else {
109 >
110 >        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
111 >
112 >        //get stream size
113 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
114 >
115 >        char* buf = new char[streamSize];
116 >        assert(buf);
117                  
118 <                  char* buf = new char[streamSize];
119 <                  assert(buf);
118 >        //receive file content
119 >        commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
120                  
121 <                  //receive file content
122 <                  commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
117 <                
118 <                  ppStream.str(buf);
119 <                  delete buf;
121 >        ppStream.str(buf);
122 >        delete [] buf;
123  
124 <            }
124 >      }
125   #endif            
126 <            // Create a scanner that reads from the input stream
127 <            MDLexer lexer(ppStream);
128 <            lexer.setFilename(mdFileName);
129 <            lexer.initDeferredLineCount();
126 >      // Create a scanner that reads from the input stream
127 >      MDLexer lexer(ppStream);
128 >      lexer.setFilename(filename);
129 >      lexer.initDeferredLineCount();
130      
131 <            // Create a parser that reads from the scanner
132 <            MDParser parser(lexer);
133 <            parser.setFilename(mdFileName);
131 >      // Create a parser that reads from the scanner
132 >      MDParser parser(lexer);
133 >      parser.setFilename(filename);
134  
135 <            // Create an observer that synchorizes file name change
136 <            FilenameObserver observer;
137 <            observer.setLexer(&lexer);
138 <            observer.setParser(&parser);
139 <            lexer.setObserver(&observer);
135 >      // Create an observer that synchorizes file name change
136 >      FilenameObserver observer;
137 >      observer.setLexer(&lexer);
138 >      observer.setParser(&parser);
139 >      lexer.setObserver(&observer);
140      
141 <            antlr::ASTFactory factory;
142 <            parser.initializeASTFactory(factory);
143 <            parser.setASTFactory(&factory);
144 <            parser.mdfile();
141 >      antlr::ASTFactory factory;
142 >      parser.initializeASTFactory(factory);
143 >      parser.setASTFactory(&factory);
144 >      parser.mdfile();
145  
146 <            // Create a tree parser that reads information into Globals
147 <            MDTreeParser treeParser;
148 <            treeParser.initializeASTFactory(factory);
149 <            treeParser.setASTFactory(&factory);
150 <             simParams = treeParser.walkTree(parser.getAST());
146 >      // Create a tree parser that reads information into Globals
147 >      MDTreeParser treeParser;
148 >      treeParser.initializeASTFactory(factory);
149 >      treeParser.setASTFactory(&factory);
150 >      simParams = treeParser.walkTree(parser.getAST());
151 >    }
152  
149        }
153        
154 <      catch(antlr::MismatchedCharException& e) {
155 <          cerr<< "parser exception: " << e.getMessage() << " " <<  e.getFilename() << ":" << e.getLine() << " " << e.getColumn() << endl;
156 <      }
157 <      catch(antlr::MismatchedTokenException &e) {
158 <          cerr<< "parser exception: " << e.getMessage() << " " <<  e.getFilename() << ":" << e.getLine() << " " << e.getColumn() << endl;
159 <      }
160 <      catch(antlr::NoViableAltForCharException &e) {
161 <          cerr<< "parser exception: " << e.getMessage() << " " <<  e.getFilename() << ":" << e.getLine() << " " << e.getColumn() << endl;
162 <      }
163 <      catch(antlr::NoViableAltException &e) {
164 <          cerr<< "parser exception: " << e.getMessage() << " " <<  e.getFilename() << ":" << e.getLine() << " " << e.getColumn() << endl;
165 <      }
166 <        catch(antlr::TokenStreamRecognitionException& e) {
167 <          cerr<< "parser exception: " << e.getMessage() << " " <<  e.getFilename() << ":" << e.getLine() << " " << e.getColumn() << endl;
168 <        }
169 <        catch(antlr::TokenStreamIOException& e) {
170 <          cerr<< "parser exception: " << e.getMessage() << endl;
171 <        }
172 <        catch(antlr::TokenStreamException& e) {
173 <          cerr<< "parser exception: " << e.getMessage() << endl;
174 <        }        
175 <       catch (antlr::RecognitionException& e) {
176 <          cerr<< "parser exception: " << e.getMessage() << " " <<  e.getFilename() << ":" << e.getLine() << " " << e.getColumn() << endl;
177 <       }
178 <       catch (antlr::CharStreamException& e) {
179 <            cerr << "parser exception: " << e.getMessage() << endl;
180 <       }
181 <        catch (exception& e) {
182 <            cerr << "parser exception: " << e.what() << endl;
183 <        }
154 >    catch(antlr::MismatchedCharException& e) {
155 >      sprintf(painCave.errMsg,
156 >              "parser exception: %s %s:%d:%d\n",
157 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
158 >      painCave.isFatal = 1;
159 >      simError();          
160 >    }
161 >    catch(antlr::MismatchedTokenException &e) {
162 >      sprintf(painCave.errMsg,
163 >              "parser exception: %s %s:%d:%d\n",
164 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
165 >      painCave.isFatal = 1;
166 >      simError();  
167 >    }
168 >    catch(antlr::NoViableAltForCharException &e) {
169 >      sprintf(painCave.errMsg,
170 >              "parser exception: %s %s:%d:%d\n",
171 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
172 >      painCave.isFatal = 1;
173 >      simError();  
174 >    }
175 >    catch(antlr::NoViableAltException &e) {
176 >      sprintf(painCave.errMsg,
177 >              "parser exception: %s %s:%d:%d\n",
178 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
179 >      painCave.isFatal = 1;
180 >      simError();  
181 >    }
182 >      
183 >    catch(antlr::TokenStreamRecognitionException& 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::TokenStreamIOException& e) {
192 >      sprintf(painCave.errMsg,
193 >              "parser exception: %s\n",
194 >              e.getMessage().c_str());
195 >      painCave.isFatal = 1;
196 >      simError();
197 >    }
198 >        
199 >    catch(antlr::TokenStreamException& e) {
200 >      sprintf(painCave.errMsg,
201 >              "parser exception: %s\n",
202 >              e.getMessage().c_str());
203 >      painCave.isFatal = 1;
204 >      simError();
205 >    }        
206 >    catch (antlr::RecognitionException& e) {
207 >      sprintf(painCave.errMsg,
208 >              "parser exception: %s %s:%d:%d\n",
209 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
210 >      painCave.isFatal = 1;
211 >      simError();          
212 >    }
213 >    catch (antlr::CharStreamException& e) {
214 >      sprintf(painCave.errMsg,
215 >              "parser exception: %s\n",
216 >              e.getMessage().c_str());
217 >      painCave.isFatal = 1;
218 >      simError();        
219 >    }
220 >    catch (OpenMDException& e) {
221 >      sprintf(painCave.errMsg,
222 >              "%s\n",
223 >              e.getMessage().c_str());
224 >      painCave.isFatal = 1;
225 >      simError();
226 >    }
227 >    catch (std::exception& e) {
228 >      sprintf(painCave.errMsg,
229 >              "parser exception: %s\n",
230 >              e.what());
231 >      painCave.isFatal = 1;
232 >      simError();
233 >    }
234  
235 <        return simParams;
235 >    simParams->setMDfileVersion(mdFileVersion);
236 >    return simParams;
237    }
238    
239    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
240                                    bool loadInitCoords) {
241 +    
242 +    const int bufferSize = 65535;
243 +    char buffer[bufferSize];
244 +    int lineNo = 0;
245 +    std::string mdRawData;
246 +    int metaDataBlockStart = -1;
247 +    int metaDataBlockEnd = -1;
248 +    int i;
249 +    int mdOffset;
250 +    int mdFileVersion;
251 +
252 + #ifdef IS_MPI            
253 +    const int masterNode = 0;
254 +    if (worldRank == masterNode) {
255 + #endif
256 +
257 +      std::ifstream mdFile_(mdFileName.c_str());
258 +      
259 +      if (mdFile_.fail()) {
260 +        sprintf(painCave.errMsg,
261 +                "SimCreator: Cannot open file: %s\n",
262 +                mdFileName.c_str());
263 +        painCave.isFatal = 1;
264 +        simError();
265 +      }
266 +
267 +      mdFile_.getline(buffer, bufferSize);
268 +      ++lineNo;
269 +      std::string line = trimLeftCopy(buffer);
270 +      i = CaseInsensitiveFind(line, "<OpenMD");
271 +      if (static_cast<size_t>(i) == string::npos) {
272 +        // try the older file strings to see if that works:
273 +        i = CaseInsensitiveFind(line, "<OOPSE");
274 +      }
275 +      
276 +      if (static_cast<size_t>(i) == string::npos) {
277 +        // still no luck!
278 +        sprintf(painCave.errMsg,
279 +                "SimCreator: File: %s is not a valid OpenMD file!\n",
280 +                mdFileName.c_str());
281 +        painCave.isFatal = 1;
282 +        simError();
283 +      }
284 +      
285 +      // found the correct opening string, now try to get the file
286 +      // format version number.
287 +
288 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
289 +      std::string fileType = tokenizer.nextToken();
290 +      toUpper(fileType);
291 +
292 +      mdFileVersion = 0;
293 +
294 +      if (fileType == "OPENMD") {
295 +        while (tokenizer.hasMoreTokens()) {
296 +          std::string token(tokenizer.nextToken());
297 +          toUpper(token);
298 +          if (token == "VERSION") {
299 +            mdFileVersion = tokenizer.nextTokenAsInt();
300 +            break;
301 +          }
302 +        }
303 +      }
304 +            
305 +      //scan through the input stream and find MetaData tag        
306 +      while(mdFile_.getline(buffer, bufferSize)) {
307 +        ++lineNo;
308 +        
309 +        std::string line = trimLeftCopy(buffer);
310 +        if (metaDataBlockStart == -1) {
311 +          i = CaseInsensitiveFind(line, "<MetaData>");
312 +          if (i != string::npos) {
313 +            metaDataBlockStart = lineNo;
314 +            mdOffset = mdFile_.tellg();
315 +          }
316 +        } else {
317 +          i = CaseInsensitiveFind(line, "</MetaData>");
318 +          if (i != string::npos) {
319 +            metaDataBlockEnd = lineNo;
320 +          }
321 +        }
322 +      }
323  
324 +      if (metaDataBlockStart == -1) {
325 +        sprintf(painCave.errMsg,
326 +                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
327 +                mdFileName.c_str());
328 +        painCave.isFatal = 1;
329 +        simError();
330 +      }
331 +      if (metaDataBlockEnd == -1) {
332 +        sprintf(painCave.errMsg,
333 +                "SimCreator: File: %s did not contain a closed MetaData block!\n",
334 +                mdFileName.c_str());
335 +        painCave.isFatal = 1;
336 +        simError();
337 +      }
338 +        
339 +      mdFile_.clear();
340 +      mdFile_.seekg(0);
341 +      mdFile_.seekg(mdOffset);
342 +
343 +      mdRawData.clear();
344 +
345 +      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
346 +        mdFile_.getline(buffer, bufferSize);
347 +        mdRawData += buffer;
348 +        mdRawData += "\n";
349 +      }
350 +
351 +      mdFile_.close();
352 +
353 + #ifdef IS_MPI
354 +    }
355 + #endif
356 +
357 +    std::stringstream rawMetaDataStream(mdRawData);
358 +
359      //parse meta-data file
360 <    Globals* simParams = parseFile(mdFileName);
360 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
361 >                                   metaDataBlockStart + 1);
362      
363      //create the force field
364 <    ForceField * ff = ForceFieldFactory::getInstance()
365 <      ->createForceField(simParams->getForceField());
194 <    
364 >    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
365 >
366      if (ff == NULL) {
367        sprintf(painCave.errMsg,
368                "ForceField Factory can not create %s force field\n",
# Line 224 | Line 395 | Globals* SimCreator::parseFile(const std::string mdFil
395      }
396      
397      ff->parse(forcefieldFileName);
227        
398      //create SimInfo
399      SimInfo * info = new SimInfo(ff, simParams);
400 +
401 +    info->setRawMetaData(mdRawData);
402      
403 <    //gather parameters (SimCreator only retrieves part of the parameters)
403 >    //gather parameters (SimCreator only retrieves part of the
404 >    //parameters)
405      gatherParameters(info, mdFileName);
406      
407      //divide the molecules and determine the global index of molecules
# Line 239 | Line 412 | Globals* SimCreator::parseFile(const std::string mdFil
412      //create the molecules
413      createMolecules(info);
414      
415 <    
416 <    //allocate memory for DataStorage(circular reference, need to break it)
415 >    //allocate memory for DataStorage(circular reference, need to
416 >    //break it)
417      info->setSnapshotManager(new SimSnapshotManager(info));
418      
419 <    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
420 <    //global index will never change again). Local indices of atoms and rigidbodies are already set by
421 <    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
419 >    //set the global index of atoms, rigidbodies and cutoffgroups
420 >    //(only need to be set once, the global index will never change
421 >    //again). Local indices of atoms and rigidbodies are already set
422 >    //by MoleculeCreator class which actually delegates the
423 >    //responsibility to LocalIndexManager.
424      setGlobalIndex(info);
425      
426 <    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
427 <    //atoms don't have the global index yet  (their global index are all initialized to -1).
428 <    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
429 <    //we can determine the beginning global indices of atoms before they get created.
426 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
427 >    //method, at that point atoms don't have the global index yet
428 >    //(their global index are all initialized to -1).  Therefore we
429 >    //have to call addInteractionPairs explicitly here. A way to work
430 >    //around is that we can determine the beginning global indices of
431 >    //atoms before they get created.
432      SimInfo::MoleculeIterator mi;
433      Molecule* mol;
434      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
435 <      info->addExcludePairs(mol);
435 >      info->addInteractionPairs(mol);
436      }
437      
438      if (loadInitCoords)
439 <      loadCoordinates(info);    
263 <    
439 >      loadCoordinates(info, mdFileName);    
440      return info;
441    }
442    
# Line 295 | Line 471 | Globals* SimCreator::parseFile(const std::string mdFil
471    
472   #ifdef IS_MPI
473    void SimCreator::divideMolecules(SimInfo *info) {
474 <    double numerator;
475 <    double denominator;
476 <    double precast;
477 <    double x;
478 <    double y;
479 <    double a;
474 >    RealType numerator;
475 >    RealType denominator;
476 >    RealType precast;
477 >    RealType x;
478 >    RealType y;
479 >    RealType a;
480      int old_atoms;
481      int add_atoms;
482      int new_atoms;
# Line 324 | Line 500 | Globals* SimCreator::parseFile(const std::string mdFil
500                "\tthe number of molecules.  This will not result in a \n"
501                "\tusable division of atoms for force decomposition.\n"
502                "\tEither try a smaller number of processors, or run the\n"
503 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
503 >              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
504        
505        painCave.isFatal = 1;
506        simError();
# Line 412 | Line 588 | Globals* SimCreator::parseFile(const std::string mdFil
588            //           Pacc(x) = exp(- a * x)
589            // where a = penalty / (average atoms per molecule)
590            
591 <          x = (double)(new_atoms - nTarget);
591 >          x = (RealType)(new_atoms - nTarget);
592            y = myRandom->rand();
593            
594            if (y < exp(- a * x)) {
# Line 442 | Line 618 | Globals* SimCreator::parseFile(const std::string mdFil
618      info->setMolToProcMap(molToProcMap);
619      sprintf(checkPointMsg,
620              "Successfully divided the molecules among the processors.\n");
621 <    MPIcheckPoint();
621 >    errorCheckPoint();
622    }
623    
624   #endif
# Line 459 | Line 635 | Globals* SimCreator::parseFile(const std::string mdFil
635   #endif
636          
637          stampId = info->getMoleculeStampId(i);
638 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
639 <                                                   stampId, i, info->getLocalIndexManager());
638 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
639 >                                                   info->getMoleculeStamp(stampId),
640 >                                                   stampId, i,
641 >                                                   info->getLocalIndexManager());
642          
643          info->addMolecule(mol);
644          
# Line 478 | Line 656 | Globals* SimCreator::parseFile(const std::string mdFil
656      Molecule::AtomIterator ai;
657      Molecule::RigidBodyIterator ri;
658      Molecule::CutoffGroupIterator ci;
659 +    Molecule::IntegrableObjectIterator  ioi;
660      Molecule * mol;
661      Atom * atom;
662      RigidBody * rb;
# Line 487 | Line 666 | Globals* SimCreator::parseFile(const std::string mdFil
666      int beginCutoffGroupIndex;
667      int nGlobalAtoms = info->getNGlobalAtoms();
668      
490 #ifndef IS_MPI
491    
669      beginAtomIndex = 0;
670      beginRigidBodyIndex = 0;
671      beginCutoffGroupIndex = 0;
672 <    
673 < #else
497 <    
498 <    int nproc;
499 <    int myNode;
500 <    
501 <    myNode = worldRank;
502 <    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
503 <    
504 <    std::vector < int > tmpAtomsInProc(nproc, 0);
505 <    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
506 <    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
507 <    std::vector < int > NumAtomsInProc(nproc, 0);
508 <    std::vector < int > NumRigidBodiesInProc(nproc, 0);
509 <    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
510 <    
511 <    tmpAtomsInProc[myNode] = info->getNAtoms();
512 <    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
513 <    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
514 <    
515 <    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
516 <    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
517 <                  MPI_SUM, MPI_COMM_WORLD);
518 <    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
519 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
520 <    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
521 <                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
522 <    
523 <    beginAtomIndex = 0;
524 <    beginRigidBodyIndex = 0;
525 <    beginCutoffGroupIndex = 0;
526 <    
527 <    for(int i = 0; i < myNode; i++) {
528 <      beginAtomIndex += NumAtomsInProc[i];
529 <      beginRigidBodyIndex += NumRigidBodiesInProc[i];
530 <      beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
531 <    }
532 <    
533 < #endif
534 <    
535 <    //rigidbody's index begins right after atom's
536 <    beginRigidBodyIndex += info->getNGlobalAtoms();
537 <    
538 <    for(mol = info->beginMolecule(mi); mol != NULL;
539 <        mol = info->nextMolecule(mi)) {
672 >
673 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
674        
675 <      //local index(index in DataStorge) of atom is important
676 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
677 <        atom->setGlobalIndex(beginAtomIndex++);
675 > #ifdef IS_MPI      
676 >      if (info->getMolToProc(i) == worldRank) {
677 > #endif        
678 >        // stuff to do if I own this molecule
679 >        mol = info->getMoleculeByGlobalIndex(i);
680 >
681 >        //local index(index in DataStorge) of atom is important
682 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
683 >          atom->setGlobalIndex(beginAtomIndex++);
684 >        }
685 >        
686 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
687 >            rb = mol->nextRigidBody(ri)) {
688 >          rb->setGlobalIndex(beginRigidBodyIndex++);
689 >        }
690 >        
691 >        //local index of cutoff group is trivial, it only depends on
692 >        //the order of travesing
693 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
694 >            cg = mol->nextCutoffGroup(ci)) {
695 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
696 >        }        
697 >        
698 > #ifdef IS_MPI        
699 >      }  else {
700 >
701 >        // stuff to do if I don't own this molecule
702 >        
703 >        int stampId = info->getMoleculeStampId(i);
704 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
705 >
706 >        beginAtomIndex += stamp->getNAtoms();
707 >        beginRigidBodyIndex += stamp->getNRigidBodies();
708 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
709        }
710 <      
711 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
712 <          rb = mol->nextRigidBody(ri)) {
713 <        rb->setGlobalIndex(beginRigidBodyIndex++);
549 <      }
550 <      
551 <      //local index of cutoff group is trivial, it only depends on the order of travesing
552 <      for(cg = mol->beginCutoffGroup(ci); cg != NULL;
553 <          cg = mol->nextCutoffGroup(ci)) {
554 <        cg->setGlobalIndex(beginCutoffGroupIndex++);
555 <      }
556 <    }
557 <    
710 > #endif          
711 >
712 >    } //end for(int i=0)  
713 >
714      //fill globalGroupMembership
715      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
716      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
# Line 566 | Line 722 | Globals* SimCreator::parseFile(const std::string mdFil
722          
723        }      
724      }
725 <    
725 >  
726   #ifdef IS_MPI    
727      // Since the globalGroupMembership has been zero filled and we've only
728      // poked values into the atoms we know, we can do an Allreduce
729      // to get the full globalGroupMembership array (We think).
730      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
731      // docs said we could.
732 <    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
732 >    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
733      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
734                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
735      info->setGlobalGroupMembership(tmpGroupMembership);
# Line 585 | Line 741 | Globals* SimCreator::parseFile(const std::string mdFil
741      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
742      
743      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
588      
744        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
745          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
746        }
747      }
748      
749   #ifdef IS_MPI
750 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
750 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
751      
752      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
753                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
# Line 601 | Line 756 | Globals* SimCreator::parseFile(const std::string mdFil
756   #else
757      info->setGlobalMolMembership(globalMolMembership);
758   #endif
759 +
760 +    // nIOPerMol holds the number of integrable objects per molecule
761 +    // here the molecules are listed by their global indices.
762 +
763 +    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
764 +    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
765 +      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
766 +    }
767      
768 + #ifdef IS_MPI
769 +    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
770 +    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
771 +                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
772 + #else
773 +    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
774 + #endif    
775 +
776 +    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
777 +    
778 +    int startingIndex = 0;
779 +    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
780 +      startingIOIndexForMol[i] = startingIndex;
781 +      startingIndex += numIntegrableObjectsPerMol[i];
782 +    }
783 +    
784 +    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
785 +    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
786 +      int myGlobalIndex = mol->getGlobalIndex();
787 +      int globalIO = startingIOIndexForMol[myGlobalIndex];
788 +      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
789 +           integrableObject = mol->nextIntegrableObject(ioi)) {
790 +        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
791 +        IOIndexToIntegrableObject[globalIO] = integrableObject;
792 +        globalIO++;
793 +      }
794 +    }
795 +      
796 +    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
797 +    
798    }
799    
800 <  void SimCreator::loadCoordinates(SimInfo* info) {
800 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
801      Globals* simParams;
802 +
803      simParams = info->getSimParams();
804      
805 <    if (!simParams->haveInitialConfig()) {
612 <      sprintf(painCave.errMsg,
613 <              "Cannot intialize a simulation without an initial configuration file.\n");
614 <      painCave.isFatal = 1;;
615 <      simError();
616 <    }
617 <    
618 <    DumpReader reader(info, simParams->getInitialConfig());
805 >    DumpReader reader(info, mdFileName);
806      int nframes = reader.getNFrames();
807 <    
807 >
808      if (nframes > 0) {
809        reader.readFrame(nframes - 1);
810      } else {
811        //invalid initial coordinate file
812        sprintf(painCave.errMsg,
813                "Initial configuration file %s should at least contain one frame\n",
814 <              simParams->getInitialConfig().c_str());
814 >              mdFileName.c_str());
815        painCave.isFatal = 1;
816        simError();
817      }
631    
818      //copy the current snapshot to previous snapshot
819      info->getSnapshotManager()->advance();
820    }
821    
822 < } //end namespace oopse
822 > } //end namespace OpenMD
823  
824  

Comparing:
trunk/src/brains/SimCreator.cpp (property svn:keywords), Revision 816 by tim, Fri Dec 16 18:26:41 2005 UTC vs.
branches/development/src/brains/SimCreator.cpp (property svn:keywords), Revision 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines