ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/integrators/RNEMD.cpp

Comparing trunk/src/integrators/RNEMD.cpp (file contents):
Revision 1396 by gezelter, Sat Dec 5 02:57:05 2009 UTC vs.
Revision 1560 by skuang, Wed May 11 17:55:32 2011 UTC

#	Line 111 | Line 111 | namespace OpenMD {
111		std::cerr << "WARNING! RNEMD Type Unknown!\n";
112		}
113
114	+	output3DTemp_ = false;
115	+	if (simParams->haveRNEMD_outputDimensionalTemperature()) {
116	+	output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature();
117	+	}
118	+
119		#ifdef IS_MPI
120		if (worldRank == 0) {
121		#endif
122
123		std::string rnemdFileName;
119	–	std::string xTempFileName;
120	–	std::string yTempFileName;
121	–	std::string zTempFileName;
124		switch(rnemdType_) {
125		case rnemdKineticSwap :
126		case rnemdKineticScale :
#	Line 129 | Line 131 | namespace OpenMD {
131		case rnemdPy :
132		case rnemdPyScale :
133		rnemdFileName = "momemtum.log";
132	–	xTempFileName = "temperatureX.log";
133	–	yTempFileName = "temperatureY.log";
134	–	zTempFileName = "temperatureZ.log";
135	–	xTempLog_.open(xTempFileName.c_str());
136	–	yTempLog_.open(yTempFileName.c_str());
137	–	zTempLog_.open(zTempFileName.c_str());
134		break;
135		case rnemdPz :
136		case rnemdPzScale :
#	Line 145 | Line 141 | namespace OpenMD {
141		}
142		rnemdLog_.open(rnemdFileName.c_str());
143
144	+	std::string xTempFileName;
145	+	std::string yTempFileName;
146	+	std::string zTempFileName;
147	+	if (output3DTemp_) {
148	+	xTempFileName = "temperatureX.log";
149	+	yTempFileName = "temperatureY.log";
150	+	zTempFileName = "temperatureZ.log";
151	+	xTempLog_.open(xTempFileName.c_str());
152	+	yTempLog_.open(yTempFileName.c_str());
153	+	zTempLog_.open(zTempFileName.c_str());
154	+	}
155	+
156		#ifdef IS_MPI
157		}
158		#endif
#	Line 152 | Line 160 | namespace OpenMD {
160		set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
161		set_RNEMD_nBins(simParams->getRNEMD_nBins());
162		midBin_ = nBins_ / 2;
163	+	if (simParams->haveRNEMD_binShift()) {
164	+	if (simParams->getRNEMD_binShift()) {
165	+	zShift_ = 0.5 / (RealType)(nBins_);
166	+	} else {
167	+	zShift_ = 0.0;
168	+	}
169	+	} else {
170	+	zShift_ = 0.0;
171	+	}
172	+	//std::cerr << "we have zShift_ = " << zShift_ << "\n";
173	+	//shift slabs by half slab width, might be useful in heterogeneous systems
174	+	//set to 0.0 if not using it; can NOT be used in status output yet
175		if (simParams->haveRNEMD_logWidth()) {
176	<	rnemdLogWidth_ = simParams->getRNEMD_logWidth();
176	>	set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
177	>	/*arbitary rnemdLogWidth_ no checking
178		if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
179		std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
180		std::cerr << "Automaically set back to default.\n";
181		rnemdLogWidth_ = nBins_;
182	<	}
182	>	}*/
183		} else {
184	<	rnemdLogWidth_ = nBins_;
184	>	set_RNEMD_logWidth(nBins_);
185		}
186		valueHist_.resize(rnemdLogWidth_, 0.0);
187		valueCount_.resize(rnemdLogWidth_, 0);
188		xTempHist_.resize(rnemdLogWidth_, 0.0);
189		yTempHist_.resize(rnemdLogWidth_, 0.0);
190		zTempHist_.resize(rnemdLogWidth_, 0.0);
191	+	xyzTempCount_.resize(rnemdLogWidth_, 0);
192
193		set_RNEMD_exchange_total(0.0);
194		if (simParams->haveRNEMD_targetFlux()) {
#	Line 202 | Line 224 | namespace OpenMD {
224		rnemdLog_.close();
225		if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale)
226		std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n";
227	<	if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) {
227	>	if (output3DTemp_) {
228		xTempLog_.close();
229		yTempLog_.close();
230		zTempLog_.close();
#	Line 246 | Line 268 | namespace OpenMD {
268		// which bin is this stuntdouble in?
269		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
270
271	<	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
271	>	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
272
273
274		// if we're in bin 0 or the middleBin
#	Line 261 | Line 283 | namespace OpenMD {
283
284		value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
285		vel[2]*vel[2]);
286	+	/*
287		if (sd->isDirectional()) {
288		Vector3d angMom = sd->getJ();
289		Mat3x3d I = sd->getI();
#	Line 276 | Line 299 | namespace OpenMD {
299		+ angMom[1]*angMom[1]/I(1, 1)
300		+ angMom[2]*angMom[2]/I(2, 2);
301		}
302	<	}
302	>	} no exchange of angular momenta
303	>	*/
304		//make exchangeSum_ comparable between swap & scale
305		//temporarily without using energyConvert
306		//value = value * 0.5 / PhysicalConstants::energyConvert;
#	Line 330 | Line 354 | namespace OpenMD {
354		bool my_min_found = min_found;
355		bool my_max_found = max_found;
356
357	<	// Even if we didn't find a minimum, did someone else?
358	<	MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
359	<	1, MPI::BOOL, MPI::LAND);
360	<
357	>	// Even if we didn't find a minimum, did someone else? debugging...
358	>	//MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
359	>	//                          1, MPI::BOOL, MPI::LAND);
360	>	MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
361		// Even if we didn't find a maximum, did someone else?
362	<	MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
363	<	1, MPI::BOOL, MPI::LAND);
364	<
362	>	//MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
363	>	//                          1, MPI::BOOL, MPI::LAND);
364	>	MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
365		struct {
366		RealType val;
367		int rank;
#	Line 373 | Line 397 | namespace OpenMD {
397		#endif
398
399		if (max_found && min_found) {
400	<	if (min_val< max_val) {
400	>	if (min_val < max_val) {
401
402		#ifdef IS_MPI
403		if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
404		// I have both maximum and minimum, so proceed like a single
405		// processor version:
406		#endif
407	<	// objects to be swapped: velocity & angular velocity
407	>	// objects to be swapped: velocity ONLY
408		Vector3d min_vel = min_sd->getVel();
409		Vector3d max_vel = max_sd->getVel();
410		RealType temp_vel;
#	Line 389 | Line 413 | namespace OpenMD {
413		case rnemdKineticSwap :
414		min_sd->setVel(max_vel);
415		max_sd->setVel(min_vel);
416	+	/*
417		if (min_sd->isDirectional() && max_sd->isDirectional()) {
418		Vector3d min_angMom = min_sd->getJ();
419		Vector3d max_angMom = max_sd->getJ();
420		min_sd->setJ(max_angMom);
421		max_sd->setJ(min_angMom);
422	<	}
422	>	} no angular momentum exchange
423	>	*/
424		break;
425		case rnemdPx :
426		temp_vel = min_vel.x();
#	Line 438 | Line 464 | namespace OpenMD {
464		switch(rnemdType_) {
465		case rnemdKineticSwap :
466		max_sd->setVel(min_vel);
467	<
467	>	/*
468		if (max_sd->isDirectional()) {
469		Vector3d min_angMom;
470		Vector3d max_angMom = max_sd->getJ();
#	Line 451 | Line 477 | namespace OpenMD {
477		status);
478
479		max_sd->setJ(min_angMom);
480	<	}
480	>	} no angular momentum exchange
481	>	*/
482		break;
483		case rnemdPx :
484		max_vel.x() = min_vel.x();
#	Line 484 | Line 511 | namespace OpenMD {
511		switch(rnemdType_) {
512		case rnemdKineticSwap :
513		min_sd->setVel(max_vel);
514	<
514	>	/*
515		if (min_sd->isDirectional()) {
516		Vector3d min_angMom = min_sd->getJ();
517		Vector3d max_angMom;
#	Line 497 | Line 524 | namespace OpenMD {
524		status);
525
526		min_sd->setJ(max_angMom);
527	<	}
527	>	} no angular momentum exchange
528	>	*/
529		break;
530		case rnemdPx :
531		min_vel.x() = max_vel.x();
#	Line 570 | Line 598 | namespace OpenMD {
598		// which bin is this stuntdouble in?
599		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
600
601	<	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
601	>	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
602
603		// if we're in bin 0 or the middleBin
604		if (binNo == 0 || binNo == midBin_) {
#	Line 776 | Line 804 | namespace OpenMD {
804		+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
805		break;
806		case rnemdPzScale :
807	+	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
808	+	+ fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
809		default :
810		break;
811		}
#	Line 879 | Line 909 | namespace OpenMD {
909		StuntDouble* sd;
910		int idx;
911
912	+	// alternative approach, track all molecules instead of only those selected for scaling/swapping
913	+	//SimInfo::MoleculeIterator miter;
914	+	//std::vector<StuntDouble*>::iterator iiter;
915	+	//Molecule* mol;
916	+	//StuntDouble* integrableObject;
917	+	//for (mol = info_->beginMolecule(miter); mol != NULL;
918	+	//     mol = info_->nextMolecule(miter))
919	+	// integrableObject is essentially sd
920	+	//for (integrableObject = mol->beginIntegrableObject(iiter);
921	+	//     integrableObject != NULL;
922	+	//     integrableObject = mol->nextIntegrableObject(iiter))
923		for (sd = seleMan_.beginSelected(selei); sd != NULL;
924		sd = seleMan_.nextSelected(selei)) {
925
#	Line 894 | Line 935 | namespace OpenMD {
935		// which bin is this stuntdouble in?
936		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
937
938	<	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
939	<
938	>	int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
939	>	rnemdLogWidth_;
940	>	/* no symmetrization allowed due to arbitary rnemdLogWidth_ value
941		if (rnemdLogWidth_ == midBin_ + 1)
942		if (binNo > midBin_)
943		binNo = nBins_ - binNo;
944	<
944	>	*/
945		RealType mass = sd->getMass();
946		Vector3d vel = sd->getVel();
947		RealType value;
#	Line 940 | Line 982 | namespace OpenMD {
982		case rnemdPxScale :
983		value = mass * vel[0];
984		valueCount_[binNo]++;
943	–	xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
944	–	/ PhysicalConstants::kb;
945	–	yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
946	–	/ PhysicalConstants::kb;
947	–	zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
948	–	/ PhysicalConstants::kb;
949	–	xTempHist_[binNo] += xVal;
950	–	yTempHist_[binNo] += yVal;
951	–	zTempHist_[binNo] += zVal;
985		break;
986		case rnemdPy :
987		case rnemdPyScale :
#	Line 957 | Line 990 | namespace OpenMD {
990		break;
991		case rnemdPz :
992		case rnemdPzScale :
993	<	value = mass * vel[2];
993	>	value = pos.z(); //temporarily for homogeneous systems ONLY
994		valueCount_[binNo]++;
995		break;
996		case rnemdUnknown :
997		default :
998	+	value = 1.0;
999	+	valueCount_[binNo]++;
1000		break;
1001		}
1002		valueHist_[binNo] += value;
968	–	}
1003
1004	+	if (output3DTemp_) {
1005	+	xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
1006	+	/ PhysicalConstants::kb;
1007	+	yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1008	+	/ PhysicalConstants::kb;
1009	+	zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1010	+	/ PhysicalConstants::kb;
1011	+	xTempHist_[binNo] += xVal;
1012	+	yTempHist_[binNo] += yVal;
1013	+	zTempHist_[binNo] += zVal;
1014	+	xyzTempCount_[binNo]++;
1015	+	}
1016	+	}
1017		}
1018
1019		void RNEMD::getStarted() {
1020	+	collectData();
1021	+	/* now should be able to output profile in step 0, but might not be useful
1022		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1023		Stats& stat = currentSnap_->statData;
1024		stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1025	+	*/
1026	+	getStatus();
1027		}
1028
1029		void RNEMD::getStatus() {
#	Line 994 | Line 1045 | namespace OpenMD {
1045		rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1046		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
1047		rnemdLogWidth_, MPI::INT, MPI::SUM);
1048	<	if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale) {
1048	>	if (output3DTemp_) {
1049		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
1050		rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1051		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1052		rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1053		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1054		rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1055	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0],
1056	+	rnemdLogWidth_, MPI::INT, MPI::SUM);
1057		}
1058		// If we're the root node, should we print out the results
1059		int worldRank = MPI::COMM_WORLD.Get_rank();
#	Line 1011 | Line 1064 | namespace OpenMD {
1064		rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
1065		}
1066		rnemdLog_ << "\n";
1067	<	if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) {
1067	>	if (output3DTemp_) {
1068		xTempLog_ << time;
1069		for (j = 0; j < rnemdLogWidth_; j++) {
1070	<	xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j];
1070	>	xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j];
1071		}
1072		xTempLog_ << "\n";
1073		yTempLog_ << time;
1074		for (j = 0; j < rnemdLogWidth_; j++) {
1075	<	yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j];
1075	>	yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1076		}
1077		yTempLog_ << "\n";
1078		zTempLog_ << time;
1079		for (j = 0; j < rnemdLogWidth_; j++) {
1080	<	zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j];
1080	>	zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1081		}
1082		zTempLog_ << "\n";
1083		}
#	Line 1035 | Line 1088 | namespace OpenMD {
1088		valueCount_[j] = 0;
1089		valueHist_[j] = 0.0;
1090		}
1091	<	if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale)
1091	>	if (output3DTemp_)
1092		for (j = 0; j < rnemdLogWidth_; j++) {
1093		xTempHist_[j] = 0.0;
1094		yTempHist_[j] = 0.0;
1095		zTempHist_[j] = 0.0;
1096	+	xyzTempCount_[j] = 0;
1097		}
1098		}
1099		}

Comparing trunk/src/integrators/RNEMD.cpp (property svn:keywords):
Revision 1396 by gezelter, Sat Dec 5 02:57:05 2009 UTC vs.
Revision 1560 by skuang, Wed May 11 17:55:32 2011 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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