| 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. | 
| 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 |  | /** | 
| 46 |  | * @time 13:51am | 
| 47 |  | * @version 1.0 | 
| 48 |  | */ | 
| 49 | + | #include <exception> | 
| 50 | + | #include <iostream> | 
| 51 | + | #include <sstream> | 
| 52 | + | #include <string> | 
| 53 |  |  | 
| 54 |  | #include "brains/MoleculeCreator.hpp" | 
| 55 |  | #include "brains/SimCreator.hpp" | 
| 56 |  | #include "brains/SimSnapshotManager.hpp" | 
| 57 |  | #include "io/DumpReader.hpp" | 
| 54 | – | #include "io/parse_me.h" | 
| 58 |  | #include "UseTheForce/ForceFieldFactory.hpp" | 
| 59 |  | #include "utils/simError.h" | 
| 60 |  | #include "utils/StringUtils.hpp" | 
| 61 |  | #include "math/SeqRandNumGen.hpp" | 
| 62 | + | #include "mdParser/MDLexer.hpp" | 
| 63 | + | #include "mdParser/MDParser.hpp" | 
| 64 | + | #include "mdParser/MDTreeParser.hpp" | 
| 65 | + | #include "mdParser/SimplePreprocessor.hpp" | 
| 66 | + | #include "antlr/ANTLRException.hpp" | 
| 67 | + | #include "antlr/TokenStreamRecognitionException.hpp" | 
| 68 | + | #include "antlr/TokenStreamIOException.hpp" | 
| 69 | + | #include "antlr/TokenStreamException.hpp" | 
| 70 | + | #include "antlr/RecognitionException.hpp" | 
| 71 | + | #include "antlr/CharStreamException.hpp" | 
| 72 | + |  | 
| 73 | + | #include "antlr/MismatchedCharException.hpp" | 
| 74 | + | #include "antlr/MismatchedTokenException.hpp" | 
| 75 | + | #include "antlr/NoViableAltForCharException.hpp" | 
| 76 | + | #include "antlr/NoViableAltException.hpp" | 
| 77 | + |  | 
| 78 |  | #ifdef IS_MPI | 
| 60 | – | #include "io/mpiBASS.h" | 
| 79 |  | #include "math/ParallelRandNumGen.hpp" | 
| 80 |  | #endif | 
| 81 |  |  | 
| 82 | < | namespace oopse { | 
| 82 | > | namespace OpenMD { | 
| 83 |  |  | 
| 84 | < | void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, | 
| 85 | < | Globals* simParams){ | 
| 84 | > | Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){ | 
| 85 | > | Globals* simParams = NULL; | 
| 86 | > | try { | 
| 87 | > |  | 
| 88 | > | // Create a preprocessor that preprocesses md file into an ostringstream | 
| 89 | > | std::stringstream ppStream; | 
| 90 | > | #ifdef IS_MPI | 
| 91 | > | int streamSize; | 
| 92 | > | const int masterNode = 0; | 
| 93 | > | int commStatus; | 
| 94 | > | if (worldRank == masterNode) { | 
| 95 | > | #endif | 
| 96 | > |  | 
| 97 | > | SimplePreprocessor preprocessor; | 
| 98 | > | preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); | 
| 99 | > |  | 
| 100 | > | #ifdef IS_MPI | 
| 101 | > | //brocasting the stream size | 
| 102 | > | streamSize = ppStream.str().size() +1; | 
| 103 | > | commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); | 
| 104 | > |  | 
| 105 | > | commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 106 | > |  | 
| 107 | > |  | 
| 108 | > | } else { | 
| 109 | > | //get stream size | 
| 110 | > | commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); | 
| 111 | > |  | 
| 112 | > | char* buf = new char[streamSize]; | 
| 113 | > | assert(buf); | 
| 114 | > |  | 
| 115 | > | //receive file content | 
| 116 | > | commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); | 
| 117 | > |  | 
| 118 | > | ppStream.str(buf); | 
| 119 | > | delete [] buf; | 
| 120 | > |  | 
| 121 | > | } | 
| 122 | > | #endif | 
| 123 | > | // Create a scanner that reads from the input stream | 
| 124 | > | MDLexer lexer(ppStream); | 
| 125 | > | lexer.setFilename(filename); | 
| 126 | > | lexer.initDeferredLineCount(); | 
| 127 |  |  | 
| 128 | < | #ifdef IS_MPI | 
| 128 | > | // Create a parser that reads from the scanner | 
| 129 | > | MDParser parser(lexer); | 
| 130 | > | parser.setFilename(filename); | 
| 131 | > |  | 
| 132 | > | // Create an observer that synchorizes file name change | 
| 133 | > | FilenameObserver observer; | 
| 134 | > | observer.setLexer(&lexer); | 
| 135 | > | observer.setParser(&parser); | 
| 136 | > | lexer.setObserver(&observer); | 
| 137 |  |  | 
| 138 | < | if (worldRank == 0) { | 
| 139 | < | #endif // is_mpi | 
| 138 | > | antlr::ASTFactory factory; | 
| 139 | > | parser.initializeASTFactory(factory); | 
| 140 | > | parser.setASTFactory(&factory); | 
| 141 | > | parser.mdfile(); | 
| 142 | > |  | 
| 143 | > | // Create a tree parser that reads information into Globals | 
| 144 | > | MDTreeParser treeParser; | 
| 145 | > | treeParser.initializeASTFactory(factory); | 
| 146 | > | treeParser.setASTFactory(&factory); | 
| 147 | > | simParams = treeParser.walkTree(parser.getAST()); | 
| 148 | > | } | 
| 149 | > |  | 
| 150 |  |  | 
| 151 | < | set_interface_stamps(stamps, simParams); | 
| 151 | > | catch(antlr::MismatchedCharException& e) { | 
| 152 | > | sprintf(painCave.errMsg, | 
| 153 | > | "parser exception: %s %s:%d:%d\n", | 
| 154 | > | e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 155 | > | painCave.isFatal = 1; | 
| 156 | > | simError(); | 
| 157 | > | } | 
| 158 | > | catch(antlr::MismatchedTokenException &e) { | 
| 159 | > | sprintf(painCave.errMsg, | 
| 160 | > | "parser exception: %s %s:%d:%d\n", | 
| 161 | > | e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 162 | > | painCave.isFatal = 1; | 
| 163 | > | simError(); | 
| 164 | > | } | 
| 165 | > | catch(antlr::NoViableAltForCharException &e) { | 
| 166 | > | sprintf(painCave.errMsg, | 
| 167 | > | "parser exception: %s %s:%d:%d\n", | 
| 168 | > | e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 169 | > | painCave.isFatal = 1; | 
| 170 | > | simError(); | 
| 171 | > | } | 
| 172 | > | catch(antlr::NoViableAltException &e) { | 
| 173 | > | sprintf(painCave.errMsg, | 
| 174 | > | "parser exception: %s %s:%d:%d\n", | 
| 175 | > | e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 176 | > | painCave.isFatal = 1; | 
| 177 | > | simError(); | 
| 178 | > | } | 
| 179 |  |  | 
| 180 | < | #ifdef IS_MPI | 
| 181 | < |  | 
| 182 | < | mpiEventInit(); | 
| 183 | < |  | 
| 184 | < | #endif | 
| 185 | < |  | 
| 82 | < | yacc_BASS(mdFileName.c_str()); | 
| 83 | < |  | 
| 84 | < | #ifdef IS_MPI | 
| 85 | < |  | 
| 86 | < | throwMPIEvent(NULL); | 
| 87 | < | } else { | 
| 88 | < | set_interface_stamps(stamps, simParams); | 
| 89 | < | mpiEventInit(); | 
| 90 | < | MPIcheckPoint(); | 
| 91 | < | mpiEventLoop(); | 
| 180 | > | catch(antlr::TokenStreamRecognitionException& e) { | 
| 181 | > | sprintf(painCave.errMsg, | 
| 182 | > | "parser exception: %s %s:%d:%d\n", | 
| 183 | > | e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 184 | > | painCave.isFatal = 1; | 
| 185 | > | simError(); | 
| 186 |  | } | 
| 187 | < |  | 
| 188 | < | #endif | 
| 189 | < |  | 
| 187 | > |  | 
| 188 | > | catch(antlr::TokenStreamIOException& e) { | 
| 189 | > | sprintf(painCave.errMsg, | 
| 190 | > | "parser exception: %s\n", | 
| 191 | > | e.getMessage().c_str()); | 
| 192 | > | painCave.isFatal = 1; | 
| 193 | > | simError(); | 
| 194 | > | } | 
| 195 | > |  | 
| 196 | > | catch(antlr::TokenStreamException& e) { | 
| 197 | > | sprintf(painCave.errMsg, | 
| 198 | > | "parser exception: %s\n", | 
| 199 | > | e.getMessage().c_str()); | 
| 200 | > | painCave.isFatal = 1; | 
| 201 | > | simError(); | 
| 202 | > | } | 
| 203 | > | catch (antlr::RecognitionException& e) { | 
| 204 | > | sprintf(painCave.errMsg, | 
| 205 | > | "parser exception: %s %s:%d:%d\n", | 
| 206 | > | e.getMessage().c_str(),e.getFilename().c_str(), e.getLine(), e.getColumn()); | 
| 207 | > | painCave.isFatal = 1; | 
| 208 | > | simError(); | 
| 209 | > | } | 
| 210 | > | catch (antlr::CharStreamException& e) { | 
| 211 | > | sprintf(painCave.errMsg, | 
| 212 | > | "parser exception: %s\n", | 
| 213 | > | e.getMessage().c_str()); | 
| 214 | > | painCave.isFatal = 1; | 
| 215 | > | simError(); | 
| 216 | > | } | 
| 217 | > | catch (OpenMDException& e) { | 
| 218 | > | sprintf(painCave.errMsg, | 
| 219 | > | "%s\n", | 
| 220 | > | e.getMessage().c_str()); | 
| 221 | > | painCave.isFatal = 1; | 
| 222 | > | simError(); | 
| 223 | > | } | 
| 224 | > | catch (std::exception& e) { | 
| 225 | > | sprintf(painCave.errMsg, | 
| 226 | > | "parser exception: %s\n", | 
| 227 | > | e.what()); | 
| 228 | > | painCave.isFatal = 1; | 
| 229 | > | simError(); | 
| 230 | > | } | 
| 231 | > |  | 
| 232 | > | return simParams; | 
| 233 |  | } | 
| 234 |  |  | 
| 235 |  | SimInfo*  SimCreator::createSim(const std::string & mdFileName, | 
| 236 |  | bool loadInitCoords) { | 
| 237 | < |  | 
| 238 | < | MakeStamps * stamps = new MakeStamps(); | 
| 239 | < |  | 
| 240 | < | Globals * simParams = new Globals(); | 
| 241 | < |  | 
| 237 | > |  | 
| 238 | > | const int bufferSize = 65535; | 
| 239 | > | char buffer[bufferSize]; | 
| 240 | > | int lineNo = 0; | 
| 241 | > | std::string mdRawData; | 
| 242 | > | int metaDataBlockStart = -1; | 
| 243 | > | int metaDataBlockEnd = -1; | 
| 244 | > | int i; | 
| 245 | > | int mdOffset; | 
| 246 | > |  | 
| 247 | > | #ifdef IS_MPI | 
| 248 | > | const int masterNode = 0; | 
| 249 | > | if (worldRank == masterNode) { | 
| 250 | > | #endif | 
| 251 | > |  | 
| 252 | > | std::ifstream mdFile_(mdFileName.c_str()); | 
| 253 | > |  | 
| 254 | > | if (mdFile_.fail()) { | 
| 255 | > | sprintf(painCave.errMsg, | 
| 256 | > | "SimCreator: Cannot open file: %s\n", | 
| 257 | > | mdFileName.c_str()); | 
| 258 | > | painCave.isFatal = 1; | 
| 259 | > | simError(); | 
| 260 | > | } | 
| 261 | > |  | 
| 262 | > | mdFile_.getline(buffer, bufferSize); | 
| 263 | > | ++lineNo; | 
| 264 | > | std::string line = trimLeftCopy(buffer); | 
| 265 | > | i = CaseInsensitiveFind(line, "<OpenMD"); | 
| 266 | > | if (static_cast<size_t>(i) == string::npos) { | 
| 267 | > | // try the older file strings to see if that works: | 
| 268 | > | i = CaseInsensitiveFind(line, "<OOPSE"); | 
| 269 | > | } | 
| 270 | > |  | 
| 271 | > | if (static_cast<size_t>(i) == string::npos) { | 
| 272 | > | // still no luck! | 
| 273 | > | sprintf(painCave.errMsg, | 
| 274 | > | "SimCreator: File: %s is not a valid OpenMD file!\n", | 
| 275 | > | mdFileName.c_str()); | 
| 276 | > | painCave.isFatal = 1; | 
| 277 | > | simError(); | 
| 278 | > | } | 
| 279 | > |  | 
| 280 | > | //scan through the input stream and find MetaData tag | 
| 281 | > | while(mdFile_.getline(buffer, bufferSize)) { | 
| 282 | > | ++lineNo; | 
| 283 | > |  | 
| 284 | > | std::string line = trimLeftCopy(buffer); | 
| 285 | > | if (metaDataBlockStart == -1) { | 
| 286 | > | i = CaseInsensitiveFind(line, "<MetaData>"); | 
| 287 | > | if (i != string::npos) { | 
| 288 | > | metaDataBlockStart = lineNo; | 
| 289 | > | mdOffset = mdFile_.tellg(); | 
| 290 | > | } | 
| 291 | > | } else { | 
| 292 | > | i = CaseInsensitiveFind(line, "</MetaData>"); | 
| 293 | > | if (i != string::npos) { | 
| 294 | > | metaDataBlockEnd = lineNo; | 
| 295 | > | } | 
| 296 | > | } | 
| 297 | > | } | 
| 298 | > |  | 
| 299 | > | if (metaDataBlockStart == -1) { | 
| 300 | > | sprintf(painCave.errMsg, | 
| 301 | > | "SimCreator: File: %s did not contain a <MetaData> tag!\n", | 
| 302 | > | mdFileName.c_str()); | 
| 303 | > | painCave.isFatal = 1; | 
| 304 | > | simError(); | 
| 305 | > | } | 
| 306 | > | if (metaDataBlockEnd == -1) { | 
| 307 | > | sprintf(painCave.errMsg, | 
| 308 | > | "SimCreator: File: %s did not contain a closed MetaData block!\n", | 
| 309 | > | mdFileName.c_str()); | 
| 310 | > | painCave.isFatal = 1; | 
| 311 | > | simError(); | 
| 312 | > | } | 
| 313 | > |  | 
| 314 | > | mdFile_.clear(); | 
| 315 | > | mdFile_.seekg(0); | 
| 316 | > | mdFile_.seekg(mdOffset); | 
| 317 | > |  | 
| 318 | > | mdRawData.clear(); | 
| 319 | > |  | 
| 320 | > | for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) { | 
| 321 | > | mdFile_.getline(buffer, bufferSize); | 
| 322 | > | mdRawData += buffer; | 
| 323 | > | mdRawData += "\n"; | 
| 324 | > | } | 
| 325 | > |  | 
| 326 | > | mdFile_.close(); | 
| 327 | > |  | 
| 328 | > | #ifdef IS_MPI | 
| 329 | > | } | 
| 330 | > | #endif | 
| 331 | > |  | 
| 332 | > | std::stringstream rawMetaDataStream(mdRawData); | 
| 333 | > |  | 
| 334 |  | //parse meta-data file | 
| 335 | < | parseFile(mdFileName, stamps, simParams); | 
| 335 | > | Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1); | 
| 336 |  |  | 
| 337 |  | //create the force field | 
| 338 | < | ForceField * ff = ForceFieldFactory::getInstance() | 
| 339 | < | ->createForceField(simParams->getForceField()); | 
| 111 | < |  | 
| 338 | > | ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField()); | 
| 339 | > |  | 
| 340 |  | if (ff == NULL) { | 
| 341 |  | sprintf(painCave.errMsg, | 
| 342 |  | "ForceField Factory can not create %s force field\n", | 
| 369 |  | } | 
| 370 |  |  | 
| 371 |  | ff->parse(forcefieldFileName); | 
| 372 | < |  | 
| 145 | < | //extract the molecule stamps | 
| 146 | < | std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs; | 
| 147 | < | compList(stamps, simParams, moleculeStampPairs); | 
| 148 | < |  | 
| 372 | > | ff->setFortranForceOptions(); | 
| 373 |  | //create SimInfo | 
| 374 | < | SimInfo * info = new SimInfo(stamps, moleculeStampPairs, ff, simParams); | 
| 374 | > | SimInfo * info = new SimInfo(ff, simParams); | 
| 375 | > |  | 
| 376 | > | info->setRawMetaData(mdRawData); | 
| 377 |  |  | 
| 378 | < | //gather parameters (SimCreator only retrieves part of the parameters) | 
| 378 | > | //gather parameters (SimCreator only retrieves part of the | 
| 379 | > | //parameters) | 
| 380 |  | gatherParameters(info, mdFileName); | 
| 381 |  |  | 
| 382 |  | //divide the molecules and determine the global index of molecules | 
| 388 |  | createMolecules(info); | 
| 389 |  |  | 
| 390 |  |  | 
| 391 | < | //allocate memory for DataStorage(circular reference, need to break it) | 
| 391 | > | //allocate memory for DataStorage(circular reference, need to | 
| 392 | > | //break it) | 
| 393 |  | info->setSnapshotManager(new SimSnapshotManager(info)); | 
| 394 |  |  | 
| 395 | < | //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the | 
| 396 | < | //global index will never change again). Local indices of atoms and rigidbodies are already set by | 
| 397 | < | //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager. | 
| 395 | > | //set the global index of atoms, rigidbodies and cutoffgroups | 
| 396 | > | //(only need to be set once, the global index will never change | 
| 397 | > | //again). Local indices of atoms and rigidbodies are already set | 
| 398 | > | //by MoleculeCreator class which actually delegates the | 
| 399 | > | //responsibility to LocalIndexManager. | 
| 400 |  | setGlobalIndex(info); | 
| 401 |  |  | 
| 402 | < | //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point | 
| 403 | < | //atoms don't have the global index yet  (their global index are all initialized to -1). | 
| 404 | < | //Therefore we have to call addExcludePairs explicitly here. A way to work around is that | 
| 405 | < | //we can determine the beginning global indices of atoms before they get created. | 
| 402 | > | //Although addInteractionPairs is called inside SimInfo's addMolecule | 
| 403 | > | //method, at that point atoms don't have the global index yet | 
| 404 | > | //(their global index are all initialized to -1).  Therefore we | 
| 405 | > | //have to call addInteractionPairs explicitly here. A way to work | 
| 406 | > | //around is that we can determine the beginning global indices of | 
| 407 | > | //atoms before they get created. | 
| 408 |  | SimInfo::MoleculeIterator mi; | 
| 409 |  | Molecule* mol; | 
| 410 |  | for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 411 | < | info->addExcludePairs(mol); | 
| 411 | > | info->addInteractionPairs(mol); | 
| 412 |  | } | 
| 413 |  |  | 
| 414 |  | if (loadInitCoords) | 
| 415 | < | loadCoordinates(info); | 
| 415 | > | loadCoordinates(info, mdFileName); | 
| 416 |  |  | 
| 417 |  | return info; | 
| 418 |  | } | 
| 448 |  |  | 
| 449 |  | #ifdef IS_MPI | 
| 450 |  | void SimCreator::divideMolecules(SimInfo *info) { | 
| 451 | < | double numerator; | 
| 452 | < | double denominator; | 
| 453 | < | double precast; | 
| 454 | < | double x; | 
| 455 | < | double y; | 
| 456 | < | double a; | 
| 451 | > | RealType numerator; | 
| 452 | > | RealType denominator; | 
| 453 | > | RealType precast; | 
| 454 | > | RealType x; | 
| 455 | > | RealType y; | 
| 456 | > | RealType a; | 
| 457 |  | int old_atoms; | 
| 458 |  | int add_atoms; | 
| 459 |  | int new_atoms; | 
| 477 |  | "\tthe number of molecules.  This will not result in a \n" | 
| 478 |  | "\tusable division of atoms for force decomposition.\n" | 
| 479 |  | "\tEither try a smaller number of processors, or run the\n" | 
| 480 | < | "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols); | 
| 480 | > | "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols); | 
| 481 |  |  | 
| 482 |  | painCave.isFatal = 1; | 
| 483 |  | simError(); | 
| 565 |  | //           Pacc(x) = exp(- a * x) | 
| 566 |  | // where a = penalty / (average atoms per molecule) | 
| 567 |  |  | 
| 568 | < | x = (double)(new_atoms - nTarget); | 
| 568 | > | x = (RealType)(new_atoms - nTarget); | 
| 569 |  | y = myRandom->rand(); | 
| 570 |  |  | 
| 571 |  | if (y < exp(- a * x)) { | 
| 595 |  | info->setMolToProcMap(molToProcMap); | 
| 596 |  | sprintf(checkPointMsg, | 
| 597 |  | "Successfully divided the molecules among the processors.\n"); | 
| 598 | < | MPIcheckPoint(); | 
| 598 | > | errorCheckPoint(); | 
| 599 |  | } | 
| 600 |  |  | 
| 601 |  | #endif | 
| 625 |  |  | 
| 626 |  | } //end for(int i=0) | 
| 627 |  | } | 
| 396 | – |  | 
| 397 | – | void SimCreator::compList(MakeStamps *stamps, Globals* simParams, | 
| 398 | – | std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) { | 
| 399 | – | int i; | 
| 400 | – | char * id; | 
| 401 | – | MoleculeStamp * currentStamp; | 
| 402 | – | Component** the_components = simParams->getComponents(); | 
| 403 | – | int n_components = simParams->getNComponents(); | 
| 628 |  |  | 
| 405 | – | if (!simParams->haveNMol()) { | 
| 406 | – | // we don't have the total number of molecules, so we assume it is | 
| 407 | – | // given in each component | 
| 408 | – |  | 
| 409 | – | for(i = 0; i < n_components; i++) { | 
| 410 | – | if (!the_components[i]->haveNMol()) { | 
| 411 | – | // we have a problem | 
| 412 | – | sprintf(painCave.errMsg, | 
| 413 | – | "SimCreator Error. No global NMol or component NMol given.\n" | 
| 414 | – | "\tCannot calculate the number of atoms.\n"); | 
| 415 | – |  | 
| 416 | – | painCave.isFatal = 1; | 
| 417 | – | simError(); | 
| 418 | – | } | 
| 419 | – |  | 
| 420 | – | id = the_components[i]->getType(); | 
| 421 | – |  | 
| 422 | – | currentStamp = stamps->getMolStamp(id); | 
| 423 | – | if (currentStamp == NULL) { | 
| 424 | – | sprintf(painCave.errMsg, | 
| 425 | – | "SimCreator error: Component \"%s\" was not found in the " | 
| 426 | – | "list of declared molecules\n", id); | 
| 427 | – |  | 
| 428 | – | painCave.isFatal = 1; | 
| 429 | – | simError(); | 
| 430 | – | } | 
| 431 | – |  | 
| 432 | – | moleculeStampPairs.push_back( | 
| 433 | – | std::make_pair(currentStamp, the_components[i]->getNMol())); | 
| 434 | – | } //end for (i = 0; i < n_components; i++) | 
| 435 | – | } else { | 
| 436 | – | sprintf(painCave.errMsg, "SimSetup error.\n" | 
| 437 | – | "\tSorry, the ability to specify total" | 
| 438 | – | " nMols and then give molfractions in the components\n" | 
| 439 | – | "\tis not currently supported." | 
| 440 | – | " Please give nMol in the components.\n"); | 
| 441 | – |  | 
| 442 | – | painCave.isFatal = 1; | 
| 443 | – | simError(); | 
| 444 | – | } | 
| 445 | – |  | 
| 446 | – | #ifdef IS_MPI | 
| 447 | – |  | 
| 448 | – | strcpy(checkPointMsg, "Component stamps successfully extracted\n"); | 
| 449 | – | MPIcheckPoint(); | 
| 450 | – |  | 
| 451 | – | #endif // is_mpi | 
| 452 | – |  | 
| 453 | – | } | 
| 454 | – |  | 
| 629 |  | void SimCreator::setGlobalIndex(SimInfo *info) { | 
| 630 |  | SimInfo::MoleculeIterator mi; | 
| 631 |  | Molecule::AtomIterator ai; | 
| 632 |  | Molecule::RigidBodyIterator ri; | 
| 633 |  | Molecule::CutoffGroupIterator ci; | 
| 634 | + | Molecule::IntegrableObjectIterator  ioi; | 
| 635 |  | Molecule * mol; | 
| 636 |  | Atom * atom; | 
| 637 |  | RigidBody * rb; | 
| 640 |  | int beginRigidBodyIndex; | 
| 641 |  | int beginCutoffGroupIndex; | 
| 642 |  | int nGlobalAtoms = info->getNGlobalAtoms(); | 
| 643 | < |  | 
| 643 | > |  | 
| 644 | > | /**@todo fixme */ | 
| 645 |  | #ifndef IS_MPI | 
| 646 |  |  | 
| 647 |  | beginAtomIndex = 0; | 
| 728 |  | // to get the full globalGroupMembership array (We think). | 
| 729 |  | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | 
| 730 |  | // docs said we could. | 
| 731 | < | std::vector<int> tmpGroupMembership(nGlobalAtoms, 0); | 
| 731 | > | std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); | 
| 732 |  | MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, | 
| 733 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 734 |  | info->setGlobalGroupMembership(tmpGroupMembership); | 
| 740 |  | std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); | 
| 741 |  |  | 
| 742 |  | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 567 | – |  | 
| 743 |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 744 |  | globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); | 
| 745 |  | } | 
| 746 |  | } | 
| 747 |  |  | 
| 748 |  | #ifdef IS_MPI | 
| 749 | < | std::vector<int> tmpMolMembership(nGlobalAtoms, 0); | 
| 749 | > | std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0); | 
| 750 |  |  | 
| 751 |  | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, | 
| 752 |  | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 755 |  | #else | 
| 756 |  | info->setGlobalMolMembership(globalMolMembership); | 
| 757 |  | #endif | 
| 758 | + |  | 
| 759 | + | // nIOPerMol holds the number of integrable objects per molecule | 
| 760 | + | // here the molecules are listed by their global indices. | 
| 761 | + |  | 
| 762 | + | std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0); | 
| 763 | + | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 764 | + | nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects(); | 
| 765 | + | } | 
| 766 |  |  | 
| 767 | + | #ifdef IS_MPI | 
| 768 | + | std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0); | 
| 769 | + | MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], | 
| 770 | + | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 771 | + | #else | 
| 772 | + | std::vector<int> numIntegrableObjectsPerMol = nIOPerMol; | 
| 773 | + | #endif | 
| 774 | + |  | 
| 775 | + | std::vector<int> startingIOIndexForMol(info->getNGlobalMolecules()); | 
| 776 | + |  | 
| 777 | + | int startingIndex = 0; | 
| 778 | + | for (int i = 0; i < info->getNGlobalMolecules(); i++) { | 
| 779 | + | startingIOIndexForMol[i] = startingIndex; | 
| 780 | + | startingIndex += numIntegrableObjectsPerMol[i]; | 
| 781 | + | } | 
| 782 | + |  | 
| 783 | + | std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | 
| 784 | + | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 785 | + | int myGlobalIndex = mol->getGlobalIndex(); | 
| 786 | + | int globalIO = startingIOIndexForMol[myGlobalIndex]; | 
| 787 | + | for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL; | 
| 788 | + | integrableObject = mol->nextIntegrableObject(ioi)) { | 
| 789 | + | integrableObject->setGlobalIntegrableObjectIndex(globalIO); | 
| 790 | + | IOIndexToIntegrableObject[globalIO] = integrableObject; | 
| 791 | + | globalIO++; | 
| 792 | + | } | 
| 793 | + | } | 
| 794 | + |  | 
| 795 | + | info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); | 
| 796 | + |  | 
| 797 |  | } | 
| 798 |  |  | 
| 799 | < | void SimCreator::loadCoordinates(SimInfo* info) { | 
| 799 | > | void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { | 
| 800 |  | Globals* simParams; | 
| 801 |  | simParams = info->getSimParams(); | 
| 802 |  |  | 
| 590 | – | if (!simParams->haveInitialConfig()) { | 
| 591 | – | sprintf(painCave.errMsg, | 
| 592 | – | "Cannot intialize a simulation without an initial configuration file.\n"); | 
| 593 | – | painCave.isFatal = 1;; | 
| 594 | – | simError(); | 
| 595 | – | } | 
| 803 |  |  | 
| 804 | < | DumpReader reader(info, simParams->getInitialConfig()); | 
| 804 | > | DumpReader reader(info, mdFileName); | 
| 805 |  | int nframes = reader.getNFrames(); | 
| 806 |  |  | 
| 807 |  | if (nframes > 0) { | 
| 810 |  | //invalid initial coordinate file | 
| 811 |  | sprintf(painCave.errMsg, | 
| 812 |  | "Initial configuration file %s should at least contain one frame\n", | 
| 813 | < | simParams->getInitialConfig().c_str()); | 
| 813 | > | mdFileName.c_str()); | 
| 814 |  | painCave.isFatal = 1; | 
| 815 |  | simError(); | 
| 816 |  | } | 
| 819 |  | info->getSnapshotManager()->advance(); | 
| 820 |  | } | 
| 821 |  |  | 
| 822 | < | } //end namespace oopse | 
| 822 | > | } //end namespace OpenMD | 
| 823 |  |  | 
| 824 |  |  |