[ViewVC] Diff of: OpenMD/branches/development/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.
Revision 1850 by gezelter, Wed Feb 20 15:39:39 2013 UTC

#	Line 35 \| Line 35
35		*
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).
38	>	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39		* [4] Vardeman & Gezelter, in progress (2009).
40		*/
41
#	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];
1521	–	den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1522	–	temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1523	–	PhysicalConstants::energyConvert);
1509
1510	<	for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1511	<	if(outputMask_[j]) {
1512	<	switch(j) {
1513	<	case Z:
1514	<	(data_[j].accumulator[i])->add(z);
1515	<	break;
1516	<	case TEMPERATURE:
1517	<	data_[j].accumulator[i]->add(temp);
1518	<	break;
1519	<	case VELOCITY:
1520	<	dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1521	<	break;
1522	<	case DENSITY:
1523	<	data_[j].accumulator[i]->add(den);
1524	<	break;
1510	>	den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1511	>	/ currentSnap_->getVolume() ;
1512	>
1513	>	if (binCount[i] > 0) {
1514	>	// only add values if there are things to add
1515	>	temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1516	>	PhysicalConstants::energyConvert);
1517	>
1518	>	for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1519	>	if(outputMask_[j]) {
1520	>	switch(j) {
1521	>	case Z:
1522	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1523	>	break;
1524	>	case TEMPERATURE:
1525	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1526	>	break;
1527	>	case VELOCITY:
1528	>	dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1529	>	break;
1530	>	case DENSITY:
1531	>	dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1532	>	break;
1533	>	}
1534		}
1535		}
1536		}
#	Line 1544 \| Line 1538 \| namespace OpenMD {
1538		}
1539
1540		void RNEMD::getStarted() {
1541	+	if (!doRNEMD_) return;
1542		collectData();
1543		writeOutputFile();
1544		}
1545
1546		void RNEMD::parseOutputFileFormat(const std::string& format) {
1547	+	if (!doRNEMD_) return;
1548		StringTokenizer tokenizer(format, " ,;\|\t\n\r");
1549
1550		while(tokenizer.hasMoreTokens()) {
#	Line 1569 \| Line 1565 \| namespace OpenMD {
1565		}
1566
1567		void RNEMD::writeOutputFile() {
1568	+	if (!doRNEMD_) return;
1569
1570		#ifdef IS_MPI
1571		// If we're the root node, should we print out the results
#	Line 1588 \| Line 1585 \| namespace OpenMD {
1585		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1586
1587		RealType time = currentSnap_->getTime();
1588	<
1589	<
1588	>	RealType avgArea;
1589	>	areaAccumulator_->getAverage(avgArea);
1590	>	RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1591	>	/ PhysicalConstants::energyConvert;
1592	>	Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);
1593	>
1594		rnemdFile_ << "#######################################################\n";
1595		rnemdFile_ << "# RNEMD {\n";
1596
1597		map<string, RNEMDMethod>::iterator mi;
1598		for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1599		if ( (*mi).second == rnemdMethod_)
1600	<	rnemdFile_ << "# exchangeMethod = " << (*mi).first << "\n";
1600	>	rnemdFile_ << "# exchangeMethod = \"" << (*mi).first << "\";\n";
1601		}
1602		map<string, RNEMDFluxType>::iterator fi;
1603		for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1604		if ( (*fi).second == rnemdFluxType_)
1605	<	rnemdFile_ << "# fluxType = " << (*fi).first << "\n";
1605	>	rnemdFile_ << "# fluxType = \"" << (*fi).first << "\";\n";
1606		}
1607
1608	<	rnemdFile_ << "# exchangeTime = " << exchangeTime_ << " fs\n";
1608	>	rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n";
1609
1610		rnemdFile_ << "# objectSelection = \""
1611	<	<< rnemdObjectSelection_ << "\"\n";
1612	<	rnemdFile_ << "# slabWidth = " << slabWidth_ << " angstroms\n";
1613	<	rnemdFile_ << "# slabAcenter = " << slabACenter_ << " angstroms\n";
1614	<	rnemdFile_ << "# slabBcenter = " << slabBCenter_ << " angstroms\n";
1611	>	<< rnemdObjectSelection_ << "\";\n";
1612	>	rnemdFile_ << "# slabWidth = " << slabWidth_ << ";\n";
1613	>	rnemdFile_ << "# slabAcenter = " << slabACenter_ << ";\n";
1614	>	rnemdFile_ << "# slabBcenter = " << slabBCenter_ << ";\n";
1615		rnemdFile_ << "# }\n";
1616		rnemdFile_ << "#######################################################\n";
1617	<
1618	<	rnemdFile_ << "# running time = " << time << " fs\n";
1619	<	rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1620	<	rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1621	<
1622	<	rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1623	<	<< "\n";
1624	<	rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1625	<	<< "\n";
1626	<
1627	<	rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1628	<	rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1629	<	<< "\n";
1630	<
1631	<	rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1632	<	rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1633	<
1634	<
1617	>	rnemdFile_ << "# RNEMD report:\n";
1618	>	rnemdFile_ << "# running time = " << time << " fs\n";
1619	>	rnemdFile_ << "# target flux:\n";
1620	>	rnemdFile_ << "# kinetic = "
1621	>	<< kineticFlux_ / PhysicalConstants::energyConvert
1622	>	<< " (kcal/mol/A^2/fs)\n";
1623	>	rnemdFile_ << "# momentum = " << momentumFluxVector_
1624	>	<< " (amu/A/fs^2)\n";
1625	>	rnemdFile_ << "# target one-time exchanges:\n";
1626	>	rnemdFile_ << "# kinetic = "
1627	>	<< kineticTarget_ / PhysicalConstants::energyConvert
1628	>	<< " (kcal/mol)\n";
1629	>	rnemdFile_ << "# momentum = " << momentumTarget_
1630	>	<< " (amu*A/fs)\n";
1631	>	rnemdFile_ << "# actual exchange totals:\n";
1632	>	rnemdFile_ << "# kinetic = "
1633	>	<< kineticExchange_ / PhysicalConstants::energyConvert
1634	>	<< " (kcal/mol)\n";
1635	>	rnemdFile_ << "# momentum = " << momentumExchange_
1636	>	<< " (amu*A/fs)\n";
1637	>	rnemdFile_ << "# actual flux:\n";
1638	>	rnemdFile_ << "# kinetic = " << Jz
1639	>	<< " (kcal/mol/A^2/fs)\n";
1640	>	rnemdFile_ << "# momentum = " << JzP
1641	>	<< " (amu/A/fs^2)\n";
1642	>	rnemdFile_ << "# exchange statistics:\n";
1643	>	rnemdFile_ << "# attempted = " << trialCount_ << "\n";
1644	>	rnemdFile_ << "# failed = " << failTrialCount_ << "\n";
1645		if (rnemdMethod_ == rnemdNIVS) {
1646	<	rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1646	>	rnemdFile_ << "# NIVS root-check errors = "
1647	>	<< failRootCount_ << "\n";
1648		}
1637	–
1649		rnemdFile_ << "#######################################################\n";
1650
1651
#	Line 1645 \| Line 1656 \| namespace OpenMD {
1656		if (outputMask_[i]) {
1657		rnemdFile_ << "\t" << data_[i].title <<
1658		"(" << data_[i].units << ")";
1659	+	// add some extra tabs for column alignment
1660	+	if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1661		}
1662		}
1663		rnemdFile_ << std::endl;
1664
1665		rnemdFile_.precision(8);
1666
1667	<	for (unsigned int j = 0; j < nBins_; j++) {
1667	>	for (int j = 0; j < nBins_; j++) {
1668
1669		for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1670		if (outputMask_[i]) {
#	Line 1671 \| Line 1684 \| namespace OpenMD {
1684		rnemdFile_ << std::endl;
1685
1686		}
1687	+
1688	+	rnemdFile_ << "#######################################################\n";
1689	+	rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1690	+	rnemdFile_ << "#######################################################\n";
1691	+
1692	+
1693	+	for (int j = 0; j < nBins_; j++) {
1694	+	rnemdFile_ << "#";
1695	+	for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1696	+	if (outputMask_[i]) {
1697	+	if (data_[i].dataType == "RealType")
1698	+	writeRealStdDev(i,j);
1699	+	else if (data_[i].dataType == "Vector3d")
1700	+	writeVectorStdDev(i,j);
1701	+	else {
1702	+	sprintf( painCave.errMsg,
1703	+	"RNEMD found an unknown data type for: %s ",
1704	+	data_[i].title.c_str());
1705	+	painCave.isFatal = 1;
1706	+	simError();
1707	+	}
1708	+	}
1709	+	}
1710	+	rnemdFile_ << std::endl;
1711	+
1712	+	}
1713
1714		rnemdFile_.flush();
1715		rnemdFile_.close();
#	Line 1682 \| Line 1721 \| namespace OpenMD {
1721		}
1722
1723		void RNEMD::writeReal(int index, unsigned int bin) {
1724	+	if (!doRNEMD_) return;
1725		assert(index >=0 && index < ENDINDEX);
1726	<	assert(bin >=0 && bin < nBins_);
1726	>	assert(int(bin) < nBins_);
1727		RealType s;
1728
1729	<	data_[index].accumulator[bin]->getAverage(s);
1729	>	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
1730
1731		if (! isinf(s) && ! isnan(s)) {
1732		rnemdFile_ << "\t" << s;
#	Line 1700 \| Line 1740 \| namespace OpenMD {
1740		}
1741
1742		void RNEMD::writeVector(int index, unsigned int bin) {
1743	+	if (!doRNEMD_) return;
1744		assert(index >=0 && index < ENDINDEX);
1745	<	assert(bin >=0 && bin < nBins_);
1745	>	assert(int(bin) < nBins_);
1746		Vector3d s;
1747		dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1748		if (isinf(s[0]) \|\| isnan(s[0]) \|\|
#	Line 1716 \| Line 1757 \| namespace OpenMD {
1757		rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1758		}
1759		}
1760	+
1761	+	void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1762	+	if (!doRNEMD_) return;
1763	+	assert(index >=0 && index < ENDINDEX);
1764	+	assert(int(bin) < nBins_);
1765	+	RealType s;
1766	+
1767	+	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
1768	+
1769	+	if (! isinf(s) && ! isnan(s)) {
1770	+	rnemdFile_ << "\t" << s;
1771	+	} else{
1772	+	sprintf( painCave.errMsg,
1773	+	"RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1774	+	data_[index].title.c_str(), bin);
1775	+	painCave.isFatal = 1;
1776	+	simError();
1777	+	}
1778	+	}
1779	+
1780	+	void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1781	+	if (!doRNEMD_) return;
1782	+	assert(index >=0 && index < ENDINDEX);
1783	+	assert(int(bin) < nBins_);
1784	+	Vector3d s;
1785	+	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1786	+	if (isinf(s[0]) \|\| isnan(s[0]) \|\|
1787	+	isinf(s[1]) \|\| isnan(s[1]) \|\|
1788	+	isinf(s[2]) \|\| isnan(s[2]) ) {
1789	+	sprintf( painCave.errMsg,
1790	+	"RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1791	+	data_[index].title.c_str(), bin);
1792	+	painCave.isFatal = 1;
1793	+	simError();
1794	+	} else {
1795	+	rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1796	+	}
1797	+	}
1798		}
1799

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. Revision 1850 by gezelter, Wed Feb 20 15:39:39 2013 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.
Revision 1850 by gezelter, Wed Feb 20 15:39:39 2013 UTC