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 945 by gezelter, Tue Apr 25 02:09:01 2006 UTC vs.
branches/development/src/brains/SimCreator.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   /**
# Line 76 | Line 77
77   #include "antlr/NoViableAltException.hpp"
78  
79   #ifdef IS_MPI
80 + #include "mpi.h"
81   #include "math/ParallelRandNumGen.hpp"
82   #endif
83  
84 < namespace oopse {
84 > namespace OpenMD {
85    
86 < Globals* SimCreator::parseFile(const std::string mdFileName){
87 <        Globals* simParams = NULL;
88 <        try {
86 >  Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){
87 >    Globals* simParams = NULL;
88 >    try {
89  
90 <            // Create a preprocessor that preprocesses md file into an ostringstream
91 <            std::stringstream ppStream;
90 >      // Create a preprocessor that preprocesses md file into an ostringstream
91 >      std::stringstream ppStream;
92   #ifdef IS_MPI            
93 <            int streamSize;
94 <            const int masterNode = 0;
95 <            int commStatus;
96 <            if (worldRank == masterNode) {
97 < #endif
93 >      int streamSize;
94 >      const int masterNode = 0;
95 >      int commStatus;
96 >      if (worldRank == masterNode) {
97 >        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
98 > #endif                
99 >        SimplePreprocessor preprocessor;
100 >        preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream);
101                  
97                SimplePreprocessor preprocessor;
98                preprocessor.preprocess(mdFileName, ppStream);
99                
102   #ifdef IS_MPI            
103 <                //brocasting the stream size
104 <                streamSize = ppStream.str().size() +1;
105 <                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
103 >        //brocasting the stream size
104 >        streamSize = ppStream.str().size() +1;
105 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);                  
106  
107 <                commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
107 >        commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
108              
109                  
110 <            } else {
111 <                //get stream size
112 <                commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
110 >      } else {
111 >
112 >        commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD);
113 >
114 >        //get stream size
115 >        commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD);  
116 >
117 >        char* buf = new char[streamSize];
118 >        assert(buf);
119                  
120 <                  char* buf = new char[streamSize];
121 <                  assert(buf);
120 >        //receive file content
121 >        commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
122                  
123 <                  //receive file content
124 <                  commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD);
117 <                
118 <                  ppStream.str(buf);
119 <                  delete buf;
123 >        ppStream.str(buf);
124 >        delete [] buf;
125  
126 <            }
126 >      }
127   #endif            
128 <            // Create a scanner that reads from the input stream
129 <            MDLexer lexer(ppStream);
130 <            lexer.setFilename(mdFileName);
131 <            lexer.initDeferredLineCount();
128 >      // Create a scanner that reads from the input stream
129 >      MDLexer lexer(ppStream);
130 >      lexer.setFilename(filename);
131 >      lexer.initDeferredLineCount();
132      
133 <            // Create a parser that reads from the scanner
134 <            MDParser parser(lexer);
135 <            parser.setFilename(mdFileName);
133 >      // Create a parser that reads from the scanner
134 >      MDParser parser(lexer);
135 >      parser.setFilename(filename);
136  
137 <            // Create an observer that synchorizes file name change
138 <            FilenameObserver observer;
139 <            observer.setLexer(&lexer);
140 <            observer.setParser(&parser);
141 <            lexer.setObserver(&observer);
137 >      // Create an observer that synchorizes file name change
138 >      FilenameObserver observer;
139 >      observer.setLexer(&lexer);
140 >      observer.setParser(&parser);
141 >      lexer.setObserver(&observer);
142      
143 <            antlr::ASTFactory factory;
144 <            parser.initializeASTFactory(factory);
145 <            parser.setASTFactory(&factory);
146 <            parser.mdfile();
143 >      antlr::ASTFactory factory;
144 >      parser.initializeASTFactory(factory);
145 >      parser.setASTFactory(&factory);
146 >      parser.mdfile();
147  
148 <            // Create a tree parser that reads information into Globals
149 <            MDTreeParser treeParser;
150 <            treeParser.initializeASTFactory(factory);
151 <            treeParser.setASTFactory(&factory);
152 <             simParams = treeParser.walkTree(parser.getAST());
148 >      // Create a tree parser that reads information into Globals
149 >      MDTreeParser treeParser;
150 >      treeParser.initializeASTFactory(factory);
151 >      treeParser.setASTFactory(&factory);
152 >      simParams = treeParser.walkTree(parser.getAST());
153 >    }
154  
149        }
150
155        
156 <      catch(antlr::MismatchedCharException& e) {
157 <          sprintf(painCave.errMsg,
158 <                  "parser exception: %s %s:%d:%d\n",
159 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
160 <          painCave.isFatal = 1;
161 <          simError();          
162 <      }
163 <      catch(antlr::MismatchedTokenException &e) {
164 <          sprintf(painCave.errMsg,
165 <                  "parser exception: %s %s:%d:%d\n",
166 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
167 <          painCave.isFatal = 1;
168 <          simError();  
169 <      }
170 <      catch(antlr::NoViableAltForCharException &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::NoViableAltException &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 <      }
156 >    catch(antlr::MismatchedCharException& e) {
157 >      sprintf(painCave.errMsg,
158 >              "parser exception: %s %s:%d:%d\n",
159 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
160 >      painCave.isFatal = 1;
161 >      simError();          
162 >    }
163 >    catch(antlr::MismatchedTokenException &e) {
164 >      sprintf(painCave.errMsg,
165 >              "parser exception: %s %s:%d:%d\n",
166 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
167 >      painCave.isFatal = 1;
168 >      simError();  
169 >    }
170 >    catch(antlr::NoViableAltForCharException &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::NoViableAltException &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        
185 <        catch(antlr::TokenStreamRecognitionException& e) {
186 <          sprintf(painCave.errMsg,
187 <                  "parser exception: %s %s:%d:%d\n",
188 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
189 <          painCave.isFatal = 1;
190 <          simError();  
191 <        }
185 >    catch(antlr::TokenStreamRecognitionException& e) {
186 >      sprintf(painCave.errMsg,
187 >              "parser exception: %s %s:%d:%d\n",
188 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
189 >      painCave.isFatal = 1;
190 >      simError();  
191 >    }
192          
193 <        catch(antlr::TokenStreamIOException& e) {
194 <          sprintf(painCave.errMsg,
195 <                  "parser exception: %s\n",
196 <                  e.getMessage().c_str());
197 <          painCave.isFatal = 1;
198 <          simError();
199 <        }
193 >    catch(antlr::TokenStreamIOException& e) {
194 >      sprintf(painCave.errMsg,
195 >              "parser exception: %s\n",
196 >              e.getMessage().c_str());
197 >      painCave.isFatal = 1;
198 >      simError();
199 >    }
200          
201 <        catch(antlr::TokenStreamException& e) {
202 <          sprintf(painCave.errMsg,
203 <                  "parser exception: %s\n",
204 <                  e.getMessage().c_str());
205 <          painCave.isFatal = 1;
206 <          simError();
207 <        }        
208 <       catch (antlr::RecognitionException& e) {
209 <          sprintf(painCave.errMsg,
210 <                  "parser exception: %s %s:%d:%d\n",
211 <                  e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
212 <          painCave.isFatal = 1;
213 <          simError();          
214 <       }
215 <       catch (antlr::CharStreamException& e) {
216 <          sprintf(painCave.errMsg,
217 <                  "parser exception: %s\n",
218 <                  e.getMessage().c_str());
219 <          painCave.isFatal = 1;
220 <          simError();        
221 <       }
222 <       catch (OOPSEException& e) {
223 <          sprintf(painCave.errMsg,
224 <                  "%s\n",
225 <                  e.getMessage().c_str());
226 <          painCave.isFatal = 1;
227 <          simError();
228 <       }
229 <       catch (std::exception& e) {
230 <          sprintf(painCave.errMsg,
231 <                  "parser exception: %s\n",
232 <                  e.what());
233 <          painCave.isFatal = 1;
234 <          simError();
235 <       }
201 >    catch(antlr::TokenStreamException& e) {
202 >      sprintf(painCave.errMsg,
203 >              "parser exception: %s\n",
204 >              e.getMessage().c_str());
205 >      painCave.isFatal = 1;
206 >      simError();
207 >    }        
208 >    catch (antlr::RecognitionException& e) {
209 >      sprintf(painCave.errMsg,
210 >              "parser exception: %s %s:%d:%d\n",
211 >              e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn());
212 >      painCave.isFatal = 1;
213 >      simError();          
214 >    }
215 >    catch (antlr::CharStreamException& e) {
216 >      sprintf(painCave.errMsg,
217 >              "parser exception: %s\n",
218 >              e.getMessage().c_str());
219 >      painCave.isFatal = 1;
220 >      simError();        
221 >    }
222 >    catch (OpenMDException& e) {
223 >      sprintf(painCave.errMsg,
224 >              "%s\n",
225 >              e.getMessage().c_str());
226 >      painCave.isFatal = 1;
227 >      simError();
228 >    }
229 >    catch (std::exception& e) {
230 >      sprintf(painCave.errMsg,
231 >              "parser exception: %s\n",
232 >              e.what());
233 >      painCave.isFatal = 1;
234 >      simError();
235 >    }
236  
237 <        return simParams;
237 >    simParams->setMDfileVersion(mdFileVersion);
238 >    return simParams;
239    }
240    
241    SimInfo*  SimCreator::createSim(const std::string & mdFileName,
242                                    bool loadInitCoords) {
243 +    
244 +    const int bufferSize = 65535;
245 +    char buffer[bufferSize];
246 +    int lineNo = 0;
247 +    std::string mdRawData;
248 +    int metaDataBlockStart = -1;
249 +    int metaDataBlockEnd = -1;
250 +    int i;
251 +    int mdOffset;
252 +    int mdFileVersion;
253  
254 + #ifdef IS_MPI            
255 +    const int masterNode = 0;
256 +    if (worldRank == masterNode) {
257 + #endif
258 +
259 +      std::ifstream mdFile_(mdFileName.c_str());
260 +      
261 +      if (mdFile_.fail()) {
262 +        sprintf(painCave.errMsg,
263 +                "SimCreator: Cannot open file: %s\n",
264 +                mdFileName.c_str());
265 +        painCave.isFatal = 1;
266 +        simError();
267 +      }
268 +
269 +      mdFile_.getline(buffer, bufferSize);
270 +      ++lineNo;
271 +      std::string line = trimLeftCopy(buffer);
272 +      i = CaseInsensitiveFind(line, "<OpenMD");
273 +      if (static_cast<size_t>(i) == string::npos) {
274 +        // try the older file strings to see if that works:
275 +        i = CaseInsensitiveFind(line, "<OOPSE");
276 +      }
277 +      
278 +      if (static_cast<size_t>(i) == string::npos) {
279 +        // still no luck!
280 +        sprintf(painCave.errMsg,
281 +                "SimCreator: File: %s is not a valid OpenMD file!\n",
282 +                mdFileName.c_str());
283 +        painCave.isFatal = 1;
284 +        simError();
285 +      }
286 +      
287 +      // found the correct opening string, now try to get the file
288 +      // format version number.
289 +
290 +      StringTokenizer tokenizer(line, "=<> \t\n\r");
291 +      std::string fileType = tokenizer.nextToken();
292 +      toUpper(fileType);
293 +
294 +      mdFileVersion = 0;
295 +
296 +      if (fileType == "OPENMD") {
297 +        while (tokenizer.hasMoreTokens()) {
298 +          std::string token(tokenizer.nextToken());
299 +          toUpper(token);
300 +          if (token == "VERSION") {
301 +            mdFileVersion = tokenizer.nextTokenAsInt();
302 +            break;
303 +          }
304 +        }
305 +      }
306 +            
307 +      //scan through the input stream and find MetaData tag        
308 +      while(mdFile_.getline(buffer, bufferSize)) {
309 +        ++lineNo;
310 +        
311 +        std::string line = trimLeftCopy(buffer);
312 +        if (metaDataBlockStart == -1) {
313 +          i = CaseInsensitiveFind(line, "<MetaData>");
314 +          if (i != string::npos) {
315 +            metaDataBlockStart = lineNo;
316 +            mdOffset = mdFile_.tellg();
317 +          }
318 +        } else {
319 +          i = CaseInsensitiveFind(line, "</MetaData>");
320 +          if (i != string::npos) {
321 +            metaDataBlockEnd = lineNo;
322 +          }
323 +        }
324 +      }
325 +
326 +      if (metaDataBlockStart == -1) {
327 +        sprintf(painCave.errMsg,
328 +                "SimCreator: File: %s did not contain a <MetaData> tag!\n",
329 +                mdFileName.c_str());
330 +        painCave.isFatal = 1;
331 +        simError();
332 +      }
333 +      if (metaDataBlockEnd == -1) {
334 +        sprintf(painCave.errMsg,
335 +                "SimCreator: File: %s did not contain a closed MetaData block!\n",
336 +                mdFileName.c_str());
337 +        painCave.isFatal = 1;
338 +        simError();
339 +      }
340 +        
341 +      mdFile_.clear();
342 +      mdFile_.seekg(0);
343 +      mdFile_.seekg(mdOffset);
344 +
345 +      mdRawData.clear();
346 +
347 +      for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) {
348 +        mdFile_.getline(buffer, bufferSize);
349 +        mdRawData += buffer;
350 +        mdRawData += "\n";
351 +      }
352 +
353 +      mdFile_.close();
354 +
355 + #ifdef IS_MPI
356 +    }
357 + #endif
358 +
359 +    std::stringstream rawMetaDataStream(mdRawData);
360 +
361      //parse meta-data file
362 <    Globals* simParams = parseFile(mdFileName);
362 >    Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion,
363 >                                   metaDataBlockStart + 1);
364      
365      //create the force field
366 <    ForceField * ff = ForceFieldFactory::getInstance()
367 <      ->createForceField(simParams->getForceField());
245 <    
366 >    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField());
367 >
368      if (ff == NULL) {
369        sprintf(painCave.errMsg,
370                "ForceField Factory can not create %s force field\n",
# Line 275 | Line 397 | Globals* SimCreator::parseFile(const std::string mdFil
397      }
398      
399      ff->parse(forcefieldFileName);
278    ff->setFortranForceOptions();
400      //create SimInfo
401      SimInfo * info = new SimInfo(ff, simParams);
402 +
403 +    info->setRawMetaData(mdRawData);
404      
405      //gather parameters (SimCreator only retrieves part of the
406      //parameters)
# Line 291 | Line 414 | Globals* SimCreator::parseFile(const std::string mdFil
414      //create the molecules
415      createMolecules(info);
416      
294    
417      //allocate memory for DataStorage(circular reference, need to
418      //break it)
419      info->setSnapshotManager(new SimSnapshotManager(info));
# Line 303 | Line 425 | Globals* SimCreator::parseFile(const std::string mdFil
425      //responsibility to LocalIndexManager.
426      setGlobalIndex(info);
427      
428 <    //Although addExcludePairs is called inside SimInfo's addMolecule
428 >    //Although addInteractionPairs is called inside SimInfo's addMolecule
429      //method, at that point atoms don't have the global index yet
430      //(their global index are all initialized to -1).  Therefore we
431 <    //have to call addExcludePairs explicitly here. A way to work
431 >    //have to call addInteractionPairs explicitly here. A way to work
432      //around is that we can determine the beginning global indices of
433      //atoms before they get created.
434      SimInfo::MoleculeIterator mi;
435      Molecule* mol;
436      for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
437 <      info->addExcludePairs(mol);
437 >      info->addInteractionPairs(mol);
438      }
439      
440      if (loadInitCoords)
441 <      loadCoordinates(info);    
320 <    
441 >      loadCoordinates(info, mdFileName);    
442      return info;
443    }
444    
# Line 352 | Line 473 | Globals* SimCreator::parseFile(const std::string mdFil
473    
474   #ifdef IS_MPI
475    void SimCreator::divideMolecules(SimInfo *info) {
476 <    double numerator;
477 <    double denominator;
478 <    double precast;
479 <    double x;
480 <    double y;
481 <    double a;
476 >    RealType numerator;
477 >    RealType denominator;
478 >    RealType precast;
479 >    RealType x;
480 >    RealType y;
481 >    RealType a;
482      int old_atoms;
483      int add_atoms;
484      int new_atoms;
# Line 381 | Line 502 | Globals* SimCreator::parseFile(const std::string mdFil
502                "\tthe number of molecules.  This will not result in a \n"
503                "\tusable division of atoms for force decomposition.\n"
504                "\tEither try a smaller number of processors, or run the\n"
505 <              "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
505 >              "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols);
506        
507        painCave.isFatal = 1;
508        simError();
# Line 469 | Line 590 | Globals* SimCreator::parseFile(const std::string mdFil
590            //           Pacc(x) = exp(- a * x)
591            // where a = penalty / (average atoms per molecule)
592            
593 <          x = (double)(new_atoms - nTarget);
593 >          x = (RealType)(new_atoms - nTarget);
594            y = myRandom->rand();
595            
596            if (y < exp(- a * x)) {
# Line 499 | Line 620 | Globals* SimCreator::parseFile(const std::string mdFil
620      info->setMolToProcMap(molToProcMap);
621      sprintf(checkPointMsg,
622              "Successfully divided the molecules among the processors.\n");
623 <    MPIcheckPoint();
623 >    errorCheckPoint();
624    }
625    
626   #endif
# Line 516 | Line 637 | Globals* SimCreator::parseFile(const std::string mdFil
637   #endif
638          
639          stampId = info->getMoleculeStampId(i);
640 <        Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
641 <                                                   stampId, i, info->getLocalIndexManager());
640 >        Molecule * mol = molCreator.createMolecule(info->getForceField(),
641 >                                                   info->getMoleculeStamp(stampId),
642 >                                                   stampId, i,
643 >                                                   info->getLocalIndexManager());
644          
645          info->addMolecule(mol);
646          
# Line 535 | Line 658 | Globals* SimCreator::parseFile(const std::string mdFil
658      Molecule::AtomIterator ai;
659      Molecule::RigidBodyIterator ri;
660      Molecule::CutoffGroupIterator ci;
661 +    Molecule::IntegrableObjectIterator  ioi;
662      Molecule * mol;
663      Atom * atom;
664      RigidBody * rb;
# Line 543 | Line 667 | Globals* SimCreator::parseFile(const std::string mdFil
667      int beginRigidBodyIndex;
668      int beginCutoffGroupIndex;
669      int nGlobalAtoms = info->getNGlobalAtoms();
546    
547 #ifndef IS_MPI
548    
549    beginAtomIndex = 0;
550    beginRigidBodyIndex = 0;
551    beginCutoffGroupIndex = 0;
552    
553 #else
554    
555    int nproc;
556    int myNode;
557    
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);
670      
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    
671      beginAtomIndex = 0;
672      beginRigidBodyIndex = 0;
673      beginCutoffGroupIndex = 0;
674 <    
675 <    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)) {
674 >
675 >    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
676        
677 <      //local index(index in DataStorge) of atom is important
678 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
679 <        atom->setGlobalIndex(beginAtomIndex++);
677 > #ifdef IS_MPI      
678 >      if (info->getMolToProc(i) == worldRank) {
679 > #endif        
680 >        // stuff to do if I own this molecule
681 >        mol = info->getMoleculeByGlobalIndex(i);
682 >
683 >        //local index(index in DataStorge) of atom is important
684 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
685 >          atom->setGlobalIndex(beginAtomIndex++);
686 >        }
687 >        
688 >        for(rb = mol->beginRigidBody(ri); rb != NULL;
689 >            rb = mol->nextRigidBody(ri)) {
690 >          rb->setGlobalIndex(beginRigidBodyIndex++);
691 >        }
692 >        
693 >        //local index of cutoff group is trivial, it only depends on
694 >        //the order of travesing
695 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
696 >            cg = mol->nextCutoffGroup(ci)) {
697 >          cg->setGlobalIndex(beginCutoffGroupIndex++);
698 >        }        
699 >        
700 > #ifdef IS_MPI        
701 >      }  else {
702 >
703 >        // stuff to do if I don't own this molecule
704 >        
705 >        int stampId = info->getMoleculeStampId(i);
706 >        MoleculeStamp* stamp = info->getMoleculeStamp(stampId);
707 >
708 >        beginAtomIndex += stamp->getNAtoms();
709 >        beginRigidBodyIndex += stamp->getNRigidBodies();
710 >        beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms();
711        }
712 <      
713 <      for(rb = mol->beginRigidBody(ri); rb != NULL;
714 <          rb = mol->nextRigidBody(ri)) {
715 <        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 <    
712 > #endif          
713 >
714 >    } //end for(int i=0)  
715 >
716      //fill globalGroupMembership
717      std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
718      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {        
# Line 623 | Line 724 | Globals* SimCreator::parseFile(const std::string mdFil
724          
725        }      
726      }
727 <    
727 >  
728   #ifdef IS_MPI    
729      // Since the globalGroupMembership has been zero filled and we've only
730      // poked values into the atoms we know, we can do an Allreduce
731      // to get the full globalGroupMembership array (We think).
732      // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
733      // docs said we could.
734 <    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
734 >    std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0);
735      MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
736                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
737      info->setGlobalGroupMembership(tmpGroupMembership);
# Line 642 | Line 743 | Globals* SimCreator::parseFile(const std::string mdFil
743      std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
744      
745      for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
645      
746        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
747          globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
748        }
749      }
750      
751   #ifdef IS_MPI
752 <    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
752 >    std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0);
753      
754      MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
755                    MPI_INT, MPI_SUM, MPI_COMM_WORLD);
# Line 658 | Line 758 | Globals* SimCreator::parseFile(const std::string mdFil
758   #else
759      info->setGlobalMolMembership(globalMolMembership);
760   #endif
761 +
762 +    // nIOPerMol holds the number of integrable objects per molecule
763 +    // here the molecules are listed by their global indices.
764 +
765 +    std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0);
766 +    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
767 +      nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects();      
768 +    }
769      
770 + #ifdef IS_MPI
771 +    std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0);
772 +    MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0],
773 +                  info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD);
774 + #else
775 +    std::vector<int> numIntegrableObjectsPerMol = nIOPerMol;
776 + #endif    
777 +
778 +    std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules());
779 +    
780 +    int startingIndex = 0;
781 +    for (int i = 0; i < info->getNGlobalMolecules(); i++) {
782 +      startingIOIndexForMol[i] = startingIndex;
783 +      startingIndex += numIntegrableObjectsPerMol[i];
784 +    }
785 +    
786 +    std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL);
787 +    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
788 +      int myGlobalIndex = mol->getGlobalIndex();
789 +      int globalIO = startingIOIndexForMol[myGlobalIndex];
790 +      for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL;
791 +           integrableObject = mol->nextIntegrableObject(ioi)) {
792 +        integrableObject->setGlobalIntegrableObjectIndex(globalIO);
793 +        IOIndexToIntegrableObject[globalIO] = integrableObject;
794 +        globalIO++;
795 +      }
796 +    }
797 +      
798 +    info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject);
799 +    
800    }
801    
802 <  void SimCreator::loadCoordinates(SimInfo* info) {
802 >  void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) {
803      Globals* simParams;
804 +
805      simParams = info->getSimParams();
806      
807 <    if (!simParams->haveInitialConfig()) {
669 <      sprintf(painCave.errMsg,
670 <              "Cannot intialize a simulation without an initial configuration file.\n");
671 <      painCave.isFatal = 1;;
672 <      simError();
673 <    }
674 <    
675 <    DumpReader reader(info, simParams->getInitialConfig());
807 >    DumpReader reader(info, mdFileName);
808      int nframes = reader.getNFrames();
809 <    
809 >
810      if (nframes > 0) {
811        reader.readFrame(nframes - 1);
812      } else {
813        //invalid initial coordinate file
814        sprintf(painCave.errMsg,
815                "Initial configuration file %s should at least contain one frame\n",
816 <              simParams->getInitialConfig().c_str());
816 >              mdFileName.c_str());
817        painCave.isFatal = 1;
818        simError();
819      }
688    
820      //copy the current snapshot to previous snapshot
821      info->getSnapshotManager()->advance();
822    }
823    
824 < } //end namespace oopse
824 > } //end namespace OpenMD
825  
826  

Comparing:
trunk/src/brains/SimCreator.cpp (property svn:keywords), Revision 945 by gezelter, Tue Apr 25 02:09:01 2006 UTC vs.
branches/development/src/brains/SimCreator.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines