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 1390 by gezelter, Wed Nov 25 20:02:06 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();
177	<	if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) {
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 194 | Line 216 | namespace OpenMD {
216
217		RNEMD::~RNEMD() {
218		delete randNumGen_;
219	<
198	<	std::cerr << "total fail trials: " << failTrialCount_ << "\n";
219	>
220		#ifdef IS_MPI
221		if (worldRank == 0) {
222		#endif
223	+	std::cerr << "total fail trials: " << failTrialCount_ << "\n";
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 642 | Line 670 | namespace OpenMD {
670		a001 = Kcx * px * px;
671		*/
672
673	<	//scale all three dimensions, let x = y
673	>	//scale all three dimensions, let c_x = c_y
674		a000 = Kcx + Kcy;
675		a110 = Kcz;
676		c0 = targetFlux_ - Kcx - Kcy - Kcz;
#	Line 678 | Line 706 | namespace OpenMD {
706		+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
707		break;
708		case rnemdPzScale ://we don't really do this, do we?
709	+	c = 1 - targetFlux_ / Pcz;
710	+	a000 = Kcx;
711	+	a110 = Kcy;
712	+	c0 = Kcz * c * c - Kcx - Kcy - Kcz;
713	+	a001 = px * px * Khx;
714	+	a111 = py * py * Khy;
715	+	b01 = -2.0 * Khx * px * (1.0 + px);
716	+	b11 = -2.0 * Khy * py * (1.0 + py);
717	+	c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
718	+	+ Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
719	+	break;
720		default :
721		break;
722		}
#	Line 721 | Line 760 | namespace OpenMD {
760		r2 = *ri;
761		//check if FindRealRoots() give the right answer
762		if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
763	<	std::cerr << "WARNING! eq solvers might have mistakes!\n";
763	>	sprintf(painCave.errMsg,
764	>	"RNEMD Warning: polynomial solve seems to have an error!");
765	>	painCave.isFatal = 0;
766	>	simError();
767		failRootCount_++;
768		}
769		//might not be useful w/o rescaling coefficients
#	Line 762 | 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 797 | Line 841 | namespace OpenMD {
841		z = bestPair.second;
842		break;
843		case rnemdPzScale :
844	+	x = bestPair.first;
845	+	y = bestPair.second;
846	+	z = c;
847	+	break;
848		default :
849		break;
850		}
#	Line 823 | Line 871 | namespace OpenMD {
871		exchangeSum_ += targetFlux_;
872		//we may want to check whether the exchange has been successful
873		} else {
874	<	std::cerr << "exchange NOT performed!\n";
874	>	std::cerr << "exchange NOT performed!\n";//MPI incompatible
875		failTrialCount_++;
876		}
877
#	Line 861 | 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 876 | 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 922 | Line 982 | namespace OpenMD {
982		case rnemdPxScale :
983		value = mass * vel[0];
984		valueCount_[binNo]++;
925	–	xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
926	–	/ PhysicalConstants::kb;
927	–	yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
928	–	/ PhysicalConstants::kb;
929	–	zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
930	–	/ PhysicalConstants::kb;
931	–	xTempHist_[binNo] += xVal;
932	–	yTempHist_[binNo] += yVal;
933	–	zTempHist_[binNo] += zVal;
985		break;
986		case rnemdPy :
987		case rnemdPyScale :
#	Line 939 | 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;
950	–	}
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 965 | Line 1034 | namespace OpenMD {
1034
1035		stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1036		//or to be more meaningful, define another item as exchangeSum_ / time
1037	+	int j;
1038
969	–
1039		#ifdef IS_MPI
1040
1041		// all processors have the same number of bins, and STL vectors pack their
#	Line 976 | 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	<
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();
1060		if (worldRank == 0) {
1061		#endif
984	–	int j;
1062		rnemdLog_ << time;
1063		for (j = 0; j < rnemdLogWidth_; j++) {
1064		rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
988	–	valueHist_[j] = 0.0;
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];
995	<	xTempHist_[j] = 0.0;
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];
1001	<	yTempHist_[j] = 0.0;
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];
1007	<	zTempHist_[j] = 0.0;
1080	>	zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1081		}
1082		zTempLog_ << "\n";
1083		}
1011	–	for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0;
1084		#ifdef IS_MPI
1085	<	}
1085	>	}
1086		#endif
1087	<
1088	<
1087	>	for (j = 0; j < rnemdLogWidth_; j++) {
1088	>	valueCount_[j] = 0;
1089	>	valueHist_[j] = 0.0;
1090	>	}
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		}
1018	–
1099		}

Comparing trunk/src/integrators/RNEMD.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 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