--- branches/development/src/io/DumpWriter.cpp 2010/07/09 23:08:25 1465
+++ branches/development/src/io/DumpWriter.cpp 2012/06/10 14:05:02 1752
@@ -36,7 +36,8 @@
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
- * [4] Vardeman & Gezelter, in progress (2009).
+ * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
+ * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
*/
#include "io/DumpWriter.hpp"
@@ -58,8 +59,18 @@ namespace OpenMD {
: info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
Globals* simParams = info->getSimParams();
- needCompression_ = simParams->getCompressDumpFile();
- needForceVector_ = simParams->getOutputForceVector();
+ needCompression_ = simParams->getCompressDumpFile();
+ needForceVector_ = simParams->getOutputForceVector();
+ needParticlePot_ = simParams->getOutputParticlePotential();
+ needFlucQ_ = simParams->getOutputFluctuatingCharges();
+ needElectricField_ = simParams->getOutputElectricField();
+
+ if (needParticlePot_ || needFlucQ_ || needElectricField_) {
+ doSiteData_ = true;
+ } else {
+ doSiteData_ = false;
+ }
+
createDumpFile_ = true;
#ifdef HAVE_LIBZ
if (needCompression_) {
@@ -97,8 +108,18 @@ namespace OpenMD {
Globals* simParams = info->getSimParams();
eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
- needCompression_ = simParams->getCompressDumpFile();
- needForceVector_ = simParams->getOutputForceVector();
+ needCompression_ = simParams->getCompressDumpFile();
+ needForceVector_ = simParams->getOutputForceVector();
+ needParticlePot_ = simParams->getOutputParticlePotential();
+ needFlucQ_ = simParams->getOutputFluctuatingCharges();
+ needElectricField_ = simParams->getOutputElectricField();
+
+ if (needParticlePot_ || needFlucQ_ || needElectricField_) {
+ doSiteData_ = true;
+ } else {
+ doSiteData_ = false;
+ }
+
createDumpFile_ = true;
#ifdef HAVE_LIBZ
if (needCompression_) {
@@ -136,9 +157,18 @@ namespace OpenMD {
Globals* simParams = info->getSimParams();
eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
- needCompression_ = simParams->getCompressDumpFile();
- needForceVector_ = simParams->getOutputForceVector();
-
+ needCompression_ = simParams->getCompressDumpFile();
+ needForceVector_ = simParams->getOutputForceVector();
+ needParticlePot_ = simParams->getOutputParticlePotential();
+ needFlucQ_ = simParams->getOutputFluctuatingCharges();
+ needElectricField_ = simParams->getOutputElectricField();
+
+ if (needParticlePot_ || needFlucQ_ || needElectricField_) {
+ doSiteData_ = true;
+ } else {
+ doSiteData_ = false;
+ }
+
#ifdef HAVE_LIBZ
if (needCompression_) {
filename_ += ".gz";
@@ -272,6 +302,8 @@ namespace OpenMD {
StuntDouble* integrableObject;
SimInfo::MoleculeIterator mi;
Molecule::IntegrableObjectIterator ii;
+ RigidBody::AtomIterator ai;
+ Atom* atom;
#ifndef IS_MPI
os << " \n";
@@ -289,7 +321,32 @@ namespace OpenMD {
}
}
os << " \n";
-
+
+ if (doSiteData_) {
+ os << " \n";
+ for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
+
+ for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
+ integrableObject = mol->nextIntegrableObject(ii)) {
+
+ int ioIndex = integrableObject->getGlobalIntegrableObjectIndex();
+ // do one for the IO itself
+ os << prepareSiteLine(integrableObject, ioIndex, 0);
+
+ if (integrableObject->isRigidBody()) {
+
+ RigidBody* rb = static_cast(integrableObject);
+ int siteIndex = 0;
+ for (atom = rb->beginAtom(ai); atom != NULL;
+ atom = rb->nextAtom(ai)) {
+ os << prepareSiteLine(atom, ioIndex, siteIndex);
+ siteIndex++;
+ }
+ }
+ }
+ }
+ os << " \n";
+ }
os << " \n";
os.flush();
@@ -306,7 +363,8 @@ namespace OpenMD {
}
const int masterNode = 0;
-
+ int nProc;
+ MPI_Comm_size(MPI_COMM_WORLD, &nProc);
if (worldRank == masterNode) {
os << " \n";
writeFrameProperties(os, info_->getSnapshotManager()->getCurrentSnapshot());
@@ -314,13 +372,12 @@ namespace OpenMD {
os << buffer;
- int nProc;
- MPI_Comm_size(MPI_COMM_WORLD, &nProc);
for (int i = 1; i < nProc; ++i) {
// receive the length of the string buffer that was
// prepared by processor i
+ MPI_Bcast(&i, 1, MPI_INT,masterNode,MPI_COMM_WORLD);
int recvLength;
MPI_Recv(&recvLength, 1, MPI_INT, i, 0, MPI_COMM_WORLD, &istatus);
char* recvBuffer = new char[recvLength];
@@ -337,8 +394,14 @@ namespace OpenMD {
os.flush();
} else {
int sendBufferLength = buffer.size() + 1;
- MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
- MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
+ int myturn = 0;
+ for (int i = 1; i < nProc; ++i){
+ MPI_Bcast(&myturn,1, MPI_INT,masterNode,MPI_COMM_WORLD);
+ if (myturn == worldRank){
+ MPI_Send(&sendBufferLength, 1, MPI_INT, masterNode, 0, MPI_COMM_WORLD);
+ MPI_Send((void *)buffer.c_str(), sendBufferLength, MPI_CHAR, masterNode, 0, MPI_COMM_WORLD);
+ }
+ }
}
#endif // is_mpi
@@ -420,10 +483,7 @@ namespace OpenMD {
if (needForceVector_) {
type += "f";
- Vector3d frc;
-
- frc = integrableObject->getFrc();
-
+ Vector3d frc = integrableObject->getFrc();
if (isinf(frc[0]) || isnan(frc[0]) ||
isinf(frc[1]) || isnan(frc[1]) ||
isinf(frc[2]) || isnan(frc[2]) ) {
@@ -439,10 +499,7 @@ namespace OpenMD {
if (integrableObject->isDirectional()) {
type += "t";
- Vector3d trq;
-
- trq = integrableObject->getTrq();
-
+ Vector3d trq = integrableObject->getTrq();
if (isinf(trq[0]) || isnan(trq[0]) ||
isinf(trq[1]) || isnan(trq[1]) ||
isinf(trq[2]) || isnan(trq[2]) ) {
@@ -451,18 +508,109 @@ namespace OpenMD {
" for object %d", index);
painCave.isFatal = 1;
simError();
- }
-
+ }
sprintf(tempBuffer, " %13e %13e %13e",
trq[0], trq[1], trq[2]);
line += tempBuffer;
}
}
-
+
sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
return std::string(tempBuffer);
}
+ std::string DumpWriter::prepareSiteLine(StuntDouble* integrableObject, int ioIndex, int siteIndex) {
+
+
+ std::string id;
+ std::string type;
+ std::string line;
+ char tempBuffer[4096];
+
+ if (integrableObject->isRigidBody()) {
+ sprintf(tempBuffer, "%10d ", ioIndex);
+ id = std::string(tempBuffer);
+ } else {
+ sprintf(tempBuffer, "%10d %10d", ioIndex, siteIndex);
+ id = std::string(tempBuffer);
+ }
+
+ if (needFlucQ_) {
+ type += "cw";
+ RealType fqPos = integrableObject->getFlucQPos();
+ if (isinf(fqPos) || isnan(fqPos) ) {
+ sprintf( painCave.errMsg,
+ "DumpWriter detected a numerical error writing the"
+ " fluctuating charge for object %s", id.c_str());
+ painCave.isFatal = 1;
+ simError();
+ }
+ sprintf(tempBuffer, " %13e ", fqPos);
+ line += tempBuffer;
+
+ RealType fqVel = integrableObject->getFlucQVel();
+ if (isinf(fqVel) || isnan(fqVel) ) {
+ sprintf( painCave.errMsg,
+ "DumpWriter detected a numerical error writing the"
+ " fluctuating charge velocity for object %s", id.c_str());
+ painCave.isFatal = 1;
+ simError();
+ }
+ sprintf(tempBuffer, " %13e ", fqVel);
+ line += tempBuffer;
+
+ if (needForceVector_) {
+ type += "g";
+ RealType fqFrc = integrableObject->getFlucQFrc();
+ if (isinf(fqFrc) || isnan(fqFrc) ) {
+ sprintf( painCave.errMsg,
+ "DumpWriter detected a numerical error writing the"
+ " fluctuating charge force for object %s", id.c_str());
+ painCave.isFatal = 1;
+ simError();
+ }
+ sprintf(tempBuffer, " %13e ", fqFrc);
+ line += tempBuffer;
+ }
+ }
+
+ if (needElectricField_) {
+ type += "e";
+ Vector3d eField= integrableObject->getElectricField();
+ if (isinf(eField[0]) || isnan(eField[0]) ||
+ isinf(eField[1]) || isnan(eField[1]) ||
+ isinf(eField[2]) || isnan(eField[2]) ) {
+ sprintf( painCave.errMsg,
+ "DumpWriter detected a numerical error writing the electric"
+ " field for object %s", id.c_str());
+ painCave.isFatal = 1;
+ simError();
+ }
+ sprintf(tempBuffer, " %13e %13e %13e",
+ eField[0], eField[1], eField[2]);
+ line += tempBuffer;
+ }
+
+
+ if (needParticlePot_) {
+ type += "u";
+ RealType particlePot = integrableObject->getParticlePot();
+ if (isinf(particlePot) || isnan(particlePot)) {
+ sprintf( painCave.errMsg,
+ "DumpWriter detected a numerical error writing the particle "
+ " potential for object %s", id.c_str());
+ painCave.isFatal = 1;
+ simError();
+ }
+ sprintf(tempBuffer, " %13e", particlePot);
+ line += tempBuffer;
+ }
+
+
+ sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
+ return std::string(tempBuffer);
+ }
+
void DumpWriter::writeDump() {
writeFrame(*dumpFile_);
}
@@ -540,7 +688,7 @@ namespace OpenMD {
newOStream = new std::ofstream(filename.c_str());
#endif
//write out MetaData first
- (*newOStream) << "" << std::endl;
+ (*newOStream) << "" << std::endl;
(*newOStream) << " " << std::endl;
(*newOStream) << info_->getRawMetaData();
(*newOStream) << " " << std::endl;