| 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). |
| 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 |
|
/** |
| 77 |
|
#include "antlr/NoViableAltException.hpp" |
| 78 |
|
|
| 79 |
|
#ifdef IS_MPI |
| 80 |
+ |
#include "mpi.h" |
| 81 |
|
#include "math/ParallelRandNumGen.hpp" |
| 82 |
|
#endif |
| 83 |
|
|
| 84 |
|
namespace OpenMD { |
| 85 |
|
|
| 86 |
< |
Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int startOfMetaDataBlock ){ |
| 86 |
> |
Globals* SimCreator::parseFile(std::istream& rawMetaDataStream, const std::string& filename, int mdFileVersion, int startOfMetaDataBlock ){ |
| 87 |
|
Globals* simParams = NULL; |
| 88 |
|
try { |
| 89 |
|
|
| 94 |
|
const int masterNode = 0; |
| 95 |
|
int commStatus; |
| 96 |
|
if (worldRank == masterNode) { |
| 97 |
< |
#endif |
| 98 |
< |
|
| 97 |
> |
commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
| 98 |
> |
#endif |
| 99 |
|
SimplePreprocessor preprocessor; |
| 100 |
|
preprocessor.preprocess(rawMetaDataStream, filename, startOfMetaDataBlock, ppStream); |
| 101 |
|
|
| 108 |
|
|
| 109 |
|
|
| 110 |
|
} else { |
| 111 |
+ |
|
| 112 |
+ |
commStatus = MPI_Bcast(&mdFileVersion, 1, MPI_INT, masterNode, MPI_COMM_WORLD); |
| 113 |
+ |
|
| 114 |
|
//get stream size |
| 115 |
|
commStatus = MPI_Bcast(&streamSize, 1, MPI_LONG, masterNode, MPI_COMM_WORLD); |
| 116 |
|
|
| 234 |
|
simError(); |
| 235 |
|
} |
| 236 |
|
|
| 237 |
+ |
simParams->setMDfileVersion(mdFileVersion); |
| 238 |
|
return simParams; |
| 239 |
|
} |
| 240 |
|
|
| 241 |
|
SimInfo* SimCreator::createSim(const std::string & mdFileName, |
| 242 |
|
bool loadInitCoords) { |
| 243 |
< |
|
| 243 |
> |
|
| 244 |
|
const int bufferSize = 65535; |
| 245 |
|
char buffer[bufferSize]; |
| 246 |
|
int lineNo = 0; |
| 249 |
|
int metaDataBlockEnd = -1; |
| 250 |
|
int i; |
| 251 |
|
int mdOffset; |
| 252 |
+ |
int mdFileVersion; |
| 253 |
|
|
| 254 |
|
#ifdef IS_MPI |
| 255 |
|
const int masterNode = 0; |
| 283 |
|
painCave.isFatal = 1; |
| 284 |
|
simError(); |
| 285 |
|
} |
| 286 |
+ |
|
| 287 |
+ |
// found the correct opening string, now try to get the file |
| 288 |
+ |
// format version number. |
| 289 |
+ |
|
| 290 |
+ |
StringTokenizer tokenizer(line, "=<> \t\n\r"); |
| 291 |
+ |
std::string fileType = tokenizer.nextToken(); |
| 292 |
+ |
toUpper(fileType); |
| 293 |
|
|
| 294 |
+ |
mdFileVersion = 0; |
| 295 |
+ |
|
| 296 |
+ |
if (fileType == "OPENMD") { |
| 297 |
+ |
while (tokenizer.hasMoreTokens()) { |
| 298 |
+ |
std::string token(tokenizer.nextToken()); |
| 299 |
+ |
toUpper(token); |
| 300 |
+ |
if (token == "VERSION") { |
| 301 |
+ |
mdFileVersion = tokenizer.nextTokenAsInt(); |
| 302 |
+ |
break; |
| 303 |
+ |
} |
| 304 |
+ |
} |
| 305 |
+ |
} |
| 306 |
+ |
|
| 307 |
|
//scan through the input stream and find MetaData tag |
| 308 |
|
while(mdFile_.getline(buffer, bufferSize)) { |
| 309 |
|
++lineNo; |
| 359 |
|
std::stringstream rawMetaDataStream(mdRawData); |
| 360 |
|
|
| 361 |
|
//parse meta-data file |
| 362 |
< |
Globals* simParams = parseFile(rawMetaDataStream, mdFileName, metaDataBlockStart+1); |
| 362 |
> |
Globals* simParams = parseFile(rawMetaDataStream, mdFileName, mdFileVersion, |
| 363 |
> |
metaDataBlockStart + 1); |
| 364 |
|
|
| 365 |
|
//create the force field |
| 366 |
|
ForceField * ff = ForceFieldFactory::getInstance()->createForceField(simParams->getForceField()); |
| 397 |
|
} |
| 398 |
|
|
| 399 |
|
ff->parse(forcefieldFileName); |
| 372 |
– |
ff->setFortranForceOptions(); |
| 400 |
|
//create SimInfo |
| 401 |
|
SimInfo * info = new SimInfo(ff, simParams); |
| 402 |
|
|
| 414 |
|
//create the molecules |
| 415 |
|
createMolecules(info); |
| 416 |
|
|
| 390 |
– |
|
| 417 |
|
//allocate memory for DataStorage(circular reference, need to |
| 418 |
|
//break it) |
| 419 |
|
info->setSnapshotManager(new SimSnapshotManager(info)); |
| 439 |
|
|
| 440 |
|
if (loadInitCoords) |
| 441 |
|
loadCoordinates(info, mdFileName); |
| 416 |
– |
|
| 442 |
|
return info; |
| 443 |
|
} |
| 444 |
|
|
| 637 |
|
#endif |
| 638 |
|
|
| 639 |
|
stampId = info->getMoleculeStampId(i); |
| 640 |
< |
Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId), |
| 641 |
< |
stampId, i, info->getLocalIndexManager()); |
| 640 |
> |
Molecule * mol = molCreator.createMolecule(info->getForceField(), |
| 641 |
> |
info->getMoleculeStamp(stampId), |
| 642 |
> |
stampId, i, |
| 643 |
> |
info->getLocalIndexManager()); |
| 644 |
|
|
| 645 |
|
info->addMolecule(mol); |
| 646 |
|
|
| 667 |
|
int beginRigidBodyIndex; |
| 668 |
|
int beginCutoffGroupIndex; |
| 669 |
|
int nGlobalAtoms = info->getNGlobalAtoms(); |
| 643 |
– |
|
| 644 |
– |
/**@todo fixme */ |
| 645 |
– |
#ifndef IS_MPI |
| 646 |
– |
|
| 647 |
– |
beginAtomIndex = 0; |
| 648 |
– |
beginRigidBodyIndex = 0; |
| 649 |
– |
beginCutoffGroupIndex = 0; |
| 650 |
– |
|
| 651 |
– |
#else |
| 652 |
– |
|
| 653 |
– |
int nproc; |
| 654 |
– |
int myNode; |
| 655 |
– |
|
| 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); |
| 670 |
|
|
| 671 |
|
beginAtomIndex = 0; |
| 672 |
|
beginRigidBodyIndex = 0; |
| 673 |
|
beginCutoffGroupIndex = 0; |
| 674 |
< |
|
| 675 |
< |
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)) { |
| 674 |
> |
|
| 675 |
> |
for(int i = 0; i < info->getNGlobalMolecules(); i++) { |
| 676 |
|
|
| 677 |
< |
//local index(index in DataStorge) of atom is important |
| 678 |
< |
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
| 679 |
< |
atom->setGlobalIndex(beginAtomIndex++); |
| 677 |
> |
#ifdef IS_MPI |
| 678 |
> |
if (info->getMolToProc(i) == worldRank) { |
| 679 |
> |
#endif |
| 680 |
> |
// stuff to do if I own this molecule |
| 681 |
> |
mol = info->getMoleculeByGlobalIndex(i); |
| 682 |
> |
|
| 683 |
> |
//local index(index in DataStorge) of atom is important |
| 684 |
> |
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
| 685 |
> |
atom->setGlobalIndex(beginAtomIndex++); |
| 686 |
> |
} |
| 687 |
> |
|
| 688 |
> |
for(rb = mol->beginRigidBody(ri); rb != NULL; |
| 689 |
> |
rb = mol->nextRigidBody(ri)) { |
| 690 |
> |
rb->setGlobalIndex(beginRigidBodyIndex++); |
| 691 |
> |
} |
| 692 |
> |
|
| 693 |
> |
//local index of cutoff group is trivial, it only depends on |
| 694 |
> |
//the order of travesing |
| 695 |
> |
for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
| 696 |
> |
cg = mol->nextCutoffGroup(ci)) { |
| 697 |
> |
cg->setGlobalIndex(beginCutoffGroupIndex++); |
| 698 |
> |
} |
| 699 |
> |
|
| 700 |
> |
#ifdef IS_MPI |
| 701 |
> |
} else { |
| 702 |
> |
|
| 703 |
> |
// stuff to do if I don't own this molecule |
| 704 |
> |
|
| 705 |
> |
int stampId = info->getMoleculeStampId(i); |
| 706 |
> |
MoleculeStamp* stamp = info->getMoleculeStamp(stampId); |
| 707 |
> |
|
| 708 |
> |
beginAtomIndex += stamp->getNAtoms(); |
| 709 |
> |
beginRigidBodyIndex += stamp->getNRigidBodies(); |
| 710 |
> |
beginCutoffGroupIndex += stamp->getNCutoffGroups() + stamp->getNFreeAtoms(); |
| 711 |
|
} |
| 712 |
< |
|
| 713 |
< |
for(rb = mol->beginRigidBody(ri); rb != NULL; |
| 714 |
< |
rb = mol->nextRigidBody(ri)) { |
| 715 |
< |
rb->setGlobalIndex(beginRigidBodyIndex++); |
| 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 |
< |
|
| 712 |
> |
#endif |
| 713 |
> |
|
| 714 |
> |
} //end for(int i=0) |
| 715 |
> |
|
| 716 |
|
//fill globalGroupMembership |
| 717 |
|
std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0); |
| 718 |
|
for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
| 724 |
|
|
| 725 |
|
} |
| 726 |
|
} |
| 727 |
< |
|
| 727 |
> |
|
| 728 |
|
#ifdef IS_MPI |
| 729 |
|
// Since the globalGroupMembership has been zero filled and we've only |
| 730 |
|
// poked values into the atoms we know, we can do an Allreduce |
| 794 |
|
globalIO++; |
| 795 |
|
} |
| 796 |
|
} |
| 797 |
< |
|
| 797 |
> |
|
| 798 |
|
info->setIOIndexToIntegrableObject(IOIndexToIntegrableObject); |
| 799 |
|
|
| 800 |
|
} |
| 801 |
|
|
| 802 |
|
void SimCreator::loadCoordinates(SimInfo* info, const std::string& mdFileName) { |
| 803 |
|
Globals* simParams; |
| 804 |
+ |
|
| 805 |
|
simParams = info->getSimParams(); |
| 806 |
|
|
| 803 |
– |
|
| 807 |
|
DumpReader reader(info, mdFileName); |
| 808 |
|
int nframes = reader.getNFrames(); |
| 809 |
< |
|
| 809 |
> |
|
| 810 |
|
if (nframes > 0) { |
| 811 |
|
reader.readFrame(nframes - 1); |
| 812 |
|
} else { |
| 817 |
|
painCave.isFatal = 1; |
| 818 |
|
simError(); |
| 819 |
|
} |
| 817 |
– |
|
| 820 |
|
//copy the current snapshot to previous snapshot |
| 821 |
|
info->getSnapshotManager()->advance(); |
| 822 |
|
} |
