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

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 1627 by gezelter, Tue Sep 13 22:05:04 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines