[ViewVC] Diff of: OpenMD/trunk/src/rnemd/RNEMD.cpp

Comparing:
branches/development/src/rnemd/RNEMD.cpp (file contents), Revision 1773 by gezelter, Tue Aug 7 18:26:40 2012 UTC vs.
trunk/src/rnemd/RNEMD.cpp (file contents), Revision 1801 by gezelter, Mon Oct 1 18:21:15 2012 UTC

#	Line 53 \| Line 53
53		#include <mpi.h>
54		#endif
55
56	+	#ifdef _MSC_VER
57	+	#define isnan(x) _isnan((x))
58	+	#define isinf(x) (!_finite(x) && !_isnan(x))
59	+	#endif
60	+
61		#define HONKING_LARGE_VALUE 1.0e10
62
63		using namespace std;
#	Line 65 \| Line 70 \| namespace OpenMD {
70		failTrialCount_ = 0;
71		failRootCount_ = 0;
72
68	–	int seedValue;
73		Globals * simParams = info->getSimParams();
74		RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
75
76	+	doRNEMD_ = rnemdParams->getUseRNEMD();
77	+	if (!doRNEMD_) return;
78	+
79		stringToMethod_["Swap"] = rnemdSwap;
80		stringToMethod_["NIVS"] = rnemdNIVS;
81		stringToMethod_["VSS"] = rnemdVSS;
#	Line 77 \| Line 84 \| namespace OpenMD {
84		stringToFluxType_["Px"] = rnemdPx;
85		stringToFluxType_["Py"] = rnemdPy;
86		stringToFluxType_["Pz"] = rnemdPz;
87	+	stringToFluxType_["Pvector"] = rnemdPvector;
88		stringToFluxType_["KE+Px"] = rnemdKePx;
89		stringToFluxType_["KE+Py"] = rnemdKePy;
90		stringToFluxType_["KE+Pvector"] = rnemdKePvector;
#	Line 98 \| Line 106 \| namespace OpenMD {
106		sprintf(painCave.errMsg,
107		"RNEMD: No fluxType was set in the md file. This parameter,\n"
108		"\twhich must be one of the following values:\n"
109	<	"\tKE, Px, Py, Pz, KE+Px, KE+Py, KE+Pvector, must be set to\n"
110	<	"\tuse RNEMD\n");
109	>	"\tKE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector\n"
110	>	"\tmust be set to use RNEMD\n");
111		painCave.isFatal = 1;
112		painCave.severity = OPENMD_ERROR;
113		simError();
#	Line 197 \| Line 205 \| namespace OpenMD {
205		break;
206		case rnemdPvector:
207		hasCorrectFlux = hasMomentumFluxVector;
208	+	break;
209		case rnemdKePx:
210		case rnemdKePy:
211		hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
#	Line 224 \| Line 233 \| namespace OpenMD {
233		}
234		if (!hasCorrectFlux) {
235		sprintf(painCave.errMsg,
236	<	"RNEMD: The current method,\n"
228	<	"\t%s, and flux type %s\n"
236	>	"RNEMD: The current method, %s, and flux type, %s,\n"
237		"\tdid not have the correct flux value specified. Options\n"
238		"\tinclude: kineticFlux, momentumFlux, and momentumFluxVector\n",
239		methStr.c_str(), fluxStr.c_str());
#	Line 235 \| Line 243 \| namespace OpenMD {
243		}
244
245		if (hasKineticFlux) {
246	<	kineticFlux_ = rnemdParams->getKineticFlux();
246	>	// convert the kcal / mol / Angstroms^2 / fs values in the md file
247	>	// into amu / fs^3:
248	>	kineticFlux_ = rnemdParams->getKineticFlux()
249	>	* PhysicalConstants::energyConvert;
250		} else {
251		kineticFlux_ = 0.0;
252		}
#	Line 270 \| Line 281 \| namespace OpenMD {
281		// do some sanity checking
282
283		int selectionCount = seleMan_.getSelectionCount();
284	+
285		int nIntegrable = info->getNGlobalIntegrableObjects();
286
287		if (selectionCount > nIntegrable) {
#	Line 288 \| Line 300 \| namespace OpenMD {
300		simError();
301		}
302
303	+	areaAccumulator_ = new Accumulator();
304	+
305		nBins_ = rnemdParams->getOutputBins();
306
307		data_.resize(RNEMD::ENDINDEX);
#	Line 296 \| Line 310 \| namespace OpenMD {
310		z.title = "Z";
311		z.dataType = "RealType";
312		z.accumulator.reserve(nBins_);
313	<	for (unsigned int i = 0; i < nBins_; i++)
313	>	for (int i = 0; i < nBins_; i++)
314		z.accumulator.push_back( new Accumulator() );
315		data_[Z] = z;
316		outputMap_["Z"] = Z;
#	Line 306 \| Line 320 \| namespace OpenMD {
320		temperature.title = "Temperature";
321		temperature.dataType = "RealType";
322		temperature.accumulator.reserve(nBins_);
323	<	for (unsigned int i = 0; i < nBins_; i++)
323	>	for (int i = 0; i < nBins_; i++)
324		temperature.accumulator.push_back( new Accumulator() );
325		data_[TEMPERATURE] = temperature;
326		outputMap_["TEMPERATURE"] = TEMPERATURE;
327
328		OutputData velocity;
329	<	velocity.units = "amu/fs";
329	>	velocity.units = "angstroms/fs";
330		velocity.title = "Velocity";
331		velocity.dataType = "Vector3d";
332		velocity.accumulator.reserve(nBins_);
333	<	for (unsigned int i = 0; i < nBins_; i++)
333	>	for (int i = 0; i < nBins_; i++)
334		velocity.accumulator.push_back( new VectorAccumulator() );
335		data_[VELOCITY] = velocity;
336		outputMap_["VELOCITY"] = VELOCITY;
#	Line 326 \| Line 340 \| namespace OpenMD {
340		density.title = "Density";
341		density.dataType = "RealType";
342		density.accumulator.reserve(nBins_);
343	<	for (unsigned int i = 0; i < nBins_; i++)
343	>	for (int i = 0; i < nBins_; i++)
344		density.accumulator.push_back( new Accumulator() );
345		data_[DENSITY] = density;
346		outputMap_["DENSITY"] = DENSITY;
#	Line 381 \| Line 395 \| namespace OpenMD {
395		// dt = exchange time interval
396		// flux = target flux
397
398	<	kineticTarget_ = 2.0kineticFlux_exchangeTime_hmat(0,0)hmat(1,1);
399	<	momentumTarget_ = 2.0momentumFluxVector_exchangeTime_hmat(0,0)hmat(1,1);
398	>	RealType area = currentSnap_->getXYarea();
399	>	kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
400	>	momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
401
402		// total exchange sums are zeroed out at the beginning:
403
#	Line 407 \| Line 422 \| namespace OpenMD {
422		}
423
424		RNEMD::~RNEMD() {
425	<
425	>	if (!doRNEMD_) return;
426		#ifdef IS_MPI
427		if (worldRank == 0) {
428		#endif
#	Line 429 \| Line 444 \| namespace OpenMD {
444		}
445
446		void RNEMD::doSwap() {
447	<
447	>	if (!doRNEMD_) return;
448		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
449		Mat3x3d hmat = currentSnap_->getHmat();
450
#	Line 437 \| Line 452 \| namespace OpenMD {
452
453		int selei;
454		StuntDouble* sd;
440	–	int idx;
455
456		RealType min_val;
457		bool min_found = false;
#	Line 450 \| Line 464 \| namespace OpenMD {
464		for (sd = seleMan_.beginSelected(selei); sd != NULL;
465		sd = seleMan_.nextSelected(selei)) {
466
453	–	idx = sd->getLocalIndex();
454	–
467		Vector3d pos = sd->getPos();
468
469		// wrap the stuntdouble's position back into the box:
#	Line 488 \| Line 500 \| namespace OpenMD {
500		+ angMom[2]*angMom[2]/I(2, 2);
501		}
502		} //angular momenta exchange enabled
491	–	//energyConvert temporarily disabled
492	–	//make kineticExchange_ comparable between swap & scale
493	–	//value = value * 0.5 / PhysicalConstants::energyConvert;
503		value *= 0.5;
504		break;
505		case rnemdPx :
#	Line 532 \| Line 541 \| namespace OpenMD {
541		}
542		}
543
544	<	#ifdef IS_MPI
545	<	int nProc, worldRank;
544	>	#ifdef IS_MPI
545	>	int worldRank = MPI::COMM_WORLD.Get_rank();
546
538	–	nProc = MPI::COMM_WORLD.Get_size();
539	–	worldRank = MPI::COMM_WORLD.Get_rank();
540	–
547		bool my_min_found = min_found;
548		bool my_max_found = max_found;
549
#	Line 728 \| Line 734 \| namespace OpenMD {
734
735		switch(rnemdFluxType_) {
736		case rnemdKE:
731	–	cerr << "KE\n";
737		kineticExchange_ += max_val - min_val;
738		break;
739		case rnemdPx:
#	Line 741 \| Line 746 \| namespace OpenMD {
746		momentumExchange_.z() += max_val - min_val;
747		break;
748		default:
744	–	cerr << "default\n";
749		break;
750		}
751		} else {
#	Line 764 \| Line 768 \| namespace OpenMD {
768		}
769
770		void RNEMD::doNIVS() {
771	<
771	>	if (!doRNEMD_) return;
772		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
773		Mat3x3d hmat = currentSnap_->getHmat();
774
#	Line 772 \| Line 776 \| namespace OpenMD {
776
777		int selei;
778		StuntDouble* sd;
775	–	int idx;
779
780		vector<StuntDouble*> hotBin, coldBin;
781
#	Line 794 \| Line 797 \| namespace OpenMD {
797		for (sd = seleMan_.beginSelected(selei); sd != NULL;
798		sd = seleMan_.nextSelected(selei)) {
799
797	–	idx = sd->getLocalIndex();
798	–
800		Vector3d pos = sd->getPos();
801
802		// wrap the stuntdouble's position back into the box:
#	Line 908 \| Line 909 \| namespace OpenMD {
909
910		if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients
911		c = sqrt(c);
912	<	//std::cerr << "cold slab scaling coefficient: " << c << endl;
912	<	//now convert to hotBin coefficient
912	>
913		RealType w = 0.0;
914		if (rnemdFluxType_ == rnemdFullKE) {
915		x = 1.0 + px * (1.0 - c);
#	Line 947 \| Line 947 \| namespace OpenMD {
947		}
948		}
949		w = sqrt(w);
950	–	// std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
951	–	// << "\twh= " << w << endl;
950		for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
951		if (rnemdFluxType_ == rnemdFullKE) {
952		vel = (*sdi)->getVel();
#	Line 1213 \| Line 1211 \| namespace OpenMD {
1211		}
1212
1213		void RNEMD::doVSS() {
1214	<
1214	>	if (!doRNEMD_) return;
1215		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1216		RealType time = currentSnap_->getTime();
1217		Mat3x3d hmat = currentSnap_->getHmat();
#	Line 1222 \| Line 1220 \| namespace OpenMD {
1220
1221		int selei;
1222		StuntDouble* sd;
1225	–	int idx;
1223
1224		vector<StuntDouble*> hotBin, coldBin;
1225
#	Line 1237 \| Line 1234 \| namespace OpenMD {
1234		for (sd = seleMan_.beginSelected(selei); sd != NULL;
1235		sd = seleMan_.nextSelected(selei)) {
1236
1240	–	idx = sd->getLocalIndex();
1241	–
1237		Vector3d pos = sd->getPos();
1238
1239		// wrap the stuntdouble's position back into the box:
#	Line 1257 \| Line 1252 \| namespace OpenMD {
1252
1253		if (inA) {
1254		hotBin.push_back(sd);
1260	–	//std::cerr << "before, velocity = " << vel << endl;
1255		Ph += mass * vel;
1262	–	//std::cerr << "after, velocity = " << vel << endl;
1256		Mh += mass;
1257		Kh += mass * vel.lengthSquare();
1258		if (rnemdFluxType_ == rnemdFullKE) {
#	Line 1307 \| Line 1300 \| namespace OpenMD {
1300
1301		Kh *= 0.5;
1302		Kc *= 0.5;
1310	–
1311	–	// std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1312	–	// << "\tKc= " << Kc << endl;
1313	–	// std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1303
1304		#ifdef IS_MPI
1305		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
#	Line 1342 \| Line 1331 \| namespace OpenMD {
1331		if (hDenominator > 0.0) {
1332		RealType h = sqrt(hNumerator / hDenominator);
1333		if ((h > 0.9) && (h < 1.1)) {
1334	<	// std::cerr << "cold slab scaling coefficient: " << c << "\n";
1346	<	// std::cerr << "hot slab scaling coefficient: " << h << "\n";
1334	>
1335		vector<StuntDouble*>::iterator sdi;
1336		Vector3d vel;
1337		for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
#	Line 1391 \| Line 1379 \| namespace OpenMD {
1379		}
1380
1381		void RNEMD::doRNEMD() {
1382	<
1382	>	if (!doRNEMD_) return;
1383		trialCount_++;
1384		switch(rnemdMethod_) {
1385		case rnemdSwap:
#	Line 1410 \| Line 1398 \| namespace OpenMD {
1398		}
1399
1400		void RNEMD::collectData() {
1401	<
1401	>	if (!doRNEMD_) return;
1402		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1403		Mat3x3d hmat = currentSnap_->getHmat();
1404
1405	+	areaAccumulator_->add(currentSnap_->getXYarea());
1406	+
1407		seleMan_.setSelectionSet(evaluator_.evaluate());
1408
1409	<	int selei;
1409	>	int selei(0);
1410		StuntDouble* sd;
1421	–	int idx;
1411
1412		vector<RealType> binMass(nBins_, 0.0);
1413		vector<RealType> binPx(nBins_, 0.0);
#	Line 1442 \| Line 1431 \| namespace OpenMD {
1431		sd != NULL;
1432		sd = mol->nextIntegrableObject(iiter))
1433		*/
1434	+
1435		for (sd = seleMan_.beginSelected(selei); sd != NULL;
1436	<	sd = seleMan_.nextSelected(selei)) {
1437	<
1448	<	idx = sd->getLocalIndex();
1449	<
1436	>	sd = seleMan_.nextSelected(selei)) {
1437	>
1438		Vector3d pos = sd->getPos();
1439
1440		// wrap the stuntdouble's position back into the box:
#	Line 1454 \| Line 1442 \| namespace OpenMD {
1442		if (usePeriodicBoundaryConditions_)
1443		currentSnap_->wrapVector(pos);
1444
1445	+
1446		// which bin is this stuntdouble in?
1447		// wrapped positions are in the range [-0.5hmat(2,2), +0.5hmat(2,2)]
1448		// Shift molecules by half a box to have bins start at 0
1449		// The modulo operator is used to wrap the case when we are
1450		// beyond the end of the bins back to the beginning.
1451		int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
1452	<
1452	>
1453		RealType mass = sd->getMass();
1454		Vector3d vel = sd->getVel();
1455
#	Line 1491 \| Line 1480 \| namespace OpenMD {
1480		}
1481		}
1482
1494	–
1483		#ifdef IS_MPI
1484		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1485		nBins_, MPI::INT, MPI::SUM);
#	Line 1518 \| Line 1506 \| namespace OpenMD {
1506		vel.x() = binPx[i] / binMass[i];
1507		vel.y() = binPy[i] / binMass[i];
1508		vel.z() = binPz[i] / binMass[i];
1509	<	den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1509	>
1510	>	den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1511	>	/ currentSnap_->getVolume() ;
1512	>
1513		temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1514		PhysicalConstants::energyConvert);
1515	<
1515	>
1516		for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1517		if(outputMask_[j]) {
1518		switch(j) {
1519		case Z:
1520	<	(data_[j].accumulator[i])->add(z);
1520	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1521		break;
1522		case TEMPERATURE:
1523	<	data_[j].accumulator[i]->add(temp);
1523	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1524		break;
1525		case VELOCITY:
1526		dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1527		break;
1528		case DENSITY:
1529	<	data_[j].accumulator[i]->add(den);
1529	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1530		break;
1531		}
1532		}
#	Line 1544 \| Line 1535 \| namespace OpenMD {
1535		}
1536
1537		void RNEMD::getStarted() {
1538	+	if (!doRNEMD_) return;
1539		collectData();
1540		writeOutputFile();
1541		}
1542
1543		void RNEMD::parseOutputFileFormat(const std::string& format) {
1544	+	if (!doRNEMD_) return;
1545		StringTokenizer tokenizer(format, " ,;\|\t\n\r");
1546
1547		while(tokenizer.hasMoreTokens()) {
#	Line 1569 \| Line 1562 \| namespace OpenMD {
1562		}
1563
1564		void RNEMD::writeOutputFile() {
1565	+	if (!doRNEMD_) return;
1566
1567		#ifdef IS_MPI
1568		// If we're the root node, should we print out the results
#	Line 1588 \| Line 1582 \| namespace OpenMD {
1582		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1583
1584		RealType time = currentSnap_->getTime();
1585	<
1586	<
1585	>	RealType avgArea;
1586	>	areaAccumulator_->getAverage(avgArea);
1587	>	RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1588	>	/ PhysicalConstants::energyConvert;
1589	>	Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);
1590	>
1591		rnemdFile_ << "#######################################################\n";
1592		rnemdFile_ << "# RNEMD {\n";
1593
1594		map<string, RNEMDMethod>::iterator mi;
1595		for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1596		if ( (*mi).second == rnemdMethod_)
1597	<	rnemdFile_ << "# exchangeMethod = " << (*mi).first << "\n";
1597	>	rnemdFile_ << "# exchangeMethod = \"" << (*mi).first << "\";\n";
1598		}
1599		map<string, RNEMDFluxType>::iterator fi;
1600		for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1601		if ( (*fi).second == rnemdFluxType_)
1602	<	rnemdFile_ << "# fluxType = " << (*fi).first << "\n";
1602	>	rnemdFile_ << "# fluxType = \"" << (*fi).first << "\";\n";
1603		}
1604
1605	<	rnemdFile_ << "# exchangeTime = " << exchangeTime_ << " fs\n";
1605	>	rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n";
1606
1607		rnemdFile_ << "# objectSelection = \""
1608	<	<< rnemdObjectSelection_ << "\"\n";
1609	<	rnemdFile_ << "# slabWidth = " << slabWidth_ << " angstroms\n";
1610	<	rnemdFile_ << "# slabAcenter = " << slabACenter_ << " angstroms\n";
1611	<	rnemdFile_ << "# slabBcenter = " << slabBCenter_ << " angstroms\n";
1608	>	<< rnemdObjectSelection_ << "\";\n";
1609	>	rnemdFile_ << "# slabWidth = " << slabWidth_ << ";\n";
1610	>	rnemdFile_ << "# slabAcenter = " << slabACenter_ << ";\n";
1611	>	rnemdFile_ << "# slabBcenter = " << slabBCenter_ << ";\n";
1612		rnemdFile_ << "# }\n";
1613		rnemdFile_ << "#######################################################\n";
1614	<
1615	<	rnemdFile_ << "# running time = " << time << " fs\n";
1616	<	rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1617	<	rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1618	<
1619	<	rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1620	<	<< "\n";
1621	<	rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1622	<	<< "\n";
1623	<
1624	<	rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1625	<	rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1626	<	<< "\n";
1627	<
1628	<	rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1629	<	rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1630	<
1631	<
1614	>	rnemdFile_ << "# RNEMD report:\n";
1615	>	rnemdFile_ << "# running time = " << time << " fs\n";
1616	>	rnemdFile_ << "# target flux:\n";
1617	>	rnemdFile_ << "# kinetic = "
1618	>	<< kineticFlux_ / PhysicalConstants::energyConvert
1619	>	<< " (kcal/mol/A^2/fs)\n";
1620	>	rnemdFile_ << "# momentum = " << momentumFluxVector_
1621	>	<< " (amu/A/fs^2)\n";
1622	>	rnemdFile_ << "# target one-time exchanges:\n";
1623	>	rnemdFile_ << "# kinetic = "
1624	>	<< kineticTarget_ / PhysicalConstants::energyConvert
1625	>	<< " (kcal/mol)\n";
1626	>	rnemdFile_ << "# momentum = " << momentumTarget_
1627	>	<< " (amu*A/fs)\n";
1628	>	rnemdFile_ << "# actual exchange totals:\n";
1629	>	rnemdFile_ << "# kinetic = "
1630	>	<< kineticExchange_ / PhysicalConstants::energyConvert
1631	>	<< " (kcal/mol)\n";
1632	>	rnemdFile_ << "# momentum = " << momentumExchange_
1633	>	<< " (amu*A/fs)\n";
1634	>	rnemdFile_ << "# actual flux:\n";
1635	>	rnemdFile_ << "# kinetic = " << Jz
1636	>	<< " (kcal/mol/A^2/fs)\n";
1637	>	rnemdFile_ << "# momentum = " << JzP
1638	>	<< " (amu/A/fs^2)\n";
1639	>	rnemdFile_ << "# exchange statistics:\n";
1640	>	rnemdFile_ << "# attempted = " << trialCount_ << "\n";
1641	>	rnemdFile_ << "# failed = " << failTrialCount_ << "\n";
1642		if (rnemdMethod_ == rnemdNIVS) {
1643	<	rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1643	>	rnemdFile_ << "# NIVS root-check errors = "
1644	>	<< failRootCount_ << "\n";
1645		}
1637	–
1646		rnemdFile_ << "#######################################################\n";
1647
1648
#	Line 1645 \| Line 1653 \| namespace OpenMD {
1653		if (outputMask_[i]) {
1654		rnemdFile_ << "\t" << data_[i].title <<
1655		"(" << data_[i].units << ")";
1656	+	// add some extra tabs for column alignment
1657	+	if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1658		}
1659		}
1660		rnemdFile_ << std::endl;
1661
1662		rnemdFile_.precision(8);
1663
1664	<	for (unsigned int j = 0; j < nBins_; j++) {
1664	>	for (int j = 0; j < nBins_; j++) {
1665
1666		for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1667		if (outputMask_[i]) {
#	Line 1671 \| Line 1681 \| namespace OpenMD {
1681		rnemdFile_ << std::endl;
1682
1683		}
1684	+
1685	+	rnemdFile_ << "#######################################################\n";
1686	+	rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1687	+	rnemdFile_ << "#######################################################\n";
1688	+
1689	+
1690	+	for (int j = 0; j < nBins_; j++) {
1691	+	rnemdFile_ << "#";
1692	+	for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1693	+	if (outputMask_[i]) {
1694	+	if (data_[i].dataType == "RealType")
1695	+	writeRealStdDev(i,j);
1696	+	else if (data_[i].dataType == "Vector3d")
1697	+	writeVectorStdDev(i,j);
1698	+	else {
1699	+	sprintf( painCave.errMsg,
1700	+	"RNEMD found an unknown data type for: %s ",
1701	+	data_[i].title.c_str());
1702	+	painCave.isFatal = 1;
1703	+	simError();
1704	+	}
1705	+	}
1706	+	}
1707	+	rnemdFile_ << std::endl;
1708	+
1709	+	}
1710
1711		rnemdFile_.flush();
1712		rnemdFile_.close();
#	Line 1682 \| Line 1718 \| namespace OpenMD {
1718		}
1719
1720		void RNEMD::writeReal(int index, unsigned int bin) {
1721	+	if (!doRNEMD_) return;
1722		assert(index >=0 && index < ENDINDEX);
1723	<	assert(bin >=0 && bin < nBins_);
1723	>	assert(bin < nBins_);
1724		RealType s;
1725
1726	<	data_[index].accumulator[bin]->getAverage(s);
1726	>	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
1727
1728		if (! isinf(s) && ! isnan(s)) {
1729		rnemdFile_ << "\t" << s;
#	Line 1700 \| Line 1737 \| namespace OpenMD {
1737		}
1738
1739		void RNEMD::writeVector(int index, unsigned int bin) {
1740	+	if (!doRNEMD_) return;
1741		assert(index >=0 && index < ENDINDEX);
1742	<	assert(bin >=0 && bin < nBins_);
1742	>	assert(bin < nBins_);
1743		Vector3d s;
1744		dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1745		if (isinf(s[0]) \|\| isnan(s[0]) \|\|
#	Line 1716 \| Line 1754 \| namespace OpenMD {
1754		rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1755		}
1756		}
1757	+
1758	+	void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1759	+	if (!doRNEMD_) return;
1760	+	assert(index >=0 && index < ENDINDEX);
1761	+	assert(bin < nBins_);
1762	+	RealType s;
1763	+
1764	+	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
1765	+
1766	+	if (! isinf(s) && ! isnan(s)) {
1767	+	rnemdFile_ << "\t" << s;
1768	+	} else{
1769	+	sprintf( painCave.errMsg,
1770	+	"RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1771	+	data_[index].title.c_str(), bin);
1772	+	painCave.isFatal = 1;
1773	+	simError();
1774	+	}
1775	+	}
1776	+
1777	+	void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1778	+	if (!doRNEMD_) return;
1779	+	assert(index >=0 && index < ENDINDEX);
1780	+	assert(bin < nBins_);
1781	+	Vector3d s;
1782	+	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1783	+	if (isinf(s[0]) \|\| isnan(s[0]) \|\|
1784	+	isinf(s[1]) \|\| isnan(s[1]) \|\|
1785	+	isinf(s[2]) \|\| isnan(s[2]) ) {
1786	+	sprintf( painCave.errMsg,
1787	+	"RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1788	+	data_[index].title.c_str(), bin);
1789	+	painCave.isFatal = 1;
1790	+	simError();
1791	+	} else {
1792	+	rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1793	+	}
1794	+	}
1795		}
1796

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing: branches/development/src/rnemd/RNEMD.cpp (file contents), Revision 1773 by gezelter, Tue Aug 7 18:26:40 2012 UTC vs. trunk/src/rnemd/RNEMD.cpp (file contents), Revision 1801 by gezelter, Mon Oct 1 18:21:15 2012 UTC

Diff Legend

Comparing:
branches/development/src/rnemd/RNEMD.cpp (file contents), Revision 1773 by gezelter, Tue Aug 7 18:26:40 2012 UTC vs.
trunk/src/rnemd/RNEMD.cpp (file contents), Revision 1801 by gezelter, Mon Oct 1 18:21:15 2012 UTC