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 1561 by gezelter, Wed May 11 19:04:40 2011 UTC

#	Line 57 | Line 57
57
58		#define HONKING_LARGE_VALUE 1.0e10
59
60	+	using namespace std;
61		namespace OpenMD {
62
63	<	RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
63	>	RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
64	>	usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
65
66		failTrialCount_ = 0;
67		failRootCount_ = 0;
#	Line 88 | Line 90 | namespace OpenMD {
90
91		if (selectionCount > nIntegrable) {
92		sprintf(painCave.errMsg,
93	<	"RNEMD warning: The current RNEMD_objectSelection,\n"
93	>	"RNEMD: The current RNEMD_objectSelection,\n"
94		"\t\t%s\n"
95		"\thas resulted in %d selected objects.  However,\n"
96		"\tthe total number of integrable objects in the system\n"
#	Line 98 | Line 100 | namespace OpenMD {
100		rnemdObjectSelection_.c_str(),
101		selectionCount, nIntegrable);
102		painCave.isFatal = 0;
103	+	painCave.severity = OPENMD_WARNING;
104		simError();
102	–
105		}
106
107	<	const std::string st = simParams->getRNEMD_exchangeType();
107	>	const string st = simParams->getRNEMD_exchangeType();
108
109	<	std::map<std::string, RNEMDTypeEnum>::iterator i;
109	>	map<string, RNEMDTypeEnum>::iterator i;
110		i = stringToEnumMap_.find(st);
111		rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
112		if (rnemdType_ == rnemdUnknown) {
113	<	std::cerr << "WARNING! RNEMD Type Unknown!\n";
113	>	sprintf(painCave.errMsg,
114	>	"RNEMD: The current RNEMD_exchangeType,\n"
115	>	"\t\t%s\n"
116	>	"\tis not one of the recognized exchange types.\n",
117	>	st.c_str());
118	>	painCave.isFatal = 1;
119	>	painCave.severity = OPENMD_ERROR;
120	>	simError();
121		}
122	+
123	+	output3DTemp_ = false;
124	+	if (simParams->haveRNEMD_outputDimensionalTemperature()) {
125	+	output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature();
126	+	}
127
128		#ifdef IS_MPI
129		if (worldRank == 0) {
130		#endif
131
132	<	std::string rnemdFileName;
119	<	std::string xTempFileName;
120	<	std::string yTempFileName;
121	<	std::string zTempFileName;
132	>	string rnemdFileName;
133		switch(rnemdType_) {
134		case rnemdKineticSwap :
135		case rnemdKineticScale :
#	Line 129 | Line 140 | namespace OpenMD {
140		case rnemdPy :
141		case rnemdPyScale :
142		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());
143		break;
144		case rnemdPz :
145		case rnemdPzScale :
#	Line 145 | Line 150 | namespace OpenMD {
150		}
151		rnemdLog_.open(rnemdFileName.c_str());
152
153	+	string xTempFileName;
154	+	string yTempFileName;
155	+	string zTempFileName;
156	+	if (output3DTemp_) {
157	+	xTempFileName = "temperatureX.log";
158	+	yTempFileName = "temperatureY.log";
159	+	zTempFileName = "temperatureZ.log";
160	+	xTempLog_.open(xTempFileName.c_str());
161	+	yTempLog_.open(yTempFileName.c_str());
162	+	zTempLog_.open(zTempFileName.c_str());
163	+	}
164	+
165		#ifdef IS_MPI
166		}
167		#endif
#	Line 152 | Line 169 | namespace OpenMD {
169		set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
170		set_RNEMD_nBins(simParams->getRNEMD_nBins());
171		midBin_ = nBins_ / 2;
172	<	if (simParams->haveRNEMD_logWidth()) {
173	<	rnemdLogWidth_ = simParams->getRNEMD_logWidth();
174	<	if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) {
175	<	std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
176	<	std::cerr << "Automaically set back to default.\n";
160	<	rnemdLogWidth_ = nBins_;
172	>	if (simParams->haveRNEMD_binShift()) {
173	>	if (simParams->getRNEMD_binShift()) {
174	>	zShift_ = 0.5 / (RealType)(nBins_);
175	>	} else {
176	>	zShift_ = 0.0;
177		}
178		} else {
179	<	rnemdLogWidth_ = nBins_;
179	>	zShift_ = 0.0;
180		}
181	+	//cerr << "we have zShift_ = " << zShift_ << "\n";
182	+	//shift slabs by half slab width, might be useful in heterogeneous systems
183	+	//set to 0.0 if not using it; can NOT be used in status output yet
184	+	if (simParams->haveRNEMD_logWidth()) {
185	+	set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
186	+	/*arbitary rnemdLogWidth_ no checking
187	+	if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
188	+	cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
189	+	cerr << "Automaically set back to default.\n";
190	+	rnemdLogWidth_ = nBins_;
191	+	}*/
192	+	} else {
193	+	set_RNEMD_logWidth(nBins_);
194	+	}
195		valueHist_.resize(rnemdLogWidth_, 0.0);
196		valueCount_.resize(rnemdLogWidth_, 0);
197		xTempHist_.resize(rnemdLogWidth_, 0.0);
198		yTempHist_.resize(rnemdLogWidth_, 0.0);
199		zTempHist_.resize(rnemdLogWidth_, 0.0);
200	+	xyzTempCount_.resize(rnemdLogWidth_, 0);
201
202		set_RNEMD_exchange_total(0.0);
203		if (simParams->haveRNEMD_targetFlux()) {
#	Line 194 | Line 225 | namespace OpenMD {
225
226		RNEMD::~RNEMD() {
227		delete randNumGen_;
228	<
198	<	std::cerr << "total fail trials: " << failTrialCount_ << "\n";
228	>
229		#ifdef IS_MPI
230		if (worldRank == 0) {
231		#endif
232	+
233	+	sprintf(painCave.errMsg,
234	+	"RNEMD: total failed trials: %d\n",
235	+	failTrialCount_);
236	+	painCave.isFatal = 0;
237	+	painCave.severity = OPENMD_INFO;
238	+	simError();
239	+
240		rnemdLog_.close();
241	<	if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale)
242	<	std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n";
243	<	if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) {
241	>	if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) {
242	>	sprintf(painCave.errMsg,
243	>	"RNEMD: total root-checking warnings: %d\n",
244	>	failRootCount_);
245	>	painCave.isFatal = 0;
246	>	painCave.severity = OPENMD_INFO;
247	>	simError();
248	>	}
249	>	if (output3DTemp_) {
250		xTempLog_.close();
251		yTempLog_.close();
252		zTempLog_.close();
#	Line 246 | Line 290 | namespace OpenMD {
290		// which bin is this stuntdouble in?
291		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
292
293	<	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
293	>	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
294
295
296		// if we're in bin 0 or the middleBin
#	Line 261 | Line 305 | namespace OpenMD {
305
306		value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
307		vel[2]*vel[2]);
308	<	if (sd->isDirectional()) {
308	>	/*
309	>	if (sd->isDirectional()) {
310		Vector3d angMom = sd->getJ();
311		Mat3x3d I = sd->getI();
312
313		if (sd->isLinear()) {
314	<	int i = sd->linearAxis();
315	<	int j = (i + 1) % 3;
316	<	int k = (i + 2) % 3;
317	<	value += angMom[j] * angMom[j] / I(j, j) +
318	<	angMom[k] * angMom[k] / I(k, k);
314	>	int i = sd->linearAxis();
315	>	int j = (i + 1) % 3;
316	>	int k = (i + 2) % 3;
317	>	value += angMom[j] * angMom[j] / I(j, j) +
318	>	angMom[k] * angMom[k] / I(k, k);
319		} else {
320	<	value += angMom[0]*angMom[0]/I(0, 0)
321	<	+ angMom[1]*angMom[1]/I(1, 1)
322	<	+ angMom[2]*angMom[2]/I(2, 2);
320	>	value += angMom[0]*angMom[0]/I(0, 0)
321	>	+ angMom[1]*angMom[1]/I(1, 1)
322	>	+ angMom[2]*angMom[2]/I(2, 2);
323		}
324	<	}
324	>	} no exchange of angular momenta
325	>	*/
326		//make exchangeSum_ comparable between swap & scale
327		//temporarily without using energyConvert
328		//value = value * 0.5 / PhysicalConstants::energyConvert;
#	Line 331 | Line 377 | namespace OpenMD {
377		bool my_max_found = max_found;
378
379		// Even if we didn't find a minimum, did someone else?
380	<	MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found,
335	<	1, MPI::BOOL, MPI::LAND);
336	<
380	>	MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
381		// Even if we didn't find a maximum, did someone else?
382	<	MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found,
339	<	1, MPI::BOOL, MPI::LAND);
340	<
382	>	MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
383		struct {
384		RealType val;
385		int rank;
#	Line 373 | Line 415 | namespace OpenMD {
415		#endif
416
417		if (max_found && min_found) {
418	<	if (min_val< max_val) {
418	>	if (min_val < max_val) {
419
420		#ifdef IS_MPI
421		if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
422		// I have both maximum and minimum, so proceed like a single
423		// processor version:
424		#endif
425	<	// objects to be swapped: velocity & angular velocity
425	>	// objects to be swapped: velocity ONLY
426		Vector3d min_vel = min_sd->getVel();
427		Vector3d max_vel = max_sd->getVel();
428		RealType temp_vel;
#	Line 389 | Line 431 | namespace OpenMD {
431		case rnemdKineticSwap :
432		min_sd->setVel(max_vel);
433		max_sd->setVel(min_vel);
434	<	if (min_sd->isDirectional() && max_sd->isDirectional()) {
434	>	/*
435	>	if (min_sd->isDirectional() && max_sd->isDirectional()) {
436		Vector3d min_angMom = min_sd->getJ();
437		Vector3d max_angMom = max_sd->getJ();
438		min_sd->setJ(max_angMom);
439		max_sd->setJ(min_angMom);
440	<	}
440	>	} no angular momentum exchange
441	>	*/
442		break;
443		case rnemdPx :
444		temp_vel = min_vel.x();
#	Line 438 | Line 482 | namespace OpenMD {
482		switch(rnemdType_) {
483		case rnemdKineticSwap :
484		max_sd->setVel(min_vel);
485	<
485	>	//no angular momentum exchange for now
486	>	/*
487		if (max_sd->isDirectional()) {
488		Vector3d min_angMom;
489		Vector3d max_angMom = max_sd->getJ();
490	<
490	>
491		// point-to-point swap of the angular momentum vector
492		MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
493		MPI::REALTYPE, min_vals.rank, 1,
494		min_angMom.getArrayPointer(), 3,
495		MPI::REALTYPE, min_vals.rank, 1,
496		status);
497	<
497	>
498		max_sd->setJ(min_angMom);
499	<	}
499	>	}
500	>	*/
501		break;
502		case rnemdPx :
503		max_vel.x() = min_vel.x();
#	Line 484 | Line 530 | namespace OpenMD {
530		switch(rnemdType_) {
531		case rnemdKineticSwap :
532		min_sd->setVel(max_vel);
533	<
533	>	// no angular momentum exchange for now
534	>	/*
535		if (min_sd->isDirectional()) {
536		Vector3d min_angMom = min_sd->getJ();
537		Vector3d max_angMom;
538	<
538	>
539		// point-to-point swap of the angular momentum vector
540		MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
541		MPI::REALTYPE, max_vals.rank, 1,
542		max_angMom.getArrayPointer(), 3,
543		MPI::REALTYPE, max_vals.rank, 1,
544		status);
545	<
545	>
546		min_sd->setJ(max_angMom);
547		}
548	+	*/
549		break;
550		case rnemdPx :
551		min_vel.x() = max_vel.x();
#	Line 517 | Line 565 | namespace OpenMD {
565		}
566		#endif
567		exchangeSum_ += max_val - min_val;
568	<	} else {
569	<	std::cerr << "exchange NOT performed!\nmin_val > max_val.\n";
568	>	} else {
569	>	sprintf(painCave.errMsg,
570	>	"RNEMD: exchange NOT performed because min_val > max_val\n");
571	>	painCave.isFatal = 0;
572	>	painCave.severity = OPENMD_INFO;
573	>	simError();
574		failTrialCount_++;
575		}
576		} else {
577	<	std::cerr << "exchange NOT performed!\n";
578	<	std::cerr << "at least one of the two slabs empty.\n";
577	>	sprintf(painCave.errMsg,
578	>	"RNEMD: exchange NOT performed because at least one\n"
579	>	"\tof the two slabs is empty\n");
580	>	painCave.isFatal = 0;
581	>	painCave.severity = OPENMD_INFO;
582	>	simError();
583		failTrialCount_++;
584		}
585
#	Line 540 | Line 596 | namespace OpenMD {
596		StuntDouble* sd;
597		int idx;
598
599	<	std::vector<StuntDouble*> hotBin, coldBin;
599	>	vector<StuntDouble*> hotBin, coldBin;
600
601		RealType Phx = 0.0;
602		RealType Phy = 0.0;
#	Line 570 | Line 626 | namespace OpenMD {
626		// which bin is this stuntdouble in?
627		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
628
629	<	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
629	>	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
630
631		// if we're in bin 0 or the middleBin
632		if (binNo == 0 || binNo == midBin_) {
#	Line 629 | Line 685 | namespace OpenMD {
685		RealType a000, a110, c0, a001, a111, b01, b11, c1, c;
686		switch(rnemdType_) {
687		case rnemdKineticScale :
688	<	/*used hotBin coeff's & only scale x & y dimensions
688	>	// used hotBin coeff's & only scale x & y dimensions
689	>	/*
690		RealType px = Phx / Pcx;
691		RealType py = Phy / Pcy;
692		a110 = Khy;
693		c0 = - Khx - Khy - targetFlux_;
694		a000 = Khx;
695	<	a111 = Kcy * py * py
695	>	a111 = Kcy * py * py;
696		b11 = -2.0 * Kcy * py * (1.0 + py);
697		c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_;
698		b01 = -2.0 * Kcx * px * (1.0 + px);
699		a001 = Kcx * px * px;
700	<	*/
701	<
645	<	//scale all three dimensions, let x = y
700	>	*/
701	>	//scale all three dimensions, let c_x = c_y
702		a000 = Kcx + Kcy;
703		a110 = Kcz;
704		c0 = targetFlux_ - Kcx - Kcy - Kcz;
#	Line 651 | Line 707 | namespace OpenMD {
707		b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
708		b11 = -2.0 * Khz * pz * (1.0 + pz);
709		c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
710	<	+ Khz * pz * (2.0 + pz) - targetFlux_;
710	>	+ Khz * pz * (2.0 + pz) - targetFlux_;
711		break;
712		case rnemdPxScale :
713		c = 1 - targetFlux_ / Pcx;
#	Line 663 | Line 719 | namespace OpenMD {
719		b01 = -2.0 * Khy * py * (1.0 + py);
720		b11 = -2.0 * Khz * pz * (1.0 + pz);
721		c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
722	<	+ Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
722	>	+ Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
723		break;
724		case rnemdPyScale :
725		c = 1 - targetFlux_ / Pcy;
#	Line 675 | Line 731 | namespace OpenMD {
731		b01 = -2.0 * Khx * px * (1.0 + px);
732		b11 = -2.0 * Khz * pz * (1.0 + pz);
733		c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
734	<	+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
734	>	+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
735		break;
736		case rnemdPzScale ://we don't really do this, do we?
737	+	c = 1 - targetFlux_ / Pcz;
738	+	a000 = Kcx;
739	+	a110 = Kcy;
740	+	c0 = Kcz * c * c - Kcx - Kcy - Kcz;
741	+	a001 = px * px * Khx;
742	+	a111 = py * py * Khy;
743	+	b01 = -2.0 * Khx * px * (1.0 + px);
744	+	b11 = -2.0 * Khy * py * (1.0 + py);
745	+	c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
746	+	+ Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
747	+	break;
748		default :
749		break;
750		}
#	Line 712 | Line 779 | namespace OpenMD {
779		poly.setCoefficient(2, u2);
780		poly.setCoefficient(1, u1);
781		poly.setCoefficient(0, u0);
782	<	std::vector<RealType> realRoots = poly.FindRealRoots();
782	>	vector<RealType> realRoots = poly.FindRealRoots();
783
784	<	std::vector<RealType>::iterator ri;
784	>	vector<RealType>::iterator ri;
785		RealType r1, r2, alpha0;
786	<	std::vector<std::pair<RealType,RealType> > rps;
786	>	vector<pair<RealType,RealType> > rps;
787		for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
788		r2 = *ri;
789		//check if FindRealRoots() give the right answer
790		if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
791	<	std::cerr << "WARNING! eq solvers might have mistakes!\n";
791	>	sprintf(painCave.errMsg,
792	>	"RNEMD Warning: polynomial solve seems to have an error!");
793	>	painCave.isFatal = 0;
794	>	simError();
795		failRootCount_++;
796		}
797		//might not be useful w/o rescaling coefficients
#	Line 729 | Line 799 | namespace OpenMD {
799		if (alpha0 >= 0.0) {
800		r1 = sqrt(alpha0 / a000);
801		if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6)
802	<	{ rps.push_back(std::make_pair(r1, r2)); }
802	>	{ rps.push_back(make_pair(r1, r2)); }
803		if (r1 > 1e-6) { //r1 non-negative
804		r1 = -r1;
805		if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6)
806	<	{ rps.push_back(std::make_pair(r1, r2)); }
806	>	{ rps.push_back(make_pair(r1, r2)); }
807		}
808		}
809		}
810	<	// Consider combininig together the solving pair part w/ the searching
810	>	// Consider combining together the solving pair part w/ the searching
811		// best solution part so that we don't need the pairs vector
812		if (!rps.empty()) {
813		RealType smallestDiff = HONKING_LARGE_VALUE;
814		RealType diff;
815	<	std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0);
816	<	std::vector<std::pair<RealType,RealType> >::iterator rpi;
815	>	pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
816	>	vector<pair<RealType,RealType> >::iterator rpi;
817		for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
818		r1 = (*rpi).first;
819		r2 = (*rpi).second;
#	Line 762 | Line 832 | namespace OpenMD {
832		+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
833		break;
834		case rnemdPzScale :
835	+	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
836	+	+ fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
837		default :
838		break;
839		}
#	Line 773 | Line 845 | namespace OpenMD {
845		#ifdef IS_MPI
846		if (worldRank == 0) {
847		#endif
848	<	std::cerr << "we choose r1 = " << bestPair.first
849	<	<< " and r2 = " << bestPair.second << "\n";
848	>	sprintf(painCave.errMsg,
849	>	"RNEMD: roots r1= %lf\tr2 = %lf\n",
850	>	bestPair.first, bestPair.second);
851	>	painCave.isFatal = 0;
852	>	painCave.severity = OPENMD_INFO;
853	>	simError();
854		#ifdef IS_MPI
855		}
856		#endif
857	<
857	>
858		RealType x, y, z;
859	<	switch(rnemdType_) {
860	<	case rnemdKineticScale :
861	<	x = bestPair.first;
862	<	y = bestPair.first;
863	<	z = bestPair.second;
864	<	break;
865	<	case rnemdPxScale :
866	<	x = c;
867	<	y = bestPair.first;
868	<	z = bestPair.second;
869	<	break;
870	<	case rnemdPyScale :
871	<	x = bestPair.first;
872	<	y = c;
873	<	z = bestPair.second;
874	<	break;
875	<	case rnemdPzScale :
876	<	default :
877	<	break;
878	<	}
879	<	std::vector<StuntDouble*>::iterator sdi;
859	>	switch(rnemdType_) {
860	>	case rnemdKineticScale :
861	>	x = bestPair.first;
862	>	y = bestPair.first;
863	>	z = bestPair.second;
864	>	break;
865	>	case rnemdPxScale :
866	>	x = c;
867	>	y = bestPair.first;
868	>	z = bestPair.second;
869	>	break;
870	>	case rnemdPyScale :
871	>	x = bestPair.first;
872	>	y = c;
873	>	z = bestPair.second;
874	>	break;
875	>	case rnemdPzScale :
876	>	x = bestPair.first;
877	>	y = bestPair.second;
878	>	z = c;
879	>	break;
880	>	default :
881	>	break;
882	>	}
883	>	vector<StuntDouble*>::iterator sdi;
884		Vector3d vel;
885		for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
886		vel = (*sdi)->getVel();
#	Line 823 | Line 903 | namespace OpenMD {
903		exchangeSum_ += targetFlux_;
904		//we may want to check whether the exchange has been successful
905		} else {
906	<	std::cerr << "exchange NOT performed!\n";
906	>	sprintf(painCave.errMsg,
907	>	"RNEMD: exchange NOT performed!\n");
908	>	painCave.isFatal = 0;
909	>	painCave.severity = OPENMD_INFO;
910	>	simError();
911		failTrialCount_++;
912		}
913
#	Line 861 | Line 945 | namespace OpenMD {
945		StuntDouble* sd;
946		int idx;
947
948	+	// alternative approach, track all molecules instead of only those
949	+	// selected for scaling/swapping:
950	+	/*
951	+	SimInfo::MoleculeIterator miter;
952	+	vector<StuntDouble*>::iterator iiter;
953	+	Molecule* mol;
954	+	StuntDouble* integrableObject;
955	+	for (mol = info_->beginMolecule(miter); mol != NULL;
956	+	mol = info_->nextMolecule(miter))
957	+	integrableObject is essentially sd
958	+	for (integrableObject = mol->beginIntegrableObject(iiter);
959	+	integrableObject != NULL;
960	+	integrableObject = mol->nextIntegrableObject(iiter))
961	+	*/
962		for (sd = seleMan_.beginSelected(selei); sd != NULL;
963		sd = seleMan_.nextSelected(selei)) {
964
#	Line 876 | Line 974 | namespace OpenMD {
974		// which bin is this stuntdouble in?
975		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
976
977	<	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
978	<
977	>	int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
978	>	rnemdLogWidth_;
979	>	// no symmetrization allowed due to arbitary rnemdLogWidth_ value
980	>	/*
981		if (rnemdLogWidth_ == midBin_ + 1)
982		if (binNo > midBin_)
983		binNo = nBins_ - binNo;
984	<
984	>	*/
985		RealType mass = sd->getMass();
986		Vector3d vel = sd->getVel();
987		RealType value;
#	Line 891 | Line 991 | namespace OpenMD {
991		case rnemdKineticSwap :
992		case rnemdKineticScale :
993
994	<	value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
895	<	vel[2]*vel[2]);
994	>	value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
995
996		valueCount_[binNo] += 3;
997		if (sd->isDirectional()) {
#	Line 905 | Line 1004 | namespace OpenMD {
1004		int k = (i + 2) % 3;
1005		value += angMom[j] * angMom[j] / I(j, j) +
1006		angMom[k] * angMom[k] / I(k, k);
1007	<
1007	>
1008		valueCount_[binNo] +=2;
1009	<
1009	>
1010		} else {
1011		value += angMom[0]*angMom[0]/I(0, 0)
1012		+ angMom[1]*angMom[1]/I(1, 1)
#	Line 916 | Line 1015 | namespace OpenMD {
1015		}
1016		}
1017		value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb;
1018	<
1018	>
1019		break;
1020		case rnemdPx :
1021		case rnemdPxScale :
1022		value = mass * vel[0];
1023		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;
1024		break;
1025		case rnemdPy :
1026		case rnemdPyScale :
#	Line 939 | Line 1029 | namespace OpenMD {
1029		break;
1030		case rnemdPz :
1031		case rnemdPzScale :
1032	<	value = mass * vel[2];
1032	>	value = pos.z(); //temporarily for homogeneous systems ONLY
1033		valueCount_[binNo]++;
1034		break;
1035		case rnemdUnknown :
1036		default :
1037	+	value = 1.0;
1038	+	valueCount_[binNo]++;
1039		break;
1040		}
1041		valueHist_[binNo] += value;
950	–	}
1042
1043	+	if (output3DTemp_) {
1044	+	xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
1045	+	/ PhysicalConstants::kb;
1046	+	yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1047	+	/ PhysicalConstants::kb;
1048	+	zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1049	+	/ PhysicalConstants::kb;
1050	+	xTempHist_[binNo] += xVal;
1051	+	yTempHist_[binNo] += yVal;
1052	+	zTempHist_[binNo] += zVal;
1053	+	xyzTempCount_[binNo]++;
1054	+	}
1055	+	}
1056		}
1057
1058		void RNEMD::getStarted() {
1059	<	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1060	<	Stats& stat = currentSnap_->statData;
1061	<	stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1059	>	collectData();
1060	>	/* now should be able to output profile in step 0, but might not be useful
1061	>	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1062	>	Stats& stat = currentSnap_->statData;
1063	>	stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1064	>	*/
1065	>	getStatus();
1066		}
1067
1068		void RNEMD::getStatus() {
#	Line 965 | Line 1073 | namespace OpenMD {
1073
1074		stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1075		//or to be more meaningful, define another item as exchangeSum_ / time
1076	+	int j;
1077
969	–
1078		#ifdef IS_MPI
1079
1080		// all processors have the same number of bins, and STL vectors pack their
#	Line 976 | Line 1084 | namespace OpenMD {
1084		rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1085		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
1086		rnemdLogWidth_, MPI::INT, MPI::SUM);
1087	<
1087	>	if (output3DTemp_) {
1088	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
1089	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1090	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1091	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1092	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1093	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1094	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0],
1095	>	rnemdLogWidth_, MPI::INT, MPI::SUM);
1096	>	}
1097		// If we're the root node, should we print out the results
1098		int worldRank = MPI::COMM_WORLD.Get_rank();
1099		if (worldRank == 0) {
1100		#endif
984	–	int j;
1101		rnemdLog_ << time;
1102		for (j = 0; j < rnemdLogWidth_; j++) {
1103		rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
988	–	valueHist_[j] = 0.0;
1104		}
1105		rnemdLog_ << "\n";
1106	<	if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) {
1106	>	if (output3DTemp_) {
1107		xTempLog_ << time;
1108		for (j = 0; j < rnemdLogWidth_; j++) {
1109	<	xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j];
995	<	xTempHist_[j] = 0.0;
1109	>	xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j];
1110		}
1111		xTempLog_ << "\n";
1112		yTempLog_ << time;
1113		for (j = 0; j < rnemdLogWidth_; j++) {
1114	<	yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j];
1001	<	yTempHist_[j] = 0.0;
1114	>	yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1115		}
1116		yTempLog_ << "\n";
1117		zTempLog_ << time;
1118		for (j = 0; j < rnemdLogWidth_; j++) {
1119	<	zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j];
1007	<	zTempHist_[j] = 0.0;
1119	>	zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1120		}
1121		zTempLog_ << "\n";
1122		}
1011	–	for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0;
1123		#ifdef IS_MPI
1124	<	}
1124	>	}
1125		#endif
1126	<
1127	<
1126	>	for (j = 0; j < rnemdLogWidth_; j++) {
1127	>	valueCount_[j] = 0;
1128	>	valueHist_[j] = 0.0;
1129	>	}
1130	>	if (output3DTemp_)
1131	>	for (j = 0; j < rnemdLogWidth_; j++) {
1132	>	xTempHist_[j] = 0.0;
1133	>	yTempHist_[j] = 0.0;
1134	>	zTempHist_[j] = 0.0;
1135	>	xyzTempCount_[j] = 0;
1136	>	}
1137		}
1018	–
1138		}

Comparing trunk/src/integrators/RNEMD.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1561 by gezelter, Wed May 11 19:04:40 2011 UTC

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

Diff Legend

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