1 |
gezelter |
403 |
/* |
2 |
|
|
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
|
|
* |
4 |
|
|
* The University of Notre Dame grants you ("Licensee") a |
5 |
|
|
* non-exclusive, royalty free, license to use, modify and |
6 |
|
|
* redistribute this software in source and binary code form, provided |
7 |
|
|
* that the following conditions are met: |
8 |
|
|
* |
9 |
gezelter |
1390 |
* 1. Redistributions of source code must retain the above copyright |
10 |
gezelter |
403 |
* notice, this list of conditions and the following disclaimer. |
11 |
|
|
* |
12 |
gezelter |
1390 |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
gezelter |
403 |
* notice, this list of conditions and the following disclaimer in the |
14 |
|
|
* documentation and/or other materials provided with the |
15 |
|
|
* distribution. |
16 |
|
|
* |
17 |
|
|
* This software is provided "AS IS," without a warranty of any |
18 |
|
|
* kind. All express or implied conditions, representations and |
19 |
|
|
* warranties, including any implied warranty of merchantability, |
20 |
|
|
* fitness for a particular purpose or non-infringement, are hereby |
21 |
|
|
* excluded. The University of Notre Dame and its licensors shall not |
22 |
|
|
* be liable for any damages suffered by licensee as a result of |
23 |
|
|
* using, modifying or distributing the software or its |
24 |
|
|
* derivatives. In no event will the University of Notre Dame or its |
25 |
|
|
* licensors be liable for any lost revenue, profit or data, or for |
26 |
|
|
* direct, indirect, special, consequential, incidental or punitive |
27 |
|
|
* damages, however caused and regardless of the theory of liability, |
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 |
gezelter |
1390 |
* |
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 |
gezelter |
1665 |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
|
|
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
gezelter |
403 |
*/ |
42 |
|
|
|
43 |
|
|
/** |
44 |
|
|
* @file SimCreator.cpp |
45 |
|
|
* @author tlin |
46 |
|
|
* @date 11/03/2004 |
47 |
|
|
* @time 13:51am |
48 |
|
|
* @version 1.0 |
49 |
|
|
*/ |
50 |
tim |
845 |
#include <exception> |
51 |
tim |
770 |
#include <iostream> |
52 |
|
|
#include <sstream> |
53 |
|
|
#include <string> |
54 |
|
|
|
55 |
gezelter |
403 |
#include "brains/MoleculeCreator.hpp" |
56 |
|
|
#include "brains/SimCreator.hpp" |
57 |
|
|
#include "brains/SimSnapshotManager.hpp" |
58 |
|
|
#include "io/DumpReader.hpp" |
59 |
|
|
#include "UseTheForce/ForceFieldFactory.hpp" |
60 |
|
|
#include "utils/simError.h" |
61 |
|
|
#include "utils/StringUtils.hpp" |
62 |
|
|
#include "math/SeqRandNumGen.hpp" |
63 |
tim |
770 |
#include "mdParser/MDLexer.hpp" |
64 |
|
|
#include "mdParser/MDParser.hpp" |
65 |
|
|
#include "mdParser/MDTreeParser.hpp" |
66 |
|
|
#include "mdParser/SimplePreprocessor.hpp" |
67 |
tim |
816 |
#include "antlr/ANTLRException.hpp" |
68 |
|
|
#include "antlr/TokenStreamRecognitionException.hpp" |
69 |
|
|
#include "antlr/TokenStreamIOException.hpp" |
70 |
|
|
#include "antlr/TokenStreamException.hpp" |
71 |
|
|
#include "antlr/RecognitionException.hpp" |
72 |
|
|
#include "antlr/CharStreamException.hpp" |
73 |
tim |
770 |
|
74 |
tim |
816 |
#include "antlr/MismatchedCharException.hpp" |
75 |
|
|
#include "antlr/MismatchedTokenException.hpp" |
76 |
|
|
#include "antlr/NoViableAltForCharException.hpp" |
77 |
|
|
#include "antlr/NoViableAltException.hpp" |
78 |
tim |
770 |
|
79 |
gezelter |
403 |
#ifdef IS_MPI |
80 |
gezelter |
1627 |
#include "mpi.h" |
81 |
gezelter |
403 |
#include "math/ParallelRandNumGen.hpp" |
82 |
|
|
#endif |
83 |
|
|
|
84 |
gezelter |
1390 |
namespace OpenMD { |
85 |
chrisfen |
417 |
|
86 |
gezelter |
1613 |
Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){ |
87 |
tim |
1024 |
Globals* simParams = NULL; |
88 |
|
|
try { |
89 |
tim |
770 |
|
90 |
tim |
1024 |
// Create a preprocessor that preprocesses md file into an ostringstream |
91 |
|
|
std::stringstream ppStream; |
92 |
tim |
770 |
#ifdef IS_MPI |
93 |
tim |
1024 |
int streamSize; |
94 |
|
|
const int masterNode = 0; |
95 |
|
|
int commStatus; |
96 |
|
|
if (worldRank == masterNode) { |
97 |
gezelter |
1613 |
commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
98 |
|
|
#endif |
99 |
tim |
1024 |
SimplePreprocessor preprocessor; |
100 |
|
|
preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); |
101 |
tim |
770 |
|
102 |
|
|
#ifdef IS_MPI |
103 |
tim |
1024 |
//brocasting the stream size |
104 |
|
|
streamSize = ppStream.str().size() +1; |
105 |
|
|
commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); |
106 |
tim |
770 |
|
107 |
tim |
1024 |
commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
108 |
tim |
770 |
|
109 |
|
|
|
110 |
tim |
1024 |
} else { |
111 |
gezelter |
1613 |
|
112 |
|
|
commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
113 |
|
|
|
114 |
tim |
1024 |
//get stream size |
115 |
|
|
commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); |
116 |
gezelter |
1313 |
|
117 |
tim |
1024 |
char* buf = new char[streamSize]; |
118 |
|
|
assert(buf); |
119 |
tim |
770 |
|
120 |
tim |
1024 |
//receive file content |
121 |
|
|
commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
122 |
tim |
770 |
|
123 |
tim |
1024 |
ppStream.str(buf); |
124 |
gezelter |
1313 |
delete [] buf; |
125 |
tim |
770 |
|
126 |
tim |
1024 |
} |
127 |
tim |
770 |
#endif |
128 |
tim |
1024 |
// Create a scanner that reads from the input stream |
129 |
|
|
MDLexer lexer(ppStream); |
130 |
|
|
lexer.setFilename(filename); |
131 |
|
|
lexer.initDeferredLineCount(); |
132 |
chrisfen |
417 |
|
133 |
tim |
1024 |
// Create a parser that reads from the scanner |
134 |
|
|
MDParser parser(lexer); |
135 |
|
|
parser.setFilename(filename); |
136 |
tim |
770 |
|
137 |
tim |
1024 |
// Create an observer that synchorizes file name change |
138 |
|
|
FilenameObserver observer; |
139 |
|
|
observer.setLexer(&lexer); |
140 |
|
|
observer.setParser(&parser); |
141 |
|
|
lexer.setObserver(&observer); |
142 |
chrisfen |
417 |
|
143 |
tim |
1024 |
antlr::ASTFactory factory; |
144 |
|
|
parser.initializeASTFactory(factory); |
145 |
|
|
parser.setASTFactory(&factory); |
146 |
|
|
parser.mdfile(); |
147 |
tim |
770 |
|
148 |
tim |
1024 |
// 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 |
tim |
845 |
|
155 |
tim |
816 |
|
156 |
tim |
1024 |
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 |
tim |
845 |
|
185 |
tim |
1024 |
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 |
tim |
845 |
|
193 |
tim |
1024 |
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 |
tim |
845 |
|
201 |
tim |
1024 |
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 |
gezelter |
1390 |
catch (OpenMDException& e) { |
223 |
tim |
1024 |
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 |
tim |
770 |
|
237 |
gezelter |
1613 |
simParams->setMDfileVersion(mdFileVersion); |
238 |
tim |
1024 |
return simParams; |
239 |
gezelter |
403 |
} |
240 |
chrisfen |
417 |
|
241 |
chrisfen |
514 |
SimInfo* SimCreator::createSim(const std::string & mdFileName, |
242 |
|
|
bool loadInitCoords) { |
243 |
gezelter |
1469 |
|
244 |
tim |
1024 |
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 |
gezelter |
1613 |
int mdFileVersion; |
253 |
tim |
1024 |
|
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 |
gezelter |
1390 |
i = CaseInsensitiveFind(line, "<OpenMD"); |
273 |
gezelter |
1287 |
if (static_cast<size_t>(i) == string::npos) { |
274 |
gezelter |
1390 |
// 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 |
tim |
1024 |
sprintf(painCave.errMsg, |
281 |
gezelter |
1390 |
"SimCreator: File: %s is not a valid OpenMD file!\n", |
282 |
tim |
1024 |
mdFileName.c_str()); |
283 |
|
|
painCave.isFatal = 1; |
284 |
|
|
simError(); |
285 |
|
|
} |
286 |
gezelter |
1613 |
|
287 |
|
|
// found the correct opening string, now try to get the file |
288 |
|
|
// format version number. |
289 |
tim |
1024 |
|
290 |
gezelter |
1613 |
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 |
tim |
1024 |
//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 |
gezelter |
403 |
//parse meta-data file |
362 |
gezelter |
1613 |
Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion, |
363 |
|
|
metaDataBlockStart + 1); |
364 |
chrisfen |
417 |
|
365 |
gezelter |
403 |
//create the force field |
366 |
gezelter |
1277 |
ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField()); |
367 |
|
|
|
368 |
gezelter |
403 |
if (ff == NULL) { |
369 |
chrisfen |
514 |
sprintf(painCave.errMsg, |
370 |
|
|
"ForceField Factory can not create %s force field\n", |
371 |
tim |
665 |
simParams->getForceField().c_str()); |
372 |
gezelter |
403 |
painCave.isFatal = 1; |
373 |
|
|
simError(); |
374 |
|
|
} |
375 |
chrisfen |
417 |
|
376 |
gezelter |
403 |
if (simParams->haveForceFieldFileName()) { |
377 |
|
|
ff->setForceFieldFileName(simParams->getForceFieldFileName()); |
378 |
|
|
} |
379 |
|
|
|
380 |
|
|
std::string forcefieldFileName; |
381 |
|
|
forcefieldFileName = ff->getForceFieldFileName(); |
382 |
chrisfen |
417 |
|
383 |
gezelter |
403 |
if (simParams->haveForceFieldVariant()) { |
384 |
|
|
//If the force field has variant, the variant force field name will be |
385 |
|
|
//Base.variant.frc. For exampel EAM.u6.frc |
386 |
chrisfen |
417 |
|
387 |
gezelter |
403 |
std::string variant = simParams->getForceFieldVariant(); |
388 |
chrisfen |
417 |
|
389 |
gezelter |
403 |
std::string::size_type pos = forcefieldFileName.rfind(".frc"); |
390 |
|
|
variant = "." + variant; |
391 |
|
|
if (pos != std::string::npos) { |
392 |
|
|
forcefieldFileName.insert(pos, variant); |
393 |
|
|
} else { |
394 |
|
|
//If the default force field file name does not containt .frc suffix, just append the .variant |
395 |
|
|
forcefieldFileName.append(variant); |
396 |
|
|
} |
397 |
|
|
} |
398 |
|
|
|
399 |
|
|
ff->parse(forcefieldFileName); |
400 |
|
|
//create SimInfo |
401 |
tim |
770 |
SimInfo * info = new SimInfo(ff, simParams); |
402 |
tim |
1024 |
|
403 |
|
|
info->setRawMetaData(mdRawData); |
404 |
tim |
490 |
|
405 |
gezelter |
945 |
//gather parameters (SimCreator only retrieves part of the |
406 |
|
|
//parameters) |
407 |
gezelter |
403 |
gatherParameters(info, mdFileName); |
408 |
chrisfen |
417 |
|
409 |
gezelter |
403 |
//divide the molecules and determine the global index of molecules |
410 |
|
|
#ifdef IS_MPI |
411 |
|
|
divideMolecules(info); |
412 |
|
|
#endif |
413 |
chrisfen |
417 |
|
414 |
gezelter |
403 |
//create the molecules |
415 |
|
|
createMolecules(info); |
416 |
chrisfen |
417 |
|
417 |
gezelter |
945 |
//allocate memory for DataStorage(circular reference, need to |
418 |
|
|
//break it) |
419 |
gezelter |
403 |
info->setSnapshotManager(new SimSnapshotManager(info)); |
420 |
|
|
|
421 |
gezelter |
945 |
//set the global index of atoms, rigidbodies and cutoffgroups |
422 |
|
|
//(only need to be set once, the global index will never change |
423 |
|
|
//again). Local indices of atoms and rigidbodies are already set |
424 |
|
|
//by MoleculeCreator class which actually delegates the |
425 |
|
|
//responsibility to LocalIndexManager. |
426 |
gezelter |
403 |
setGlobalIndex(info); |
427 |
chrisfen |
417 |
|
428 |
gezelter |
1287 |
//Although addInteractionPairs is called inside SimInfo's addMolecule |
429 |
gezelter |
945 |
//method, at that point atoms don't have the global index yet |
430 |
|
|
//(their global index are all initialized to -1). Therefore we |
431 |
gezelter |
1287 |
//have to call addInteractionPairs explicitly here. A way to work |
432 |
gezelter |
945 |
//around is that we can determine the beginning global indices of |
433 |
|
|
//atoms before they get created. |
434 |
gezelter |
403 |
SimInfo::MoleculeIterator mi; |
435 |
|
|
Molecule* mol; |
436 |
|
|
for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
437 |
gezelter |
1287 |
info->addInteractionPairs(mol); |
438 |
gezelter |
403 |
} |
439 |
|
|
|
440 |
|
|
if (loadInitCoords) |
441 |
tim |
1024 |
loadCoordinates(info, mdFileName); |
442 |
gezelter |
403 |
return info; |
443 |
|
|
} |
444 |
chrisfen |
417 |
|
445 |
gezelter |
403 |
void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) { |
446 |
chrisfen |
417 |
|
447 |
tim |
749 |
//figure out the output file names |
448 |
gezelter |
403 |
std::string prefix; |
449 |
chrisfen |
417 |
|
450 |
gezelter |
403 |
#ifdef IS_MPI |
451 |
chrisfen |
417 |
|
452 |
gezelter |
403 |
if (worldRank == 0) { |
453 |
|
|
#endif // is_mpi |
454 |
|
|
Globals * simParams = info->getSimParams(); |
455 |
|
|
if (simParams->haveFinalConfig()) { |
456 |
|
|
prefix = getPrefix(simParams->getFinalConfig()); |
457 |
|
|
} else { |
458 |
|
|
prefix = getPrefix(mdfile); |
459 |
|
|
} |
460 |
chrisfen |
417 |
|
461 |
gezelter |
403 |
info->setFinalConfigFileName(prefix + ".eor"); |
462 |
|
|
info->setDumpFileName(prefix + ".dump"); |
463 |
|
|
info->setStatFileName(prefix + ".stat"); |
464 |
chrisfen |
417 |
info->setRestFileName(prefix + ".zang"); |
465 |
|
|
|
466 |
gezelter |
403 |
#ifdef IS_MPI |
467 |
chrisfen |
417 |
|
468 |
gezelter |
403 |
} |
469 |
chrisfen |
417 |
|
470 |
gezelter |
403 |
#endif |
471 |
chrisfen |
417 |
|
472 |
gezelter |
403 |
} |
473 |
chrisfen |
417 |
|
474 |
gezelter |
403 |
#ifdef IS_MPI |
475 |
|
|
void SimCreator::divideMolecules(SimInfo *info) { |
476 |
tim |
963 |
RealType numerator; |
477 |
|
|
RealType denominator; |
478 |
|
|
RealType precast; |
479 |
|
|
RealType x; |
480 |
|
|
RealType y; |
481 |
|
|
RealType a; |
482 |
gezelter |
403 |
int old_atoms; |
483 |
|
|
int add_atoms; |
484 |
|
|
int new_atoms; |
485 |
|
|
int nTarget; |
486 |
|
|
int done; |
487 |
|
|
int i; |
488 |
|
|
int j; |
489 |
|
|
int loops; |
490 |
|
|
int which_proc; |
491 |
|
|
int nProcessors; |
492 |
|
|
std::vector<int> atomsPerProc; |
493 |
|
|
int nGlobalMols = info->getNGlobalMolecules(); |
494 |
|
|
std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition: |
495 |
|
|
|
496 |
|
|
MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); |
497 |
chrisfen |
417 |
|
498 |
gezelter |
403 |
if (nProcessors > nGlobalMols) { |
499 |
|
|
sprintf(painCave.errMsg, |
500 |
|
|
"nProcessors (%d) > nMol (%d)\n" |
501 |
|
|
"\tThe number of processors is larger than\n" |
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 |
gezelter |
1390 |
"\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols); |
506 |
chrisfen |
417 |
|
507 |
gezelter |
403 |
painCave.isFatal = 1; |
508 |
|
|
simError(); |
509 |
|
|
} |
510 |
chrisfen |
417 |
|
511 |
gezelter |
403 |
int seedValue; |
512 |
|
|
Globals * simParams = info->getSimParams(); |
513 |
|
|
SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator |
514 |
|
|
if (simParams->haveSeed()) { |
515 |
|
|
seedValue = simParams->getSeed(); |
516 |
|
|
myRandom = new SeqRandNumGen(seedValue); |
517 |
|
|
}else { |
518 |
|
|
myRandom = new SeqRandNumGen(); |
519 |
|
|
} |
520 |
chrisfen |
417 |
|
521 |
|
|
|
522 |
gezelter |
403 |
a = 3.0 * nGlobalMols / info->getNGlobalAtoms(); |
523 |
chrisfen |
417 |
|
524 |
gezelter |
403 |
//initialize atomsPerProc |
525 |
|
|
atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0); |
526 |
chrisfen |
417 |
|
527 |
gezelter |
403 |
if (worldRank == 0) { |
528 |
|
|
numerator = info->getNGlobalAtoms(); |
529 |
|
|
denominator = nProcessors; |
530 |
|
|
precast = numerator / denominator; |
531 |
|
|
nTarget = (int)(precast + 0.5); |
532 |
chrisfen |
417 |
|
533 |
gezelter |
403 |
for(i = 0; i < nGlobalMols; i++) { |
534 |
|
|
done = 0; |
535 |
|
|
loops = 0; |
536 |
chrisfen |
417 |
|
537 |
gezelter |
403 |
while (!done) { |
538 |
|
|
loops++; |
539 |
chrisfen |
417 |
|
540 |
gezelter |
403 |
// Pick a processor at random |
541 |
chrisfen |
417 |
|
542 |
gezelter |
403 |
which_proc = (int) (myRandom->rand() * nProcessors); |
543 |
chrisfen |
417 |
|
544 |
gezelter |
403 |
//get the molecule stamp first |
545 |
|
|
int stampId = info->getMoleculeStampId(i); |
546 |
|
|
MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId); |
547 |
chrisfen |
417 |
|
548 |
gezelter |
403 |
// How many atoms does this processor have so far? |
549 |
|
|
old_atoms = atomsPerProc[which_proc]; |
550 |
|
|
add_atoms = moleculeStamp->getNAtoms(); |
551 |
|
|
new_atoms = old_atoms + add_atoms; |
552 |
chrisfen |
417 |
|
553 |
gezelter |
403 |
// If we've been through this loop too many times, we need |
554 |
|
|
// to just give up and assign the molecule to this processor |
555 |
|
|
// and be done with it. |
556 |
chrisfen |
417 |
|
557 |
gezelter |
403 |
if (loops > 100) { |
558 |
|
|
sprintf(painCave.errMsg, |
559 |
|
|
"I've tried 100 times to assign molecule %d to a " |
560 |
|
|
" processor, but can't find a good spot.\n" |
561 |
|
|
"I'm assigning it at random to processor %d.\n", |
562 |
|
|
i, which_proc); |
563 |
chrisfen |
417 |
|
564 |
gezelter |
403 |
painCave.isFatal = 0; |
565 |
|
|
simError(); |
566 |
chrisfen |
417 |
|
567 |
gezelter |
403 |
molToProcMap[i] = which_proc; |
568 |
|
|
atomsPerProc[which_proc] += add_atoms; |
569 |
chrisfen |
417 |
|
570 |
gezelter |
403 |
done = 1; |
571 |
|
|
continue; |
572 |
|
|
} |
573 |
chrisfen |
417 |
|
574 |
gezelter |
403 |
// If we can add this molecule to this processor without sending |
575 |
|
|
// it above nTarget, then go ahead and do it: |
576 |
chrisfen |
417 |
|
577 |
gezelter |
403 |
if (new_atoms <= nTarget) { |
578 |
|
|
molToProcMap[i] = which_proc; |
579 |
|
|
atomsPerProc[which_proc] += add_atoms; |
580 |
chrisfen |
417 |
|
581 |
gezelter |
403 |
done = 1; |
582 |
|
|
continue; |
583 |
|
|
} |
584 |
chrisfen |
417 |
|
585 |
gezelter |
403 |
// The only situation left is when new_atoms > nTarget. We |
586 |
|
|
// want to accept this with some probability that dies off the |
587 |
|
|
// farther we are from nTarget |
588 |
chrisfen |
417 |
|
589 |
gezelter |
403 |
// roughly: x = new_atoms - nTarget |
590 |
|
|
// Pacc(x) = exp(- a * x) |
591 |
|
|
// where a = penalty / (average atoms per molecule) |
592 |
chrisfen |
417 |
|
593 |
tim |
963 |
x = (RealType)(new_atoms - nTarget); |
594 |
gezelter |
403 |
y = myRandom->rand(); |
595 |
chrisfen |
417 |
|
596 |
gezelter |
403 |
if (y < exp(- a * x)) { |
597 |
|
|
molToProcMap[i] = which_proc; |
598 |
|
|
atomsPerProc[which_proc] += add_atoms; |
599 |
chrisfen |
417 |
|
600 |
gezelter |
403 |
done = 1; |
601 |
|
|
continue; |
602 |
|
|
} else { |
603 |
|
|
continue; |
604 |
|
|
} |
605 |
|
|
} |
606 |
|
|
} |
607 |
chrisfen |
417 |
|
608 |
gezelter |
403 |
delete myRandom; |
609 |
chrisfen |
417 |
|
610 |
gezelter |
403 |
// Spray out this nonsense to all other processors: |
611 |
chrisfen |
417 |
|
612 |
gezelter |
403 |
MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); |
613 |
|
|
} else { |
614 |
chrisfen |
417 |
|
615 |
gezelter |
403 |
// Listen to your marching orders from processor 0: |
616 |
chrisfen |
417 |
|
617 |
gezelter |
403 |
MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); |
618 |
|
|
} |
619 |
chrisfen |
417 |
|
620 |
gezelter |
403 |
info->setMolToProcMap(molToProcMap); |
621 |
|
|
sprintf(checkPointMsg, |
622 |
|
|
"Successfully divided the molecules among the processors.\n"); |
623 |
gezelter |
1241 |
errorCheckPoint(); |
624 |
gezelter |
403 |
} |
625 |
chrisfen |
417 |
|
626 |
gezelter |
403 |
#endif |
627 |
chrisfen |
417 |
|
628 |
gezelter |
403 |
void SimCreator::createMolecules(SimInfo *info) { |
629 |
|
|
MoleculeCreator molCreator; |
630 |
|
|
int stampId; |
631 |
chrisfen |
417 |
|
632 |
gezelter |
403 |
for(int i = 0; i < info->getNGlobalMolecules(); i++) { |
633 |
chrisfen |
417 |
|
634 |
gezelter |
403 |
#ifdef IS_MPI |
635 |
chrisfen |
417 |
|
636 |
gezelter |
403 |
if (info->getMolToProc(i) == worldRank) { |
637 |
|
|
#endif |
638 |
chrisfen |
417 |
|
639 |
gezelter |
403 |
stampId = info->getMoleculeStampId(i); |
640 |
gezelter |
1600 |
Molecule * mol = molCreator.createMolecule(info->getForceField(), |
641 |
|
|
info->getMoleculeStamp(stampId), |
642 |
|
|
stampId, i, |
643 |
|
|
info->getLocalIndexManager()); |
644 |
chrisfen |
417 |
|
645 |
gezelter |
403 |
info->addMolecule(mol); |
646 |
chrisfen |
417 |
|
647 |
gezelter |
403 |
#ifdef IS_MPI |
648 |
chrisfen |
417 |
|
649 |
gezelter |
403 |
} |
650 |
chrisfen |
417 |
|
651 |
gezelter |
403 |
#endif |
652 |
chrisfen |
417 |
|
653 |
gezelter |
403 |
} //end for(int i=0) |
654 |
|
|
} |
655 |
chrisfen |
417 |
|
656 |
gezelter |
403 |
void SimCreator::setGlobalIndex(SimInfo *info) { |
657 |
|
|
SimInfo::MoleculeIterator mi; |
658 |
|
|
Molecule::AtomIterator ai; |
659 |
|
|
Molecule::RigidBodyIterator ri; |
660 |
|
|
Molecule::CutoffGroupIterator ci; |
661 |
tim |
1024 |
Molecule::IntegrableObjectIterator ioi; |
662 |
gezelter |
403 |
Molecule * mol; |
663 |
|
|
Atom * atom; |
664 |
|
|
RigidBody * rb; |
665 |
|
|
CutoffGroup * cg; |
666 |
|
|
int beginAtomIndex; |
667 |
|
|
int beginRigidBodyIndex; |
668 |
|
|
int beginCutoffGroupIndex; |
669 |
|
|
int nGlobalAtoms = info->getNGlobalAtoms(); |
670 |
chrisfen |
417 |
|
671 |
gezelter |
403 |
beginAtomIndex = 0; |
672 |
|
|
beginRigidBodyIndex = 0; |
673 |
|
|
beginCutoffGroupIndex = 0; |
674 |
gezelter |
1600 |
|
675 |
|
|
for(int i = 0; i < info->getNGlobalMolecules(); i++) { |
676 |
chrisfen |
417 |
|
677 |
gezelter |
1600 |
#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 |
gezelter |
403 |
} |
712 |
gezelter |
1600 |
#endif |
713 |
|
|
|
714 |
|
|
} //end for(int i=0) |
715 |
|
|
|
716 |
gezelter |
403 |
//fill globalGroupMembership |
717 |
gezelter |
1600 |
std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); |
718 |
gezelter |
403 |
for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
719 |
|
|
for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
720 |
chrisfen |
417 |
|
721 |
gezelter |
403 |
for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { |
722 |
|
|
globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex(); |
723 |
|
|
} |
724 |
chrisfen |
417 |
|
725 |
gezelter |
403 |
} |
726 |
|
|
} |
727 |
gezelter |
1577 |
|
728 |
gezelter |
403 |
#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 |
gezelter |
1313 |
std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); |
735 |
gezelter |
403 |
MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, |
736 |
gezelter |
1600 |
MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
737 |
gezelter |
403 |
info->setGlobalGroupMembership(tmpGroupMembership); |
738 |
|
|
#else |
739 |
|
|
info->setGlobalGroupMembership(globalGroupMembership); |
740 |
|
|
#endif |
741 |
chrisfen |
417 |
|
742 |
gezelter |
403 |
//fill molMembership |
743 |
|
|
std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); |
744 |
|
|
|
745 |
|
|
for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
746 |
|
|
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
747 |
|
|
globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); |
748 |
|
|
} |
749 |
|
|
} |
750 |
chrisfen |
417 |
|
751 |
gezelter |
403 |
#ifdef IS_MPI |
752 |
gezelter |
1313 |
std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0); |
753 |
chrisfen |
417 |
|
754 |
gezelter |
403 |
MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, |
755 |
|
|
MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
756 |
|
|
|
757 |
|
|
info->setGlobalMolMembership(tmpMolMembership); |
758 |
|
|
#else |
759 |
|
|
info->setGlobalMolMembership(globalMolMembership); |
760 |
|
|
#endif |
761 |
tim |
1024 |
|
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 |
chrisfen |
417 |
|
770 |
tim |
1024 |
#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 |
gezelter |
1313 |
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 |
tim |
1024 |
int myGlobalIndex = mol->getGlobalIndex(); |
789 |
|
|
int globalIO = startingIOIndexForMol[myGlobalIndex]; |
790 |
|
|
for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL; |
791 |
|
|
integrableObject = mol->nextIntegrableObject(ioi)) { |
792 |
gezelter |
1313 |
integrableObject->setGlobalIntegrableObjectIndex(globalIO); |
793 |
|
|
IOIndexToIntegrableObject[globalIO] = integrableObject; |
794 |
|
|
globalIO++; |
795 |
tim |
1024 |
} |
796 |
|
|
} |
797 |
gezelter |
1597 |
|
798 |
gezelter |
1313 |
info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); |
799 |
|
|
|
800 |
gezelter |
403 |
} |
801 |
chrisfen |
417 |
|
802 |
tim |
1024 |
void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { |
803 |
gezelter |
403 |
Globals* simParams; |
804 |
gezelter |
1540 |
|
805 |
gezelter |
403 |
simParams = info->getSimParams(); |
806 |
|
|
|
807 |
tim |
1024 |
DumpReader reader(info, mdFileName); |
808 |
gezelter |
403 |
int nframes = reader.getNFrames(); |
809 |
gezelter |
1540 |
|
810 |
gezelter |
403 |
if (nframes > 0) { |
811 |
|
|
reader.readFrame(nframes - 1); |
812 |
|
|
} else { |
813 |
|
|
//invalid initial coordinate file |
814 |
chrisfen |
417 |
sprintf(painCave.errMsg, |
815 |
|
|
"Initial configuration file %s should at least contain one frame\n", |
816 |
tim |
1024 |
mdFileName.c_str()); |
817 |
gezelter |
403 |
painCave.isFatal = 1; |
818 |
|
|
simError(); |
819 |
|
|
} |
820 |
|
|
//copy the current snapshot to previous snapshot |
821 |
|
|
info->getSnapshotManager()->advance(); |
822 |
|
|
} |
823 |
chrisfen |
417 |
|
824 |
gezelter |
1390 |
} //end namespace OpenMD |
825 |
gezelter |
403 |
|
826 |
|
|
|