[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 1793 by gezelter, Fri Aug 31 21:16:10 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 288 \| Line 299 \| namespace OpenMD {
299		simError();
300		}
301
302	+	areaAccumulator_ = new Accumulator();
303	+
304		nBins_ = rnemdParams->getOutputBins();
305
306		data_.resize(RNEMD::ENDINDEX);
#	Line 296 \| Line 309 \| namespace OpenMD {
309		z.title = "Z";
310		z.dataType = "RealType";
311		z.accumulator.reserve(nBins_);
312	<	for (unsigned int i = 0; i < nBins_; i++)
312	>	for (int i = 0; i < nBins_; i++)
313		z.accumulator.push_back( new Accumulator() );
314		data_[Z] = z;
315		outputMap_["Z"] = Z;
#	Line 306 \| Line 319 \| namespace OpenMD {
319		temperature.title = "Temperature";
320		temperature.dataType = "RealType";
321		temperature.accumulator.reserve(nBins_);
322	<	for (unsigned int i = 0; i < nBins_; i++)
322	>	for (int i = 0; i < nBins_; i++)
323		temperature.accumulator.push_back( new Accumulator() );
324		data_[TEMPERATURE] = temperature;
325		outputMap_["TEMPERATURE"] = TEMPERATURE;
326
327		OutputData velocity;
328	<	velocity.units = "amu/fs";
328	>	velocity.units = "angstroms/fs";
329		velocity.title = "Velocity";
330		velocity.dataType = "Vector3d";
331		velocity.accumulator.reserve(nBins_);
332	<	for (unsigned int i = 0; i < nBins_; i++)
332	>	for (int i = 0; i < nBins_; i++)
333		velocity.accumulator.push_back( new VectorAccumulator() );
334		data_[VELOCITY] = velocity;
335		outputMap_["VELOCITY"] = VELOCITY;
#	Line 326 \| Line 339 \| namespace OpenMD {
339		density.title = "Density";
340		density.dataType = "RealType";
341		density.accumulator.reserve(nBins_);
342	<	for (unsigned int i = 0; i < nBins_; i++)
342	>	for (int i = 0; i < nBins_; i++)
343		density.accumulator.push_back( new Accumulator() );
344		data_[DENSITY] = density;
345		outputMap_["DENSITY"] = DENSITY;
#	Line 381 \| Line 394 \| namespace OpenMD {
394		// dt = exchange time interval
395		// flux = target flux
396
397	<	kineticTarget_ = 2.0kineticFlux_exchangeTime_hmat(0,0)hmat(1,1);
398	<	momentumTarget_ = 2.0momentumFluxVector_exchangeTime_hmat(0,0)hmat(1,1);
397	>	RealType area = currentSnap_->getXYarea();
398	>	kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
399	>	momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
400
401		// total exchange sums are zeroed out at the beginning:
402
#	Line 407 \| Line 421 \| namespace OpenMD {
421		}
422
423		RNEMD::~RNEMD() {
424	<
424	>	if (!doRNEMD_) return;
425		#ifdef IS_MPI
426		if (worldRank == 0) {
427		#endif
#	Line 429 \| Line 443 \| namespace OpenMD {
443		}
444
445		void RNEMD::doSwap() {
446	<
446	>	if (!doRNEMD_) return;
447		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
448		Mat3x3d hmat = currentSnap_->getHmat();
449
#	Line 437 \| Line 451 \| namespace OpenMD {
451
452		int selei;
453		StuntDouble* sd;
440	–	int idx;
454
455		RealType min_val;
456		bool min_found = false;
#	Line 449 \| Line 462 \| namespace OpenMD {
462
463		for (sd = seleMan_.beginSelected(selei); sd != NULL;
464		sd = seleMan_.nextSelected(selei)) {
452	–
453	–	idx = sd->getLocalIndex();
465
466		Vector3d pos = sd->getPos();
467
#	Line 488 \| Line 499 \| namespace OpenMD {
499		+ angMom[2]*angMom[2]/I(2, 2);
500		}
501		} //angular momenta exchange enabled
491	–	//energyConvert temporarily disabled
492	–	//make kineticExchange_ comparable between swap & scale
493	–	//value = value * 0.5 / PhysicalConstants::energyConvert;
502		value *= 0.5;
503		break;
504		case rnemdPx :
#	Line 532 \| Line 540 \| namespace OpenMD {
540		}
541		}
542
543	<	#ifdef IS_MPI
544	<	int nProc, worldRank;
543	>	#ifdef IS_MPI
544	>	int worldRank = MPI::COMM_WORLD.Get_rank();
545
538	–	nProc = MPI::COMM_WORLD.Get_size();
539	–	worldRank = MPI::COMM_WORLD.Get_rank();
540	–
546		bool my_min_found = min_found;
547		bool my_max_found = max_found;
548
#	Line 728 \| Line 733 \| namespace OpenMD {
733
734		switch(rnemdFluxType_) {
735		case rnemdKE:
731	–	cerr << "KE\n";
736		kineticExchange_ += max_val - min_val;
737		break;
738		case rnemdPx:
#	Line 741 \| Line 745 \| namespace OpenMD {
745		momentumExchange_.z() += max_val - min_val;
746		break;
747		default:
744	–	cerr << "default\n";
748		break;
749		}
750		} else {
#	Line 764 \| Line 767 \| namespace OpenMD {
767		}
768
769		void RNEMD::doNIVS() {
770	<
770	>	if (!doRNEMD_) return;
771		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
772		Mat3x3d hmat = currentSnap_->getHmat();
773
#	Line 772 \| Line 775 \| namespace OpenMD {
775
776		int selei;
777		StuntDouble* sd;
775	–	int idx;
778
779		vector<StuntDouble*> hotBin, coldBin;
780
#	Line 793 \| Line 795 \| namespace OpenMD {
795
796		for (sd = seleMan_.beginSelected(selei); sd != NULL;
797		sd = seleMan_.nextSelected(selei)) {
796	–
797	–	idx = sd->getLocalIndex();
798
799		Vector3d pos = sd->getPos();
800
#	Line 1213 \| Line 1213 \| namespace OpenMD {
1213		}
1214
1215		void RNEMD::doVSS() {
1216	<
1216	>	if (!doRNEMD_) return;
1217		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1218		RealType time = currentSnap_->getTime();
1219		Mat3x3d hmat = currentSnap_->getHmat();
#	Line 1222 \| Line 1222 \| namespace OpenMD {
1222
1223		int selei;
1224		StuntDouble* sd;
1225	–	int idx;
1225
1226		vector<StuntDouble*> hotBin, coldBin;
1227
#	Line 1237 \| Line 1236 \| namespace OpenMD {
1236		for (sd = seleMan_.beginSelected(selei); sd != NULL;
1237		sd = seleMan_.nextSelected(selei)) {
1238
1240	–	idx = sd->getLocalIndex();
1241	–
1239		Vector3d pos = sd->getPos();
1240
1241		// wrap the stuntdouble's position back into the box:
#	Line 1391 \| Line 1388 \| namespace OpenMD {
1388		}
1389
1390		void RNEMD::doRNEMD() {
1391	<
1391	>	if (!doRNEMD_) return;
1392		trialCount_++;
1393		switch(rnemdMethod_) {
1394		case rnemdSwap:
#	Line 1410 \| Line 1407 \| namespace OpenMD {
1407		}
1408
1409		void RNEMD::collectData() {
1410	<
1410	>	if (!doRNEMD_) return;
1411		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1412		Mat3x3d hmat = currentSnap_->getHmat();
1413
1414	+	areaAccumulator_->add(currentSnap_->getXYarea());
1415	+
1416		seleMan_.setSelectionSet(evaluator_.evaluate());
1417
1418		int selei;
1419		StuntDouble* sd;
1421	–	int idx;
1420
1421		vector<RealType> binMass(nBins_, 0.0);
1422		vector<RealType> binPx(nBins_, 0.0);
#	Line 1443 \| Line 1441 \| namespace OpenMD {
1441		sd = mol->nextIntegrableObject(iiter))
1442		*/
1443		for (sd = seleMan_.beginSelected(selei); sd != NULL;
1444	<	sd = seleMan_.nextSelected(selei)) {
1444	>	sd = seleMan_.nextSelected(selei)) {
1445
1448	–	idx = sd->getLocalIndex();
1449	–
1446		Vector3d pos = sd->getPos();
1447
1448		// wrap the stuntdouble's position back into the box:
#	Line 1454 \| Line 1450 \| namespace OpenMD {
1450		if (usePeriodicBoundaryConditions_)
1451		currentSnap_->wrapVector(pos);
1452
1453	+
1454		// which bin is this stuntdouble in?
1455		// wrapped positions are in the range [-0.5hmat(2,2), +0.5hmat(2,2)]
1456		// Shift molecules by half a box to have bins start at 0
#	Line 1518 \| Line 1515 \| namespace OpenMD {
1515		vel.x() = binPx[i] / binMass[i];
1516		vel.y() = binPy[i] / binMass[i];
1517		vel.z() = binPz[i] / binMass[i];
1518	<	den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1518	>
1519	>	den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1520	>	/ currentSnap_->getVolume() ;
1521	>
1522		temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1523		PhysicalConstants::energyConvert);
1524
#	Line 1526 \| Line 1526 \| namespace OpenMD {
1526		if(outputMask_[j]) {
1527		switch(j) {
1528		case Z:
1529	<	(data_[j].accumulator[i])->add(z);
1529	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1530		break;
1531		case TEMPERATURE:
1532	<	data_[j].accumulator[i]->add(temp);
1532	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1533		break;
1534		case VELOCITY:
1535		dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1536		break;
1537		case DENSITY:
1538	<	data_[j].accumulator[i]->add(den);
1538	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1539		break;
1540		}
1541		}
#	Line 1544 \| Line 1544 \| namespace OpenMD {
1544		}
1545
1546		void RNEMD::getStarted() {
1547	+	if (!doRNEMD_) return;
1548		collectData();
1549		writeOutputFile();
1550		}
1551
1552		void RNEMD::parseOutputFileFormat(const std::string& format) {
1553	+	if (!doRNEMD_) return;
1554		StringTokenizer tokenizer(format, " ,;\|\t\n\r");
1555
1556		while(tokenizer.hasMoreTokens()) {
#	Line 1569 \| Line 1571 \| namespace OpenMD {
1571		}
1572
1573		void RNEMD::writeOutputFile() {
1574	+	if (!doRNEMD_) return;
1575
1576		#ifdef IS_MPI
1577		// If we're the root node, should we print out the results
#	Line 1588 \| Line 1591 \| namespace OpenMD {
1591		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1592
1593		RealType time = currentSnap_->getTime();
1594	<
1595	<
1594	>	RealType avgArea;
1595	>	areaAccumulator_->getAverage(avgArea);
1596	>	RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1597	>	/ PhysicalConstants::energyConvert;
1598	>	Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);
1599	>
1600		rnemdFile_ << "#######################################################\n";
1601		rnemdFile_ << "# RNEMD {\n";
1602
1603		map<string, RNEMDMethod>::iterator mi;
1604		for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1605		if ( (*mi).second == rnemdMethod_)
1606	<	rnemdFile_ << "# exchangeMethod = " << (*mi).first << "\n";
1606	>	rnemdFile_ << "# exchangeMethod = \"" << (*mi).first << "\";\n";
1607		}
1608		map<string, RNEMDFluxType>::iterator fi;
1609		for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1610		if ( (*fi).second == rnemdFluxType_)
1611	<	rnemdFile_ << "# fluxType = " << (*fi).first << "\n";
1611	>	rnemdFile_ << "# fluxType = \"" << (*fi).first << "\";\n";
1612		}
1613
1614	<	rnemdFile_ << "# exchangeTime = " << exchangeTime_ << " fs\n";
1614	>	rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n";
1615
1616		rnemdFile_ << "# objectSelection = \""
1617	<	<< rnemdObjectSelection_ << "\"\n";
1618	<	rnemdFile_ << "# slabWidth = " << slabWidth_ << " angstroms\n";
1619	<	rnemdFile_ << "# slabAcenter = " << slabACenter_ << " angstroms\n";
1620	<	rnemdFile_ << "# slabBcenter = " << slabBCenter_ << " angstroms\n";
1617	>	<< rnemdObjectSelection_ << "\";\n";
1618	>	rnemdFile_ << "# slabWidth = " << slabWidth_ << ";\n";
1619	>	rnemdFile_ << "# slabAcenter = " << slabACenter_ << ";\n";
1620	>	rnemdFile_ << "# slabBcenter = " << slabBCenter_ << ";\n";
1621		rnemdFile_ << "# }\n";
1622		rnemdFile_ << "#######################################################\n";
1623	<
1624	<	rnemdFile_ << "# running time = " << time << " fs\n";
1625	<	rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1626	<	rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1627	<
1628	<	rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1629	<	<< "\n";
1630	<	rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1631	<	<< "\n";
1632	<
1633	<	rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1634	<	rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1635	<	<< "\n";
1636	<
1637	<	rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1638	<	rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1639	<
1640	<
1623	>	rnemdFile_ << "# RNEMD report:\n";
1624	>	rnemdFile_ << "# running time = " << time << " fs\n";
1625	>	rnemdFile_ << "# target flux:\n";
1626	>	rnemdFile_ << "# kinetic = "
1627	>	<< kineticFlux_ / PhysicalConstants::energyConvert
1628	>	<< " (kcal/mol/A^2/fs)\n";
1629	>	rnemdFile_ << "# momentum = " << momentumFluxVector_
1630	>	<< " (amu/A/fs^2)\n";
1631	>	rnemdFile_ << "# target one-time exchanges:\n";
1632	>	rnemdFile_ << "# kinetic = "
1633	>	<< kineticTarget_ / PhysicalConstants::energyConvert
1634	>	<< " (kcal/mol)\n";
1635	>	rnemdFile_ << "# momentum = " << momentumTarget_
1636	>	<< " (amu*A/fs)\n";
1637	>	rnemdFile_ << "# actual exchange totals:\n";
1638	>	rnemdFile_ << "# kinetic = "
1639	>	<< kineticExchange_ / PhysicalConstants::energyConvert
1640	>	<< " (kcal/mol)\n";
1641	>	rnemdFile_ << "# momentum = " << momentumExchange_
1642	>	<< " (amu*A/fs)\n";
1643	>	rnemdFile_ << "# actual flux:\n";
1644	>	rnemdFile_ << "# kinetic = " << Jz
1645	>	<< " (kcal/mol/A^2/fs)\n";
1646	>	rnemdFile_ << "# momentum = " << JzP
1647	>	<< " (amu/A/fs^2)\n";
1648	>	rnemdFile_ << "# exchange statistics:\n";
1649	>	rnemdFile_ << "# attempted = " << trialCount_ << "\n";
1650	>	rnemdFile_ << "# failed = " << failTrialCount_ << "\n";
1651		if (rnemdMethod_ == rnemdNIVS) {
1652	<	rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1652	>	rnemdFile_ << "# NIVS root-check errors = "
1653	>	<< failRootCount_ << "\n";
1654		}
1637	–
1655		rnemdFile_ << "#######################################################\n";
1656
1657
#	Line 1645 \| Line 1662 \| namespace OpenMD {
1662		if (outputMask_[i]) {
1663		rnemdFile_ << "\t" << data_[i].title <<
1664		"(" << data_[i].units << ")";
1665	+	// add some extra tabs for column alignment
1666	+	if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1667		}
1668		}
1669		rnemdFile_ << std::endl;
1670
1671		rnemdFile_.precision(8);
1672
1673	<	for (unsigned int j = 0; j < nBins_; j++) {
1673	>	for (int j = 0; j < nBins_; j++) {
1674
1675		for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1676		if (outputMask_[i]) {
#	Line 1671 \| Line 1690 \| namespace OpenMD {
1690		rnemdFile_ << std::endl;
1691
1692		}
1693	+
1694	+	rnemdFile_ << "#######################################################\n";
1695	+	rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1696	+	rnemdFile_ << "#######################################################\n";
1697	+
1698	+
1699	+	for (int j = 0; j < nBins_; j++) {
1700	+	rnemdFile_ << "#";
1701	+	for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1702	+	if (outputMask_[i]) {
1703	+	if (data_[i].dataType == "RealType")
1704	+	writeRealStdDev(i,j);
1705	+	else if (data_[i].dataType == "Vector3d")
1706	+	writeVectorStdDev(i,j);
1707	+	else {
1708	+	sprintf( painCave.errMsg,
1709	+	"RNEMD found an unknown data type for: %s ",
1710	+	data_[i].title.c_str());
1711	+	painCave.isFatal = 1;
1712	+	simError();
1713	+	}
1714	+	}
1715	+	}
1716	+	rnemdFile_ << std::endl;
1717	+
1718	+	}
1719
1720		rnemdFile_.flush();
1721		rnemdFile_.close();
#	Line 1682 \| Line 1727 \| namespace OpenMD {
1727		}
1728
1729		void RNEMD::writeReal(int index, unsigned int bin) {
1730	+	if (!doRNEMD_) return;
1731		assert(index >=0 && index < ENDINDEX);
1732	<	assert(bin >=0 && bin < nBins_);
1732	>	assert(bin < nBins_);
1733		RealType s;
1734
1735	<	data_[index].accumulator[bin]->getAverage(s);
1735	>	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
1736
1737		if (! isinf(s) && ! isnan(s)) {
1738		rnemdFile_ << "\t" << s;
#	Line 1700 \| Line 1746 \| namespace OpenMD {
1746		}
1747
1748		void RNEMD::writeVector(int index, unsigned int bin) {
1749	+	if (!doRNEMD_) return;
1750		assert(index >=0 && index < ENDINDEX);
1751	<	assert(bin >=0 && bin < nBins_);
1751	>	assert(bin < nBins_);
1752		Vector3d s;
1753		dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1754		if (isinf(s[0]) \|\| isnan(s[0]) \|\|
#	Line 1716 \| Line 1763 \| namespace OpenMD {
1763		rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1764		}
1765		}
1766	+
1767	+	void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1768	+	if (!doRNEMD_) return;
1769	+	assert(index >=0 && index < ENDINDEX);
1770	+	assert(bin < nBins_);
1771	+	RealType s;
1772	+
1773	+	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
1774	+
1775	+	if (! isinf(s) && ! isnan(s)) {
1776	+	rnemdFile_ << "\t" << s;
1777	+	} else{
1778	+	sprintf( painCave.errMsg,
1779	+	"RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1780	+	data_[index].title.c_str(), bin);
1781	+	painCave.isFatal = 1;
1782	+	simError();
1783	+	}
1784	+	}
1785	+
1786	+	void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1787	+	if (!doRNEMD_) return;
1788	+	assert(index >=0 && index < ENDINDEX);
1789	+	assert(bin < nBins_);
1790	+	Vector3d s;
1791	+	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1792	+	if (isinf(s[0]) \|\| isnan(s[0]) \|\|
1793	+	isinf(s[1]) \|\| isnan(s[1]) \|\|
1794	+	isinf(s[2]) \|\| isnan(s[2]) ) {
1795	+	sprintf( painCave.errMsg,
1796	+	"RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1797	+	data_[index].title.c_str(), bin);
1798	+	painCave.isFatal = 1;
1799	+	simError();
1800	+	} else {
1801	+	rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1802	+	}
1803	+	}
1804		}
1805

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 1793 by gezelter, Fri Aug 31 21:16:10 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 1793 by gezelter, Fri Aug 31 21:16:10 2012 UTC