# | Line 35 | Line 35 | |
---|---|---|
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). |
38 | > | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 | > | * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 | > | * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 | */ | |
42 | ||
43 | /** | |
44 | * @file SimCreator.cpp | |
45 | * @author tlin | |
46 | * @date 11/03/2004 | |
46 | – | * @time 13:51am |
47 | * @version 1.0 | |
48 | */ | |
49 | + | |
50 | + | #ifdef IS_MPI |
51 | + | #include "mpi.h" |
52 | + | #include "math/ParallelRandNumGen.hpp" |
53 | + | #endif |
54 | + | |
55 | #include <exception> | |
56 | #include <iostream> | |
57 | #include <sstream> | |
# | Line 55 | Line 61 | |
61 | #include "brains/SimCreator.hpp" | |
62 | #include "brains/SimSnapshotManager.hpp" | |
63 | #include "io/DumpReader.hpp" | |
64 | < | #include "UseTheForce/ForceFieldFactory.hpp" |
64 | > | #include "brains/ForceField.hpp" |
65 | #include "utils/simError.h" | |
66 | #include "utils/StringUtils.hpp" | |
67 | + | #include "utils/Revision.hpp" |
68 | #include "math/SeqRandNumGen.hpp" | |
69 | #include "mdParser/MDLexer.hpp" | |
70 | #include "mdParser/MDParser.hpp" | |
# | Line 75 | Line 82 | |
82 | #include "antlr/NoViableAltForCharException.hpp" | |
83 | #include "antlr/NoViableAltException.hpp" | |
84 | ||
85 | < | #ifdef IS_MPI |
86 | < | #include "math/ParallelRandNumGen.hpp" |
87 | < | #endif |
85 | > | #include "types/DirectionalAdapter.hpp" |
86 | > | #include "types/MultipoleAdapter.hpp" |
87 | > | #include "types/EAMAdapter.hpp" |
88 | > | #include "types/SuttonChenAdapter.hpp" |
89 | > | #include "types/PolarizableAdapter.hpp" |
90 | > | #include "types/FixedChargeAdapter.hpp" |
91 | > | #include "types/FluctuatingChargeAdapter.hpp" |
92 | ||
93 | + | |
94 | namespace OpenMD { | |
95 | ||
96 | < | Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){ |
96 | > | Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){ |
97 | Globals* simParams = NULL; | |
98 | try { | |
99 | ||
# | Line 90 | Line 102 | namespace OpenMD { | |
102 | #ifdef IS_MPI | |
103 | int streamSize; | |
104 | const int masterNode = 0; | |
105 | < | int commStatus; |
105 | > | |
106 | if (worldRank == masterNode) { | |
107 | < | #endif |
108 | < | |
107 | > | MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
108 | > | // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); |
109 | > | #endif |
110 | SimplePreprocessor preprocessor; | |
111 | < | preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); |
111 | > | preprocessor.preprocess(rawMetaDataStream, filename, |
112 | > | startOfMetaDataBlock, ppStream); |
113 | ||
114 | #ifdef IS_MPI | |
115 | < | //brocasting the stream size |
115 | > | //broadcasting the stream size |
116 | streamSize = ppStream.str().size() +1; | |
117 | < | commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); |
117 | > | MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
118 | > | MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), |
119 | > | streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
120 | ||
121 | < | commStatus = MPI_Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
122 | < | |
123 | < | |
121 | > | // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); |
122 | > | // MPI::COMM_WORLD.Bcast(static_cast<void*>(const_cast<char*>(ppStream.str().c_str())), |
123 | > | // streamSize, MPI::CHAR, masterNode); |
124 | > | |
125 | } else { | |
109 | – | //get stream size |
110 | – | commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); |
126 | ||
127 | + | MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
128 | + | // MPI::COMM_WORLD.Bcast(&mdFileVersion, 1, MPI::INT, masterNode); |
129 | + | |
130 | + | //get stream size |
131 | + | MPI_Bcast(&streamSize, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
132 | + | // MPI::COMM_WORLD.Bcast(&streamSize, 1, MPI::LONG, masterNode); |
133 | char* buf = new char[streamSize]; | |
134 | assert(buf); | |
135 | ||
136 | //receive file content | |
137 | < | commStatus = MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
138 | < | |
137 | > | MPI_Bcast(buf, streamSize, MPI_CHAR, masterNode, MPI_COMM_WORLD); |
138 | > | // MPI::COMM_WORLD.Bcast(buf, streamSize, MPI::CHAR, masterNode); |
139 | > | |
140 | ppStream.str(buf); | |
141 | delete [] buf; | |
120 | – | |
142 | } | |
143 | #endif | |
144 | // Create a scanner that reads from the input stream | |
# | Line 139 | Line 160 | namespace OpenMD { | |
160 | parser.initializeASTFactory(factory); | |
161 | parser.setASTFactory(&factory); | |
162 | parser.mdfile(); | |
142 | – | |
163 | // Create a tree parser that reads information into Globals | |
164 | MDTreeParser treeParser; | |
165 | treeParser.initializeASTFactory(factory); | |
# | Line 229 | Line 249 | namespace OpenMD { | |
249 | simError(); | |
250 | } | |
251 | ||
252 | + | simParams->setMDfileVersion(mdFileVersion); |
253 | return simParams; | |
254 | } | |
255 | ||
256 | SimInfo* SimCreator::createSim(const std::string & mdFileName, | |
257 | bool loadInitCoords) { | |
258 | < | |
258 | > | |
259 | const int bufferSize = 65535; | |
260 | char buffer[bufferSize]; | |
261 | int lineNo = 0; | |
262 | std::string mdRawData; | |
263 | int metaDataBlockStart = -1; | |
264 | int metaDataBlockEnd = -1; | |
265 | < | int i; |
266 | < | int mdOffset; |
265 | > | int i, j; |
266 | > | streamoff mdOffset; |
267 | > | int mdFileVersion; |
268 | ||
269 | + | // Create a string for embedding the version information in the MetaData |
270 | + | std::string version; |
271 | + | version.assign("## Last run using OpenMD Version: "); |
272 | + | version.append(OPENMD_VERSION_MAJOR); |
273 | + | version.append("."); |
274 | + | version.append(OPENMD_VERSION_MINOR); |
275 | + | |
276 | + | std::string svnrev(g_REVISION, strnlen(g_REVISION, 20)); |
277 | + | //convert a macro from compiler to a string in c++ |
278 | + | // STR_DEFINE(svnrev, SVN_REV ); |
279 | + | version.append(" Revision: "); |
280 | + | // If there's no SVN revision, just call this the RELEASE revision. |
281 | + | if (!svnrev.empty()) { |
282 | + | version.append(svnrev); |
283 | + | } else { |
284 | + | version.append("RELEASE"); |
285 | + | } |
286 | + | |
287 | #ifdef IS_MPI | |
288 | const int masterNode = 0; | |
289 | if (worldRank == masterNode) { | |
290 | #endif | |
291 | ||
292 | < | std::ifstream mdFile_(mdFileName.c_str()); |
292 | > | std::ifstream mdFile_; |
293 | > | mdFile_.open(mdFileName.c_str(), ifstream::in | ifstream::binary); |
294 | ||
295 | if (mdFile_.fail()) { | |
296 | sprintf(painCave.errMsg, | |
# | Line 276 | Line 317 | namespace OpenMD { | |
317 | painCave.isFatal = 1; | |
318 | simError(); | |
319 | } | |
320 | + | |
321 | + | // found the correct opening string, now try to get the file |
322 | + | // format version number. |
323 | ||
324 | + | StringTokenizer tokenizer(line, "=<> \t\n\r"); |
325 | + | std::string fileType = tokenizer.nextToken(); |
326 | + | toUpper(fileType); |
327 | + | |
328 | + | mdFileVersion = 0; |
329 | + | |
330 | + | if (fileType == "OPENMD") { |
331 | + | while (tokenizer.hasMoreTokens()) { |
332 | + | std::string token(tokenizer.nextToken()); |
333 | + | toUpper(token); |
334 | + | if (token == "VERSION") { |
335 | + | mdFileVersion = tokenizer.nextTokenAsInt(); |
336 | + | break; |
337 | + | } |
338 | + | } |
339 | + | } |
340 | + | |
341 | //scan through the input stream and find MetaData tag | |
342 | while(mdFile_.getline(buffer, bufferSize)) { | |
343 | ++lineNo; | |
# | Line 317 | Line 378 | namespace OpenMD { | |
378 | ||
379 | mdRawData.clear(); | |
380 | ||
381 | + | bool foundVersion = false; |
382 | + | |
383 | for (int i = 0; i < metaDataBlockEnd - metaDataBlockStart - 1; ++i) { | |
384 | mdFile_.getline(buffer, bufferSize); | |
385 | < | mdRawData += buffer; |
385 | > | std::string line = trimLeftCopy(buffer); |
386 | > | j = CaseInsensitiveFind(line, "## Last run using OpenMD Version"); |
387 | > | if (static_cast<size_t>(j) != string::npos) { |
388 | > | foundVersion = true; |
389 | > | mdRawData += version; |
390 | > | } else { |
391 | > | mdRawData += buffer; |
392 | > | } |
393 | mdRawData += "\n"; | |
394 | } | |
395 | < | |
395 | > | |
396 | > | if (!foundVersion) mdRawData += version + "\n"; |
397 | > | |
398 | mdFile_.close(); | |
399 | ||
400 | #ifdef IS_MPI | |
# | Line 332 | Line 404 | namespace OpenMD { | |
404 | std::stringstream rawMetaDataStream(mdRawData); | |
405 | ||
406 | //parse meta-data file | |
407 | < | Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1); |
407 | > | Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion, |
408 | > | metaDataBlockStart + 1); |
409 | ||
410 | //create the force field | |
411 | < | ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField()); |
411 | > | ForceField * ff = new ForceField(simParams->getForceField()); |
412 | ||
413 | if (ff == NULL) { | |
414 | sprintf(painCave.errMsg, | |
# | Line 369 | Line 442 | namespace OpenMD { | |
442 | } | |
443 | ||
444 | ff->parse(forcefieldFileName); | |
372 | – | ff->setFortranForceOptions(); |
445 | //create SimInfo | |
446 | SimInfo * info = new SimInfo(ff, simParams); | |
447 | ||
# | Line 387 | Line 459 | namespace OpenMD { | |
459 | //create the molecules | |
460 | createMolecules(info); | |
461 | ||
462 | < | |
462 | > | //find the storage layout |
463 | > | |
464 | > | int storageLayout = computeStorageLayout(info); |
465 | > | |
466 | //allocate memory for DataStorage(circular reference, need to | |
467 | //break it) | |
468 | < | info->setSnapshotManager(new SimSnapshotManager(info)); |
468 | > | info->setSnapshotManager(new SimSnapshotManager(info, storageLayout)); |
469 | ||
470 | //set the global index of atoms, rigidbodies and cutoffgroups | |
471 | //(only need to be set once, the global index will never change | |
# | Line 413 | Line 488 | namespace OpenMD { | |
488 | ||
489 | if (loadInitCoords) | |
490 | loadCoordinates(info, mdFileName); | |
416 | – | |
491 | return info; | |
492 | } | |
493 | ||
# | Line 448 | Line 522 | namespace OpenMD { | |
522 | ||
523 | #ifdef IS_MPI | |
524 | void SimCreator::divideMolecules(SimInfo *info) { | |
451 | – | RealType numerator; |
452 | – | RealType denominator; |
453 | – | RealType precast; |
454 | – | RealType x; |
455 | – | RealType y; |
525 | RealType a; | |
457 | – | int old_atoms; |
458 | – | int add_atoms; |
459 | – | int new_atoms; |
460 | – | int nTarget; |
461 | – | int done; |
462 | – | int i; |
463 | – | int j; |
464 | – | int loops; |
465 | – | int which_proc; |
526 | int nProcessors; | |
527 | std::vector<int> atomsPerProc; | |
528 | int nGlobalMols = info->getNGlobalMolecules(); | |
529 | < | std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition: |
529 | > | std::vector<int> molToProcMap(nGlobalMols, -1); // default to an |
530 | > | // error |
531 | > | // condition: |
532 | ||
533 | < | MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); |
533 | > | MPI_Comm_size( MPI_COMM_WORLD, &nProcessors); |
534 | > | //nProcessors = MPI::COMM_WORLD.Get_size(); |
535 | ||
536 | if (nProcessors > nGlobalMols) { | |
537 | sprintf(painCave.errMsg, | |
# | Line 477 | Line 540 | namespace OpenMD { | |
540 | "\tthe number of molecules. This will not result in a \n" | |
541 | "\tusable division of atoms for force decomposition.\n" | |
542 | "\tEither try a smaller number of processors, or run the\n" | |
543 | < | "\tsingle-processor version of OpenMD.\n", nProcessors, nGlobalMols); |
543 | > | "\tsingle-processor version of OpenMD.\n", nProcessors, |
544 | > | nGlobalMols); |
545 | ||
546 | painCave.isFatal = 1; | |
547 | simError(); | |
548 | } | |
549 | ||
486 | – | int seedValue; |
550 | Globals * simParams = info->getSimParams(); | |
551 | < | SeqRandNumGen* myRandom; //divide labor does not need Parallel random number generator |
551 | > | SeqRandNumGen* myRandom; //divide labor does not need Parallel |
552 | > | //random number generator |
553 | if (simParams->haveSeed()) { | |
554 | < | seedValue = simParams->getSeed(); |
554 | > | int seedValue = simParams->getSeed(); |
555 | myRandom = new SeqRandNumGen(seedValue); | |
556 | }else { | |
557 | myRandom = new SeqRandNumGen(); | |
# | Line 500 | Line 564 | namespace OpenMD { | |
564 | atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0); | |
565 | ||
566 | if (worldRank == 0) { | |
567 | < | numerator = info->getNGlobalAtoms(); |
568 | < | denominator = nProcessors; |
569 | < | precast = numerator / denominator; |
570 | < | nTarget = (int)(precast + 0.5); |
567 | > | RealType numerator = info->getNGlobalAtoms(); |
568 | > | RealType denominator = nProcessors; |
569 | > | RealType precast = numerator / denominator; |
570 | > | int nTarget = (int)(precast + 0.5); |
571 | ||
572 | < | for(i = 0; i < nGlobalMols; i++) { |
573 | < | done = 0; |
574 | < | loops = 0; |
572 | > | for(int i = 0; i < nGlobalMols; i++) { |
573 | > | |
574 | > | int done = 0; |
575 | > | int loops = 0; |
576 | ||
577 | while (!done) { | |
578 | loops++; | |
579 | ||
580 | // Pick a processor at random | |
581 | ||
582 | < | which_proc = (int) (myRandom->rand() * nProcessors); |
582 | > | int which_proc = (int) (myRandom->rand() * nProcessors); |
583 | ||
584 | //get the molecule stamp first | |
585 | int stampId = info->getMoleculeStampId(i); | |
586 | MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId); | |
587 | ||
588 | // How many atoms does this processor have so far? | |
589 | < | old_atoms = atomsPerProc[which_proc]; |
590 | < | add_atoms = moleculeStamp->getNAtoms(); |
591 | < | new_atoms = old_atoms + add_atoms; |
589 | > | int old_atoms = atomsPerProc[which_proc]; |
590 | > | int add_atoms = moleculeStamp->getNAtoms(); |
591 | > | int new_atoms = old_atoms + add_atoms; |
592 | ||
593 | // If we've been through this loop too many times, we need | |
594 | // to just give up and assign the molecule to this processor | |
595 | // and be done with it. | |
596 | ||
597 | if (loops > 100) { | |
598 | + | |
599 | sprintf(painCave.errMsg, | |
600 | < | "I've tried 100 times to assign molecule %d to a " |
601 | < | " processor, but can't find a good spot.\n" |
602 | < | "I'm assigning it at random to processor %d.\n", |
600 | > | "There have been 100 attempts to assign molecule %d to an\n" |
601 | > | "\tunderworked processor, but there's no good place to\n" |
602 | > | "\tleave it. OpenMD is assigning it at random to processor %d.\n", |
603 | i, which_proc); | |
604 | < | |
604 | > | |
605 | painCave.isFatal = 0; | |
606 | + | painCave.severity = OPENMD_INFO; |
607 | simError(); | |
608 | ||
609 | molToProcMap[i] = which_proc; | |
# | Line 565 | Line 632 | namespace OpenMD { | |
632 | // Pacc(x) = exp(- a * x) | |
633 | // where a = penalty / (average atoms per molecule) | |
634 | ||
635 | < | x = (RealType)(new_atoms - nTarget); |
636 | < | y = myRandom->rand(); |
635 | > | RealType x = (RealType)(new_atoms - nTarget); |
636 | > | RealType y = myRandom->rand(); |
637 | ||
638 | if (y < exp(- a * x)) { | |
639 | molToProcMap[i] = which_proc; | |
# | Line 581 | Line 648 | namespace OpenMD { | |
648 | } | |
649 | ||
650 | delete myRandom; | |
651 | < | |
651 | > | |
652 | // Spray out this nonsense to all other processors: | |
586 | – | |
653 | MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | |
654 | + | // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); |
655 | } else { | |
656 | ||
657 | // Listen to your marching orders from processor 0: | |
591 | – | |
658 | MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD); | |
659 | + | // MPI::COMM_WORLD.Bcast(&molToProcMap[0], nGlobalMols, MPI::INT, 0); |
660 | + | |
661 | } | |
662 | ||
663 | info->setMolToProcMap(molToProcMap); | |
# | Line 612 | Line 680 | namespace OpenMD { | |
680 | #endif | |
681 | ||
682 | stampId = info->getMoleculeStampId(i); | |
683 | < | Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId), |
684 | < | stampId, i, info->getLocalIndexManager()); |
683 | > | Molecule * mol = molCreator.createMolecule(info->getForceField(), |
684 | > | info->getMoleculeStamp(stampId), |
685 | > | stampId, i, |
686 | > | info->getLocalIndexManager()); |
687 | ||
688 | info->addMolecule(mol); | |
689 | ||
# | Line 626 | Line 696 | namespace OpenMD { | |
696 | } //end for(int i=0) | |
697 | } | |
698 | ||
699 | + | int SimCreator::computeStorageLayout(SimInfo* info) { |
700 | + | |
701 | + | Globals* simParams = info->getSimParams(); |
702 | + | int nRigidBodies = info->getNGlobalRigidBodies(); |
703 | + | set<AtomType*> atomTypes = info->getSimulatedAtomTypes(); |
704 | + | set<AtomType*>::iterator i; |
705 | + | bool hasDirectionalAtoms = false; |
706 | + | bool hasFixedCharge = false; |
707 | + | bool hasDipoles = false; |
708 | + | bool hasQuadrupoles = false; |
709 | + | bool hasPolarizable = false; |
710 | + | bool hasFluctuatingCharge = false; |
711 | + | bool hasMetallic = false; |
712 | + | int storageLayout = 0; |
713 | + | storageLayout |= DataStorage::dslPosition; |
714 | + | storageLayout |= DataStorage::dslVelocity; |
715 | + | storageLayout |= DataStorage::dslForce; |
716 | + | |
717 | + | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
718 | + | |
719 | + | DirectionalAdapter da = DirectionalAdapter( (*i) ); |
720 | + | MultipoleAdapter ma = MultipoleAdapter( (*i) ); |
721 | + | EAMAdapter ea = EAMAdapter( (*i) ); |
722 | + | SuttonChenAdapter sca = SuttonChenAdapter( (*i) ); |
723 | + | PolarizableAdapter pa = PolarizableAdapter( (*i) ); |
724 | + | FixedChargeAdapter fca = FixedChargeAdapter( (*i) ); |
725 | + | FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter( (*i) ); |
726 | + | |
727 | + | if (da.isDirectional()){ |
728 | + | hasDirectionalAtoms = true; |
729 | + | } |
730 | + | if (ma.isDipole()){ |
731 | + | hasDipoles = true; |
732 | + | } |
733 | + | if (ma.isQuadrupole()){ |
734 | + | hasQuadrupoles = true; |
735 | + | } |
736 | + | if (ea.isEAM() || sca.isSuttonChen()){ |
737 | + | hasMetallic = true; |
738 | + | } |
739 | + | if ( fca.isFixedCharge() ){ |
740 | + | hasFixedCharge = true; |
741 | + | } |
742 | + | if ( fqa.isFluctuatingCharge() ){ |
743 | + | hasFluctuatingCharge = true; |
744 | + | } |
745 | + | if ( pa.isPolarizable() ){ |
746 | + | hasPolarizable = true; |
747 | + | } |
748 | + | } |
749 | + | |
750 | + | if (nRigidBodies > 0 || hasDirectionalAtoms) { |
751 | + | storageLayout |= DataStorage::dslAmat; |
752 | + | if(storageLayout & DataStorage::dslVelocity) { |
753 | + | storageLayout |= DataStorage::dslAngularMomentum; |
754 | + | } |
755 | + | if (storageLayout & DataStorage::dslForce) { |
756 | + | storageLayout |= DataStorage::dslTorque; |
757 | + | } |
758 | + | } |
759 | + | if (hasDipoles) { |
760 | + | storageLayout |= DataStorage::dslDipole; |
761 | + | } |
762 | + | if (hasQuadrupoles) { |
763 | + | storageLayout |= DataStorage::dslQuadrupole; |
764 | + | } |
765 | + | if (hasFixedCharge || hasFluctuatingCharge) { |
766 | + | storageLayout |= DataStorage::dslSkippedCharge; |
767 | + | } |
768 | + | if (hasMetallic) { |
769 | + | storageLayout |= DataStorage::dslDensity; |
770 | + | storageLayout |= DataStorage::dslFunctional; |
771 | + | storageLayout |= DataStorage::dslFunctionalDerivative; |
772 | + | } |
773 | + | if (hasPolarizable) { |
774 | + | storageLayout |= DataStorage::dslElectricField; |
775 | + | } |
776 | + | if (hasFluctuatingCharge){ |
777 | + | storageLayout |= DataStorage::dslFlucQPosition; |
778 | + | if(storageLayout & DataStorage::dslVelocity) { |
779 | + | storageLayout |= DataStorage::dslFlucQVelocity; |
780 | + | } |
781 | + | if (storageLayout & DataStorage::dslForce) { |
782 | + | storageLayout |= DataStorage::dslFlucQForce; |
783 | + | } |
784 | + | } |
785 | + | |
786 | + | // if the user has asked for them, make sure we've got the memory for the |
787 | + | // objects defined. |
788 | + | |
789 | + | if (simParams->getOutputParticlePotential()) { |
790 | + | storageLayout |= DataStorage::dslParticlePot; |
791 | + | } |
792 | + | |
793 | + | if (simParams->havePrintHeatFlux()) { |
794 | + | if (simParams->getPrintHeatFlux()) { |
795 | + | storageLayout |= DataStorage::dslParticlePot; |
796 | + | } |
797 | + | } |
798 | + | |
799 | + | if (simParams->getOutputElectricField() | simParams->haveElectricField()) { |
800 | + | storageLayout |= DataStorage::dslElectricField; |
801 | + | } |
802 | + | |
803 | + | if (simParams->getOutputFluctuatingCharges()) { |
804 | + | storageLayout |= DataStorage::dslFlucQPosition; |
805 | + | storageLayout |= DataStorage::dslFlucQVelocity; |
806 | + | storageLayout |= DataStorage::dslFlucQForce; |
807 | + | } |
808 | + | |
809 | + | info->setStorageLayout(storageLayout); |
810 | + | |
811 | + | return storageLayout; |
812 | + | } |
813 | + | |
814 | void SimCreator::setGlobalIndex(SimInfo *info) { | |
815 | SimInfo::MoleculeIterator mi; | |
816 | Molecule::AtomIterator ai; | |
817 | Molecule::RigidBodyIterator ri; | |
818 | Molecule::CutoffGroupIterator ci; | |
819 | + | Molecule::BondIterator boi; |
820 | + | Molecule::BendIterator bei; |
821 | + | Molecule::TorsionIterator ti; |
822 | + | Molecule::InversionIterator ii; |
823 | Molecule::IntegrableObjectIterator ioi; | |
824 | < | Molecule * mol; |
825 | < | Atom * atom; |
826 | < | RigidBody * rb; |
827 | < | CutoffGroup * cg; |
824 | > | Molecule* mol; |
825 | > | Atom* atom; |
826 | > | RigidBody* rb; |
827 | > | CutoffGroup* cg; |
828 | > | Bond* bond; |
829 | > | Bend* bend; |
830 | > | Torsion* torsion; |
831 | > | Inversion* inversion; |
832 | int beginAtomIndex; | |
833 | int beginRigidBodyIndex; | |
834 | int beginCutoffGroupIndex; | |
835 | + | int beginBondIndex; |
836 | + | int beginBendIndex; |
837 | + | int beginTorsionIndex; |
838 | + | int beginInversionIndex; |
839 | int nGlobalAtoms = info->getNGlobalAtoms(); | |
840 | < | |
644 | < | /**@todo fixme */ |
645 | < | #ifndef IS_MPI |
840 | > | int nGlobalRigidBodies = info->getNGlobalRigidBodies(); |
841 | ||
842 | beginAtomIndex = 0; | |
843 | < | beginRigidBodyIndex = 0; |
843 | > | // The rigid body indices begin immediately after the atom indices: |
844 | > | beginRigidBodyIndex = info->getNGlobalAtoms(); |
845 | beginCutoffGroupIndex = 0; | |
846 | < | |
847 | < | #else |
848 | < | |
849 | < | int nproc; |
850 | < | int myNode; |
851 | < | |
656 | < | myNode = worldRank; |
657 | < | MPI_Comm_size(MPI_COMM_WORLD, &nproc); |
658 | < | |
659 | < | std::vector < int > tmpAtomsInProc(nproc, 0); |
660 | < | std::vector < int > tmpRigidBodiesInProc(nproc, 0); |
661 | < | std::vector < int > tmpCutoffGroupsInProc(nproc, 0); |
662 | < | std::vector < int > NumAtomsInProc(nproc, 0); |
663 | < | std::vector < int > NumRigidBodiesInProc(nproc, 0); |
664 | < | std::vector < int > NumCutoffGroupsInProc(nproc, 0); |
665 | < | |
666 | < | tmpAtomsInProc[myNode] = info->getNAtoms(); |
667 | < | tmpRigidBodiesInProc[myNode] = info->getNRigidBodies(); |
668 | < | tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups(); |
669 | < | |
670 | < | //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups |
671 | < | MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT, |
672 | < | MPI_SUM, MPI_COMM_WORLD); |
673 | < | MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, |
674 | < | MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
675 | < | MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, |
676 | < | MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
677 | < | |
678 | < | beginAtomIndex = 0; |
679 | < | beginRigidBodyIndex = 0; |
680 | < | beginCutoffGroupIndex = 0; |
681 | < | |
682 | < | for(int i = 0; i < myNode; i++) { |
683 | < | beginAtomIndex += NumAtomsInProc[i]; |
684 | < | beginRigidBodyIndex += NumRigidBodiesInProc[i]; |
685 | < | beginCutoffGroupIndex += NumCutoffGroupsInProc[i]; |
686 | < | } |
687 | < | |
688 | < | #endif |
689 | < | |
690 | < | //rigidbody's index begins right after atom's |
691 | < | beginRigidBodyIndex += info->getNGlobalAtoms(); |
692 | < | |
693 | < | for(mol = info->beginMolecule(mi); mol != NULL; |
694 | < | mol = info->nextMolecule(mi)) { |
846 | > | beginBondIndex = 0; |
847 | > | beginBendIndex = 0; |
848 | > | beginTorsionIndex = 0; |
849 | > | beginInversionIndex = 0; |
850 | > | |
851 | > | for(int i = 0; i < info->getNGlobalMolecules(); i++) { |
852 | ||
853 | < | //local index(index in DataStorge) of atom is important |
854 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
855 | < | atom->setGlobalIndex(beginAtomIndex++); |
853 | > | #ifdef IS_MPI |
854 | > | if (info->getMolToProc(i) == worldRank) { |
855 | > | #endif |
856 | > | // stuff to do if I own this molecule |
857 | > | mol = info->getMoleculeByGlobalIndex(i); |
858 | > | |
859 | > | // The local index(index in DataStorge) of the atom is important: |
860 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
861 | > | atom = mol->nextAtom(ai)) { |
862 | > | atom->setGlobalIndex(beginAtomIndex++); |
863 | > | } |
864 | > | |
865 | > | for(rb = mol->beginRigidBody(ri); rb != NULL; |
866 | > | rb = mol->nextRigidBody(ri)) { |
867 | > | rb->setGlobalIndex(beginRigidBodyIndex++); |
868 | > | } |
869 | > | |
870 | > | // The local index of other objects only depends on the order |
871 | > | // of traversal: |
872 | > | for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
873 | > | cg = mol->nextCutoffGroup(ci)) { |
874 | > | cg->setGlobalIndex(beginCutoffGroupIndex++); |
875 | > | } |
876 | > | for(bond = mol->beginBond(boi); bond != NULL; |
877 | > | bond = mol->nextBond(boi)) { |
878 | > | bond->setGlobalIndex(beginBondIndex++); |
879 | > | } |
880 | > | for(bend = mol->beginBend(bei); bend != NULL; |
881 | > | bend = mol->nextBend(bei)) { |
882 | > | bend->setGlobalIndex(beginBendIndex++); |
883 | > | } |
884 | > | for(torsion = mol->beginTorsion(ti); torsion != NULL; |
885 | > | torsion = mol->nextTorsion(ti)) { |
886 | > | torsion->setGlobalIndex(beginTorsionIndex++); |
887 | > | } |
888 | > | for(inversion = mol->beginInversion(ii); inversion != NULL; |
889 | > | inversion = mol->nextInversion(ii)) { |
890 | > | inversion->setGlobalIndex(beginInversionIndex++); |
891 | > | } |
892 | > | |
893 | > | #ifdef IS_MPI |
894 | > | } else { |
895 | > | |
896 | > | // stuff to do if I don't own this molecule |
897 | > | |
898 | > | int stampId = info->getMoleculeStampId(i); |
899 | > | MoleculeStamp* stamp = info->getMoleculeStamp(stampId); |
900 | > | |
901 | > | beginAtomIndex += stamp->getNAtoms(); |
902 | > | beginRigidBodyIndex += stamp->getNRigidBodies(); |
903 | > | beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); |
904 | > | beginBondIndex += stamp->getNBonds(); |
905 | > | beginBendIndex += stamp->getNBends(); |
906 | > | beginTorsionIndex += stamp->getNTorsions(); |
907 | > | beginInversionIndex += stamp->getNInversions(); |
908 | } | |
909 | < | |
910 | < | for(rb = mol->beginRigidBody(ri); rb != NULL; |
911 | < | rb = mol->nextRigidBody(ri)) { |
912 | < | rb->setGlobalIndex(beginRigidBodyIndex++); |
704 | < | } |
705 | < | |
706 | < | //local index of cutoff group is trivial, it only depends on the order of travesing |
707 | < | for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
708 | < | cg = mol->nextCutoffGroup(ci)) { |
709 | < | cg->setGlobalIndex(beginCutoffGroupIndex++); |
710 | < | } |
711 | < | } |
712 | < | |
909 | > | #endif |
910 | > | |
911 | > | } //end for(int i=0) |
912 | > | |
913 | //fill globalGroupMembership | |
914 | std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); | |
915 | < | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
916 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
917 | < | |
915 | > | for(mol = info->beginMolecule(mi); mol != NULL; |
916 | > | mol = info->nextMolecule(mi)) { |
917 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
918 | > | cg = mol->nextCutoffGroup(ci)) { |
919 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
920 | globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex(); | |
921 | } | |
922 | ||
923 | } | |
924 | } | |
925 | < | |
925 | > | |
926 | #ifdef IS_MPI | |
927 | // Since the globalGroupMembership has been zero filled and we've only | |
928 | // poked values into the atoms we know, we can do an Allreduce | |
# | Line 729 | Line 930 | namespace OpenMD { | |
930 | // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | |
931 | // docs said we could. | |
932 | std::vector<int> tmpGroupMembership(info->getNGlobalAtoms(), 0); | |
933 | < | MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms, |
933 | > | MPI_Allreduce(&globalGroupMembership[0], |
934 | > | &tmpGroupMembership[0], nGlobalAtoms, |
935 | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | |
936 | + | // MPI::COMM_WORLD.Allreduce(&globalGroupMembership[0], |
937 | + | // &tmpGroupMembership[0], nGlobalAtoms, |
938 | + | // MPI::INT, MPI::SUM); |
939 | info->setGlobalGroupMembership(tmpGroupMembership); | |
940 | #else | |
941 | info->setGlobalGroupMembership(globalGroupMembership); | |
942 | #endif | |
943 | ||
944 | //fill molMembership | |
945 | < | std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0); |
945 | > | std::vector<int> globalMolMembership(info->getNGlobalAtoms() + |
946 | > | info->getNGlobalRigidBodies(), 0); |
947 | ||
948 | < | for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
948 | > | for(mol = info->beginMolecule(mi); mol != NULL; |
949 | > | mol = info->nextMolecule(mi)) { |
950 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
951 | globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex(); | |
952 | } | |
953 | + | for (rb = mol->beginRigidBody(ri); rb != NULL; |
954 | + | rb = mol->nextRigidBody(ri)) { |
955 | + | globalMolMembership[rb->getGlobalIndex()] = mol->getGlobalIndex(); |
956 | + | } |
957 | } | |
958 | ||
959 | #ifdef IS_MPI | |
960 | < | std::vector<int> tmpMolMembership(info->getNGlobalAtoms(), 0); |
961 | < | |
962 | < | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms, |
960 | > | std::vector<int> tmpMolMembership(info->getNGlobalAtoms() + |
961 | > | info->getNGlobalRigidBodies(), 0); |
962 | > | MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], |
963 | > | nGlobalAtoms + nGlobalRigidBodies, |
964 | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | |
965 | + | // MPI::COMM_WORLD.Allreduce(&globalMolMembership[0], &tmpMolMembership[0], |
966 | + | // nGlobalAtoms + nGlobalRigidBodies, |
967 | + | // MPI::INT, MPI::SUM); |
968 | ||
969 | info->setGlobalMolMembership(tmpMolMembership); | |
970 | #else | |
# | Line 760 | Line 975 | namespace OpenMD { | |
975 | // here the molecules are listed by their global indices. | |
976 | ||
977 | std::vector<int> nIOPerMol(info->getNGlobalMolecules(), 0); | |
978 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
978 | > | for (mol = info->beginMolecule(mi); mol != NULL; |
979 | > | mol = info->nextMolecule(mi)) { |
980 | nIOPerMol[mol->getGlobalIndex()] = mol->getNIntegrableObjects(); | |
981 | } | |
982 | ||
983 | #ifdef IS_MPI | |
984 | std::vector<int> numIntegrableObjectsPerMol(info->getNGlobalMolecules(), 0); | |
985 | MPI_Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], | |
986 | < | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
986 | > | info->getNGlobalMolecules(), MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
987 | > | // MPI::COMM_WORLD.Allreduce(&nIOPerMol[0], &numIntegrableObjectsPerMol[0], |
988 | > | // info->getNGlobalMolecules(), MPI::INT, MPI::SUM); |
989 | #else | |
990 | std::vector<int> numIntegrableObjectsPerMol = nIOPerMol; | |
991 | #endif | |
# | Line 781 | Line 999 | namespace OpenMD { | |
999 | } | |
1000 | ||
1001 | std::vector<StuntDouble*> IOIndexToIntegrableObject(info->getNGlobalIntegrableObjects(), (StuntDouble*)NULL); | |
1002 | < | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
1002 | > | for (mol = info->beginMolecule(mi); mol != NULL; |
1003 | > | mol = info->nextMolecule(mi)) { |
1004 | int myGlobalIndex = mol->getGlobalIndex(); | |
1005 | int globalIO = startingIOIndexForMol[myGlobalIndex]; | |
1006 | < | for (StuntDouble* integrableObject = mol->beginIntegrableObject(ioi); integrableObject != NULL; |
1007 | < | integrableObject = mol->nextIntegrableObject(ioi)) { |
1008 | < | integrableObject->setGlobalIntegrableObjectIndex(globalIO); |
1009 | < | IOIndexToIntegrableObject[globalIO] = integrableObject; |
1006 | > | for (StuntDouble* sd = mol->beginIntegrableObject(ioi); sd != NULL; |
1007 | > | sd = mol->nextIntegrableObject(ioi)) { |
1008 | > | sd->setGlobalIntegrableObjectIndex(globalIO); |
1009 | > | IOIndexToIntegrableObject[globalIO] = sd; |
1010 | globalIO++; | |
1011 | } | |
1012 | } | |
1013 | < | |
1013 | > | |
1014 | info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); | |
1015 | ||
1016 | } | |
1017 | ||
1018 | void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { | |
800 | – | Globals* simParams; |
801 | – | simParams = info->getSimParams(); |
1019 | ||
803 | – | |
1020 | DumpReader reader(info, mdFileName); | |
1021 | int nframes = reader.getNFrames(); | |
1022 | ||
# | Line 814 | Line 1030 | namespace OpenMD { | |
1030 | painCave.isFatal = 1; | |
1031 | simError(); | |
1032 | } | |
817 | – | |
1033 | //copy the current snapshot to previous snapshot | |
1034 | info->getSnapshotManager()->advance(); | |
1035 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |