--- branches/development/src/io/DumpWriter.cpp 2011/09/14 21:15:17 1629
+++ 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,9 +59,18 @@ namespace OpenMD {
: info_(info), filename_(info->getDumpFileName()), eorFilename_(info->getFinalConfigFileName()){
Globals* simParams = info->getSimParams();
- needCompression_ = simParams->getCompressDumpFile();
- needForceVector_ = simParams->getOutputForceVector();
- needParticlePot_ = simParams->getOutputParticlePotential();
+ 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_) {
@@ -98,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_) {
@@ -137,10 +157,18 @@ namespace OpenMD {
Globals* simParams = info->getSimParams();
eorFilename_ = filename_.substr(0, filename_.rfind(".")) + ".eor";
- needCompression_ = simParams->getCompressDumpFile();
- needForceVector_ = simParams->getOutputForceVector();
- needParticlePot_ = simParams->getOutputParticlePotential();
-
+ 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";
@@ -274,6 +302,8 @@ namespace OpenMD {
StuntDouble* integrableObject;
SimInfo::MoleculeIterator mi;
Molecule::IntegrableObjectIterator ii;
+ RigidBody::AtomIterator ai;
+ Atom* atom;
#ifndef IS_MPI
os << " \n";
@@ -291,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();
@@ -428,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]) ) {
@@ -447,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]) ) {
@@ -459,23 +508,97 @@ namespace OpenMD {
" for object %d", index);
painCave.isFatal = 1;
simError();
- }
-
+ }
sprintf(tempBuffer, " %13e %13e %13e",
trq[0], trq[1], trq[2]);
line += tempBuffer;
}
}
- if (needParticlePot_) {
- type += "u";
- RealType particlePot;
- particlePot = integrableObject->getParticlePot();
+ 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 %d", index);
+ " potential for object %s", id.c_str());
painCave.isFatal = 1;
simError();
}
@@ -483,7 +606,8 @@ namespace OpenMD {
line += tempBuffer;
}
- sprintf(tempBuffer, "%10d %7s %s\n", index, type.c_str(), line.c_str());
+
+ sprintf(tempBuffer, "%s %7s %s\n", id.c_str(), type.c_str(), line.c_str());
return std::string(tempBuffer);
}
@@ -564,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;