ViewVC Help

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

root/OpenMD/branches/development/src/rnemd/RNEMD.cpp

Comparing branches/development/src/integrators/RNEMD.cpp (file contents):
Revision 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC vs.
Revision 1722 by gezelter, Thu May 24 14:23:40 2012 UTC

#	Line 42 | Line 42
42		#include <cmath>
43		#include "integrators/RNEMD.hpp"
44		#include "math/Vector3.hpp"
45	+	#include "math/Vector.hpp"
46		#include "math/SquareMatrix3.hpp"
47		#include "math/Polynomial.hpp"
48		#include "primitives/Molecule.hpp"
#	Line 52 | Line 53
53		#ifndef IS_MPI
54		#include "math/SeqRandNumGen.hpp"
55		#else
55	–	#include <mpi.h>
56		#include "math/ParallelRandNumGen.hpp"
57		#endif
58
#	Line 72 | Line 72 | namespace OpenMD {
72
73		stringToEnumMap_["KineticSwap"] = rnemdKineticSwap;
74		stringToEnumMap_["KineticScale"] = rnemdKineticScale;
75	+	stringToEnumMap_["KineticScaleVAM"] = rnemdKineticScaleVAM;
76	+	stringToEnumMap_["KineticScaleAM"] = rnemdKineticScaleAM;
77		stringToEnumMap_["PxScale"] = rnemdPxScale;
78		stringToEnumMap_["PyScale"] = rnemdPyScale;
79		stringToEnumMap_["PzScale"] = rnemdPzScale;
80		stringToEnumMap_["Px"] = rnemdPx;
81		stringToEnumMap_["Py"] = rnemdPy;
82		stringToEnumMap_["Pz"] = rnemdPz;
83	+	stringToEnumMap_["ShiftScaleV"] = rnemdShiftScaleV;
84	+	stringToEnumMap_["ShiftScaleVAM"] = rnemdShiftScaleVAM;
85		stringToEnumMap_["Unknown"] = rnemdUnknown;
86
87		rnemdObjectSelection_ = simParams->getRNEMD_objectSelection();
#	Line 121 | Line 125 | namespace OpenMD {
125		simError();
126		}
127
128	+	outputTemp_ = false;
129	+	if (simParams->haveRNEMD_outputTemperature()) {
130	+	outputTemp_ = simParams->getRNEMD_outputTemperature();
131	+	} else if ((rnemdType_ == rnemdKineticSwap) ||
132	+	(rnemdType_ == rnemdKineticScale) ||
133	+	(rnemdType_ == rnemdKineticScaleVAM) ||
134	+	(rnemdType_ == rnemdKineticScaleAM)) {
135	+	outputTemp_ = true;
136	+	}
137	+	outputVx_ = false;
138	+	if (simParams->haveRNEMD_outputVx()) {
139	+	outputVx_ = simParams->getRNEMD_outputVx();
140	+	} else if ((rnemdType_ == rnemdPx) || (rnemdType_ == rnemdPxScale)) {
141	+	outputVx_ = true;
142	+	}
143	+	outputVy_ = false;
144	+	if (simParams->haveRNEMD_outputVy()) {
145	+	outputVy_ = simParams->getRNEMD_outputVy();
146	+	} else if ((rnemdType_ == rnemdPy) || (rnemdType_ == rnemdPyScale)) {
147	+	outputVy_ = true;
148	+	}
149		output3DTemp_ = false;
150	<	if (simParams->haveRNEMD_outputDimensionalTemperature()) {
151	<	output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature();
150	>	if (simParams->haveRNEMD_outputXyzTemperature()) {
151	>	output3DTemp_ = simParams->getRNEMD_outputXyzTemperature();
152		}
153	+	outputRotTemp_ = false;
154	+	if (simParams->haveRNEMD_outputRotTemperature()) {
155	+	outputRotTemp_ = simParams->getRNEMD_outputRotTemperature();
156	+	}
157
158		#ifdef IS_MPI
159		if (worldRank == 0) {
160		#endif
161
162	+	//may have rnemdWriter separately
163		string rnemdFileName;
164	<	switch(rnemdType_) {
165	<	case rnemdKineticSwap :
136	<	case rnemdKineticScale :
164	>
165	>	if (outputTemp_) {
166		rnemdFileName = "temperature.log";
167	<	break;
139	<	case rnemdPx :
140	<	case rnemdPxScale :
141	<	case rnemdPy :
142	<	case rnemdPyScale :
143	<	rnemdFileName = "momemtum.log";
144	<	break;
145	<	case rnemdPz :
146	<	case rnemdPzScale :
147	<	case rnemdUnknown :
148	<	default :
149	<	rnemdFileName = "rnemd.log";
150	<	break;
167	>	tempLog_.open(rnemdFileName.c_str());
168		}
169	<	rnemdLog_.open(rnemdFileName.c_str());
169	>	if (outputVx_) {
170	>	rnemdFileName = "velocityX.log";
171	>	vxzLog_.open(rnemdFileName.c_str());
172	>	}
173	>	if (outputVy_) {
174	>	rnemdFileName = "velocityY.log";
175	>	vyzLog_.open(rnemdFileName.c_str());
176	>	}
177
154	–	string xTempFileName;
155	–	string yTempFileName;
156	–	string zTempFileName;
178		if (output3DTemp_) {
179	<	xTempFileName = "temperatureX.log";
180	<	yTempFileName = "temperatureY.log";
181	<	zTempFileName = "temperatureZ.log";
182	<	xTempLog_.open(xTempFileName.c_str());
183	<	yTempLog_.open(yTempFileName.c_str());
184	<	zTempLog_.open(zTempFileName.c_str());
179	>	rnemdFileName = "temperatureX.log";
180	>	xTempLog_.open(rnemdFileName.c_str());
181	>	rnemdFileName = "temperatureY.log";
182	>	yTempLog_.open(rnemdFileName.c_str());
183	>	rnemdFileName = "temperatureZ.log";
184	>	zTempLog_.open(rnemdFileName.c_str());
185		}
186	+	if (outputRotTemp_) {
187	+	rnemdFileName = "temperatureR.log";
188	+	rotTempLog_.open(rnemdFileName.c_str());
189	+	}
190
191		#ifdef IS_MPI
192		}
#	Line 179 | Line 204 | namespace OpenMD {
204		} else {
205		zShift_ = 0.0;
206		}
207	<	//cerr << "we have zShift_ = " << zShift_ << "\n";
208	<	//shift slabs by half slab width, might be useful in heterogeneous systems
209	<	//set to 0.0 if not using it; can NOT be used in status output yet
207	>	//cerr << "I shift slabs by " << zShift_ << " Lz\n";
208	>	//shift slabs by half slab width, maybe useful in heterogeneous systems
209	>	//set to 0.0 if not using it; N/A in status output yet
210		if (simParams->haveRNEMD_logWidth()) {
211		set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
212	<	/*arbitary rnemdLogWidth_ no checking
213	<	if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
212	>	/*arbitary rnemdLogWidth_, no checking;
213	>	if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
214		cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
215		cerr << "Automaically set back to default.\n";
216		rnemdLogWidth_ = nBins_;
217	<	}*/
217	>	}*/
218		} else {
219		set_RNEMD_logWidth(nBins_);
220		}
221	<	valueHist_.resize(rnemdLogWidth_, 0.0);
222	<	valueCount_.resize(rnemdLogWidth_, 0);
221	>	tempHist_.resize(rnemdLogWidth_, 0.0);
222	>	tempCount_.resize(rnemdLogWidth_, 0);
223	>	pxzHist_.resize(rnemdLogWidth_, 0.0);
224	>	//vxzCount_.resize(rnemdLogWidth_, 0);
225	>	pyzHist_.resize(rnemdLogWidth_, 0.0);
226	>	//vyzCount_.resize(rnemdLogWidth_, 0);
227	>
228	>	mHist_.resize(rnemdLogWidth_, 0.0);
229		xTempHist_.resize(rnemdLogWidth_, 0.0);
230		yTempHist_.resize(rnemdLogWidth_, 0.0);
231		zTempHist_.resize(rnemdLogWidth_, 0.0);
232		xyzTempCount_.resize(rnemdLogWidth_, 0);
233	+	rotTempHist_.resize(rnemdLogWidth_, 0.0);
234	+	rotTempCount_.resize(rnemdLogWidth_, 0);
235
236		set_RNEMD_exchange_total(0.0);
237		if (simParams->haveRNEMD_targetFlux()) {
#	Line 206 | Line 239 | namespace OpenMD {
239		} else {
240		set_RNEMD_target_flux(0.0);
241		}
242	+	if (simParams->haveRNEMD_targetJzKE()) {
243	+	set_RNEMD_target_JzKE(simParams->getRNEMD_targetJzKE());
244	+	} else {
245	+	set_RNEMD_target_JzKE(0.0);
246	+	}
247	+	if (simParams->haveRNEMD_targetJzpx()) {
248	+	set_RNEMD_target_jzpx(simParams->getRNEMD_targetJzpx());
249	+	} else {
250	+	set_RNEMD_target_jzpx(0.0);
251	+	}
252	+	jzp_.x() = targetJzpx_;
253	+	njzp_.x() = -targetJzpx_;
254	+	if (simParams->haveRNEMD_targetJzpy()) {
255	+	set_RNEMD_target_jzpy(simParams->getRNEMD_targetJzpy());
256	+	} else {
257	+	set_RNEMD_target_jzpy(0.0);
258	+	}
259	+	jzp_.y() = targetJzpy_;
260	+	njzp_.y() = -targetJzpy_;
261	+	if (simParams->haveRNEMD_targetJzpz()) {
262	+	set_RNEMD_target_jzpz(simParams->getRNEMD_targetJzpz());
263	+	} else {
264	+	set_RNEMD_target_jzpz(0.0);
265	+	}
266	+	jzp_.z() = targetJzpz_;
267	+	njzp_.z() = -targetJzpz_;
268
269		#ifndef IS_MPI
270		if (simParams->haveSeed()) {
#	Line 238 | Line 297 | namespace OpenMD {
297		painCave.severity = OPENMD_INFO;
298		simError();
299
300	<	rnemdLog_.close();
301	<	if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) {
300	>	if (outputTemp_) tempLog_.close();
301	>	if (outputVx_)   vxzLog_.close();
302	>	if (outputVy_)   vyzLog_.close();
303	>
304	>	if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale ||
305	>	rnemdType_ == rnemdPyScale) {
306		sprintf(painCave.errMsg,
307		"RNEMD: total root-checking warnings: %d\n",
308		failRootCount_);
#	Line 252 | Line 315 | namespace OpenMD {
315		yTempLog_.close();
316		zTempLog_.close();
317		}
318	+	if (outputRotTemp_) rotTempLog_.close();
319	+
320		#ifdef IS_MPI
321		}
322		#endif
#	Line 304 | Line 369 | namespace OpenMD {
369		switch(rnemdType_) {
370		case rnemdKineticSwap :
371
372	<	value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
373	<	vel[2]*vel[2]);
374	<	/*
310	<	if (sd->isDirectional()) {
372	>	value = mass * vel.lengthSquare();
373	>
374	>	if (sd->isDirectional()) {
375		Vector3d angMom = sd->getJ();
376		Mat3x3d I = sd->getI();
377
378		if (sd->isLinear()) {
379	<	int i = sd->linearAxis();
380	<	int j = (i + 1) % 3;
381	<	int k = (i + 2) % 3;
382	<	value += angMom[j] * angMom[j] / I(j, j) +
383	<	angMom[k] * angMom[k] / I(k, k);
384	<	} else {
385	<	value += angMom[0]*angMom[0]/I(0, 0)
386	<	+ angMom[1]*angMom[1]/I(1, 1)
387	<	+ angMom[2]*angMom[2]/I(2, 2);
379	>	int i = sd->linearAxis();
380	>	int j = (i + 1) % 3;
381	>	int k = (i + 2) % 3;
382	>	value += angMom[j] * angMom[j] / I(j, j) +
383	>	angMom[k] * angMom[k] / I(k, k);
384	>	} else {
385	>	value += angMom[0]*angMom[0]/I(0, 0)
386	>	+ angMom[1]*angMom[1]/I(1, 1)
387	>	+ angMom[2]*angMom[2]/I(2, 2);
388		}
389	<	} no exchange of angular momenta
390	<	*/
389	>	} //angular momenta exchange enabled
390	>	//energyConvert temporarily disabled
391		//make exchangeSum_ comparable between swap & scale
328	–	//temporarily without using energyConvert
392		//value = value * 0.5 / PhysicalConstants::energyConvert;
393		value *= 0.5;
394		break;
#	Line 381 | Line 444 | namespace OpenMD {
444		MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
445		// Even if we didn't find a maximum, did someone else?
446		MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
447	<	struct {
448	<	RealType val;
449	<	int rank;
450	<	} max_vals, min_vals;
447	>	#endif
448	>
449	>	if (max_found && min_found) {
450	>
451	>	#ifdef IS_MPI
452	>	struct {
453	>	RealType val;
454	>	int rank;
455	>	} max_vals, min_vals;
456
457	<	if (min_found) {
390	<	if (my_min_found)
457	>	if (my_min_found) {
458		min_vals.val = min_val;
459	<	else
459	>	} else {
460		min_vals.val = HONKING_LARGE_VALUE;
461	<
461	>	}
462		min_vals.rank = worldRank;
463
464		// Who had the minimum?
465		MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
466		1, MPI::REALTYPE_INT, MPI::MINLOC);
467		min_val = min_vals.val;
401	–	}
468
469	<	if (max_found) {
404	<	if (my_max_found)
469	>	if (my_max_found) {
470		max_vals.val = max_val;
471	<	else
471	>	} else {
472		max_vals.val = -HONKING_LARGE_VALUE;
473	<
473	>	}
474		max_vals.rank = worldRank;
475
476		// Who had the maximum?
477		MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
478		1, MPI::REALTYPE_INT, MPI::MAXLOC);
479		max_val = max_vals.val;
415	–	}
480		#endif
481	<
418	<	if (max_found && min_found) {
481	>
482		if (min_val < max_val) {
483	<
483	>
484		#ifdef IS_MPI
485		if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
486		// I have both maximum and minimum, so proceed like a single
487		// processor version:
488		#endif
489	<	// objects to be swapped: velocity ONLY
489	>
490		Vector3d min_vel = min_sd->getVel();
491		Vector3d max_vel = max_sd->getVel();
492		RealType temp_vel;
#	Line 432 | Line 495 | namespace OpenMD {
495		case rnemdKineticSwap :
496		min_sd->setVel(max_vel);
497		max_sd->setVel(min_vel);
498	<	/*
436	<	if (min_sd->isDirectional() && max_sd->isDirectional()) {
498	>	if (min_sd->isDirectional() && max_sd->isDirectional()) {
499		Vector3d min_angMom = min_sd->getJ();
500		Vector3d max_angMom = max_sd->getJ();
501		min_sd->setJ(max_angMom);
502		max_sd->setJ(min_angMom);
503	<	} no angular momentum exchange
504	<	*/
503	>	}//angular momenta exchange enabled
504	>	//assumes same rigid body identity
505		break;
506		case rnemdPx :
507		temp_vel = min_vel.x();
#	Line 465 | Line 527 | namespace OpenMD {
527		default :
528		break;
529		}
530	+
531		#ifdef IS_MPI
532		// the rest of the cases only apply in parallel simulations:
533		} else if (max_vals.rank == worldRank) {
#	Line 483 | Line 546 | namespace OpenMD {
546		switch(rnemdType_) {
547		case rnemdKineticSwap :
548		max_sd->setVel(min_vel);
549	<	//no angular momentum exchange for now
487	<	/*
549	>	//angular momenta exchange enabled
550		if (max_sd->isDirectional()) {
551		Vector3d min_angMom;
552		Vector3d max_angMom = max_sd->getJ();
#	Line 497 | Line 559 | namespace OpenMD {
559		status);
560
561		max_sd->setJ(min_angMom);
562	<	}
501	<	*/
562	>	}
563		break;
564		case rnemdPx :
565		max_vel.x() = min_vel.x();
#	Line 531 | Line 592 | namespace OpenMD {
592		switch(rnemdType_) {
593		case rnemdKineticSwap :
594		min_sd->setVel(max_vel);
595	<	// no angular momentum exchange for now
535	<	/*
595	>	//angular momenta exchange enabled
596		if (min_sd->isDirectional()) {
597		Vector3d min_angMom = min_sd->getJ();
598		Vector3d max_angMom;
#	Line 546 | Line 606 | namespace OpenMD {
606
607		min_sd->setJ(max_angMom);
608		}
549	–	*/
609		break;
610		case rnemdPx :
611		min_vel.x() = max_vel.x();
#	Line 576 | Line 635 | namespace OpenMD {
635		}
636		} else {
637		sprintf(painCave.errMsg,
638	<	"RNEMD: exchange NOT performed because at least one\n"
639	<	"\tof the two slabs is empty\n");
638	>	"RNEMD: exchange NOT performed because selected object\n"
639	>	"\tnot present in at least one of the two slabs.\n");
640		painCave.isFatal = 0;
641		painCave.severity = OPENMD_INFO;
642		simError();
#	Line 605 | Line 664 | namespace OpenMD {
664		RealType Khx = 0.0;
665		RealType Khy = 0.0;
666		RealType Khz = 0.0;
667	+	RealType Khw = 0.0;
668		RealType Pcx = 0.0;
669		RealType Pcy = 0.0;
670		RealType Pcz = 0.0;
671		RealType Kcx = 0.0;
672		RealType Kcy = 0.0;
673		RealType Kcz = 0.0;
674	+	RealType Kcw = 0.0;
675
676		for (sd = seleMan_.beginSelected(selei); sd != NULL;
677		sd = seleMan_.nextSelected(selei)) {
#	Line 643 | Line 704 | namespace OpenMD {
704		Khx += mass * vel.x() * vel.x();
705		Khy += mass * vel.y() * vel.y();
706		Khz += mass * vel.z() * vel.z();
707	+	//if (rnemdType_ == rnemdKineticScaleVAM) {
708	+	if (sd->isDirectional()) {
709	+	Vector3d angMom = sd->getJ();
710	+	Mat3x3d I = sd->getI();
711	+	if (sd->isLinear()) {
712	+	int i = sd->linearAxis();
713	+	int j = (i + 1) % 3;
714	+	int k = (i + 2) % 3;
715	+	Khw += angMom[j] * angMom[j] / I(j, j) +
716	+	angMom[k] * angMom[k] / I(k, k);
717	+	} else {
718	+	Khw += angMom[0]*angMom[0]/I(0, 0)
719	+	+ angMom[1]*angMom[1]/I(1, 1)
720	+	+ angMom[2]*angMom[2]/I(2, 2);
721	+	}
722	+	}
723	+	//}
724		} else { //midBin_
725		coldBin.push_back(sd);
726		Pcx += mass * vel.x();
#	Line 651 | Line 729 | namespace OpenMD {
729		Kcx += mass * vel.x() * vel.x();
730		Kcy += mass * vel.y() * vel.y();
731		Kcz += mass * vel.z() * vel.z();
732	+	//if (rnemdType_ == rnemdKineticScaleVAM) {
733	+	if (sd->isDirectional()) {
734	+	Vector3d angMom = sd->getJ();
735	+	Mat3x3d I = sd->getI();
736	+	if (sd->isLinear()) {
737	+	int i = sd->linearAxis();
738	+	int j = (i + 1) % 3;
739	+	int k = (i + 2) % 3;
740	+	Kcw += angMom[j] * angMom[j] / I(j, j) +
741	+	angMom[k] * angMom[k] / I(k, k);
742	+	} else {
743	+	Kcw += angMom[0]*angMom[0]/I(0, 0)
744	+	+ angMom[1]*angMom[1]/I(1, 1)
745	+	+ angMom[2]*angMom[2]/I(2, 2);
746	+	}
747	+	}
748	+	//}
749		}
750		}
751		}
752	<
752	>
753		Khx *= 0.5;
754		Khy *= 0.5;
755		Khz *= 0.5;
756	+	Khw *= 0.5;
757		Kcx *= 0.5;
758		Kcy *= 0.5;
759		Kcz *= 0.5;
760	+	Kcw *= 0.5;
761
762	+	std::cerr << "Khx= " << Khx << "\tKhy= " << Khy << "\tKhz= " << Khz
763	+	<< "\tKhw= " << Khw << "\tKcx= " << Kcx << "\tKcy= " << Kcy
764	+	<< "\tKcz= " << Kcz << "\tKcw= " << Kcw << "\n";
765	+	std::cerr << "Phx= " << Phx << "\tPhy= " << Phy << "\tPhz= " << Phz
766	+	<< "\tPcx= " << Pcx << "\tPcy= " << Pcy << "\tPcz= " <<Pcz<<"\n";
767	+
768		#ifdef IS_MPI
769		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
770		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
#	Line 673 | Line 776 | namespace OpenMD {
776		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
777		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
778		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
779	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM);
780	+
781		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
782		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
783		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
784	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM);
785		#endif
786
787	<	//use coldBin coeff's
787	>	//solve coldBin coeff's first
788		RealType px = Pcx / Phx;
789		RealType py = Pcy / Phy;
790		RealType pz = Pcz / Phz;
791	+	RealType c, x, y, z;
792	+	bool successfulScale = false;
793	+	if ((rnemdType_ == rnemdKineticScaleVAM) ||
794	+	(rnemdType_ == rnemdKineticScaleAM)) {
795	+	//may need sanity check Khw & Kcw > 0
796
797	<	RealType a000, a110, c0, a001, a111, b01, b11, c1, c;
798	<	switch(rnemdType_) {
799	<	case rnemdKineticScale :
800	<	// used hotBin coeff's & only scale x & y dimensions
801	<	/*
691	<	RealType px = Phx / Pcx;
692	<	RealType py = Phy / Pcy;
693	<	a110 = Khy;
694	<	c0 = - Khx - Khy - targetFlux_;
695	<	a000 = Khx;
696	<	a111 = Kcy * py * py;
697	<	b11 = -2.0 * Kcy * py * (1.0 + py);
698	<	c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_;
699	<	b01 = -2.0 * Kcx * px * (1.0 + px);
700	<	a001 = Kcx * px * px;
701	<	*/
702	<	//scale all three dimensions, let c_x = c_y
703	<	a000 = Kcx + Kcy;
704	<	a110 = Kcz;
705	<	c0 = targetFlux_ - Kcx - Kcy - Kcz;
706	<	a001 = Khx * px * px + Khy * py * py;
707	<	a111 = Khz * pz * pz;
708	<	b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
709	<	b11 = -2.0 * Khz * pz * (1.0 + pz);
710	<	c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
711	<	+ Khz * pz * (2.0 + pz) - targetFlux_;
712	<	break;
713	<	case rnemdPxScale :
714	<	c = 1 - targetFlux_ / Pcx;
715	<	a000 = Kcy;
716	<	a110 = Kcz;
717	<	c0 = Kcx * c * c - Kcx - Kcy - Kcz;
718	<	a001 = py * py * Khy;
719	<	a111 = pz * pz * Khz;
720	<	b01 = -2.0 * Khy * py * (1.0 + py);
721	<	b11 = -2.0 * Khz * pz * (1.0 + pz);
722	<	c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
723	<	+ Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
724	<	break;
725	<	case rnemdPyScale :
726	<	c = 1 - targetFlux_ / Pcy;
727	<	a000 = Kcx;
728	<	a110 = Kcz;
729	<	c0 = Kcy * c * c - Kcx - Kcy - Kcz;
730	<	a001 = px * px * Khx;
731	<	a111 = pz * pz * Khz;
732	<	b01 = -2.0 * Khx * px * (1.0 + px);
733	<	b11 = -2.0 * Khz * pz * (1.0 + pz);
734	<	c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
735	<	+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
736	<	break;
737	<	case rnemdPzScale ://we don't really do this, do we?
738	<	c = 1 - targetFlux_ / Pcz;
739	<	a000 = Kcx;
740	<	a110 = Kcy;
741	<	c0 = Kcz * c * c - Kcx - Kcy - Kcz;
742	<	a001 = px * px * Khx;
743	<	a111 = py * py * Khy;
744	<	b01 = -2.0 * Khx * px * (1.0 + px);
745	<	b11 = -2.0 * Khy * py * (1.0 + py);
746	<	c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
747	<	+ Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
748	<	break;
749	<	default :
750	<	break;
751	<	}
797	>	if (rnemdType_ == rnemdKineticScaleVAM) {
798	>	c = 1.0 - targetFlux_ / (Kcx + Kcy + Kcz + Kcw);
799	>	} else {
800	>	c = 1.0 - targetFlux_ / Kcw;
801	>	}
802
803	<	RealType v1 = a000 * a111 - a001 * a110;
804	<	RealType v2 = a000 * b01;
805	<	RealType v3 = a000 * b11;
806	<	RealType v4 = a000 * c1 - a001 * c0;
807	<	RealType v8 = a110 * b01;
808	<	RealType v10 = - b01 * c0;
809	<
810	<	RealType u0 = v2 * v10 - v4 * v4;
811	<	RealType u1 = -2.0 * v3 * v4;
812	<	RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
813	<	RealType u3 = -2.0 * v1 * v3;
814	<	RealType u4 = - v1 * v1;
815	<	//rescale coefficients
816	<	RealType maxAbs = fabs(u0);
817	<	if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
818	<	if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
819	<	if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
820	<	if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
821	<	u0 /= maxAbs;
822	<	u1 /= maxAbs;
823	<	u2 /= maxAbs;
824	<	u3 /= maxAbs;
825	<	u4 /= maxAbs;
826	<	//max_element(start, end) is also available.
827	<	Polynomial<RealType> poly; //same as DoublePolynomial poly;
828	<	poly.setCoefficient(4, u4);
829	<	poly.setCoefficient(3, u3);
830	<	poly.setCoefficient(2, u2);
831	<	poly.setCoefficient(1, u1);
832	<	poly.setCoefficient(0, u0);
833	<	vector<RealType> realRoots = poly.FindRealRoots();
834	<
835	<	vector<RealType>::iterator ri;
836	<	RealType r1, r2, alpha0;
837	<	vector<pair<RealType,RealType> > rps;
838	<	for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
839	<	r2 = *ri;
840	<	//check if FindRealRoots() give the right answer
841	<	if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
842	<	sprintf(painCave.errMsg,
843	<	"RNEMD Warning: polynomial solve seems to have an error!");
844	<	painCave.isFatal = 0;
845	<	simError();
846	<	failRootCount_++;
803	>	if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients
804	>	c = sqrt(c);
805	>	std::cerr << "cold slab scaling coefficient: " << c << endl;
806	>	//now convert to hotBin coefficient
807	>	RealType w = 0.0;
808	>	if (rnemdType_ ==  rnemdKineticScaleVAM) {
809	>	x = 1.0 + px * (1.0 - c);
810	>	y = 1.0 + py * (1.0 - c);
811	>	z = 1.0 + pz * (1.0 - c);
812	>	/* more complicated way
813	>	w = 1.0 + (Kcw - Kcw * c * c - (c * c * (Kcx + Kcy + Kcz
814	>	+ Khx * px * px + Khy * py * py + Khz * pz * pz)
815	>	- 2.0 * c * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py)
816	>	+ Khz * pz * (1.0 + pz)) + Khx * px * (2.0 + px)
817	>	+ Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
818	>	- Kcx - Kcy - Kcz)) / Khw; the following is simpler
819	>	*/
820	>	if ((fabs(x - 1.0) < 0.1) && (fabs(y - 1.0) < 0.1) &&
821	>	(fabs(z - 1.0) < 0.1)) {
822	>	w = 1.0 + (targetFlux_ + Khx * (1.0 - x * x) + Khy * (1.0 - y * y)
823	>	+ Khz * (1.0 - z * z)) / Khw;
824	>	}//no need to calculate w if x, y or z is out of range
825	>	} else {
826	>	w = 1.0 + targetFlux_ / Khw;
827	>	}
828	>	if ((w > 0.81) && (w < 1.21)) {//restrict scaling coefficients
829	>	//if w is in the right range, so should be x, y, z.
830	>	vector<StuntDouble*>::iterator sdi;
831	>	Vector3d vel;
832	>	for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
833	>	if (rnemdType_ == rnemdKineticScaleVAM) {
834	>	vel = (*sdi)->getVel() * c;
835	>	//vel.x() *= c;
836	>	//vel.y() *= c;
837	>	//vel.z() *= c;
838	>	(*sdi)->setVel(vel);
839	>	}
840	>	if ((*sdi)->isDirectional()) {
841	>	Vector3d angMom = (*sdi)->getJ() * c;
842	>	//angMom[0] *= c;
843	>	//angMom[1] *= c;
844	>	//angMom[2] *= c;
845	>	(*sdi)->setJ(angMom);
846	>	}
847	>	}
848	>	w = sqrt(w);
849	>	std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
850	>	<< "\twh= " << w << endl;
851	>	for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
852	>	if (rnemdType_ == rnemdKineticScaleVAM) {
853	>	vel = (*sdi)->getVel();
854	>	vel.x() *= x;
855	>	vel.y() *= y;
856	>	vel.z() *= z;
857	>	(*sdi)->setVel(vel);
858	>	}
859	>	if ((*sdi)->isDirectional()) {
860	>	Vector3d angMom = (*sdi)->getJ() * w;
861	>	//angMom[0] *= w;
862	>	//angMom[1] *= w;
863	>	//angMom[2] *= w;
864	>	(*sdi)->setJ(angMom);
865	>	}
866	>	}
867	>	successfulScale = true;
868	>	exchangeSum_ += targetFlux_;
869	>	}
870		}
871	<	//might not be useful w/o rescaling coefficients
872	<	alpha0 = -c0 - a110 * r2 * r2;
800	<	if (alpha0 >= 0.0) {
801	<	r1 = sqrt(alpha0 / a000);
802	<	if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6)
803	<	{ rps.push_back(make_pair(r1, r2)); }
804	<	if (r1 > 1e-6) { //r1 non-negative
805	<	r1 = -r1;
806	<	if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6)
807	<	{ rps.push_back(make_pair(r1, r2)); }
808	<	}
809	<	}
810	<	}
811	<	// Consider combining together the solving pair part w/ the searching
812	<	// best solution part so that we don't need the pairs vector
813	<	if (!rps.empty()) {
814	<	RealType smallestDiff = HONKING_LARGE_VALUE;
815	<	RealType diff;
816	<	pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
817	<	vector<pair<RealType,RealType> >::iterator rpi;
818	<	for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
819	<	r1 = (*rpi).first;
820	<	r2 = (*rpi).second;
821	<	switch(rnemdType_) {
822	<	case rnemdKineticScale :
823	<	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
824	<	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2)
825	<	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
826	<	break;
827	<	case rnemdPxScale :
828	<	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
829	<	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
830	<	break;
831	<	case rnemdPyScale :
832	<	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
833	<	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
834	<	break;
835	<	case rnemdPzScale :
836	<	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
837	<	+ fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
838	<	default :
839	<	break;
840	<	}
841	<	if (diff < smallestDiff) {
842	<	smallestDiff = diff;
843	<	bestPair = *rpi;
844	<	}
845	<	}
846	<	#ifdef IS_MPI
847	<	if (worldRank == 0) {
848	<	#endif
849	<	sprintf(painCave.errMsg,
850	<	"RNEMD: roots r1= %lf\tr2 = %lf\n",
851	<	bestPair.first, bestPair.second);
852	<	painCave.isFatal = 0;
853	<	painCave.severity = OPENMD_INFO;
854	<	simError();
855	<	#ifdef IS_MPI
856	<	}
857	<	#endif
858	<
859	<	RealType x, y, z;
871	>	} else {
872	>	RealType a000, a110, c0, a001, a111, b01, b11, c1;
873		switch(rnemdType_) {
874		case rnemdKineticScale :
875	<	x = bestPair.first;
876	<	y = bestPair.first;
877	<	z = bestPair.second;
878	<	break;
875	>	/* used hotBin coeff's & only scale x & y dimensions
876	>	RealType px = Phx / Pcx;
877	>	RealType py = Phy / Pcy;
878	>	a110 = Khy;
879	>	c0 = - Khx - Khy - targetFlux_;
880	>	a000 = Khx;
881	>	a111 = Kcy * py * py;
882	>	b11 = -2.0 * Kcy * py * (1.0 + py);
883	>	c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_;
884	>	b01 = -2.0 * Kcx * px * (1.0 + px);
885	>	a001 = Kcx * px * px;
886	>	*/
887	>	//scale all three dimensions, let c_x = c_y
888	>	a000 = Kcx + Kcy;
889	>	a110 = Kcz;
890	>	c0 = targetFlux_ - Kcx - Kcy - Kcz;
891	>	a001 = Khx * px * px + Khy * py * py;
892	>	a111 = Khz * pz * pz;
893	>	b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
894	>	b11 = -2.0 * Khz * pz * (1.0 + pz);
895	>	c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
896	>	+ Khz * pz * (2.0 + pz) - targetFlux_;
897	>	break;
898		case rnemdPxScale :
899	<	x = c;
900	<	y = bestPair.first;
901	<	z = bestPair.second;
902	<	break;
899	>	c = 1 - targetFlux_ / Pcx;
900	>	a000 = Kcy;
901	>	a110 = Kcz;
902	>	c0 = Kcx * c * c - Kcx - Kcy - Kcz;
903	>	a001 = py * py * Khy;
904	>	a111 = pz * pz * Khz;
905	>	b01 = -2.0 * Khy * py * (1.0 + py);
906	>	b11 = -2.0 * Khz * pz * (1.0 + pz);
907	>	c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
908	>	+ Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
909	>	break;
910		case rnemdPyScale :
911	<	x = bestPair.first;
912	<	y = c;
913	<	z = bestPair.second;
914	<	break;
915	<	case rnemdPzScale :
916	<	x = bestPair.first;
917	<	y = bestPair.second;
918	<	z = c;
919	<	break;
911	>	c = 1 - targetFlux_ / Pcy;
912	>	a000 = Kcx;
913	>	a110 = Kcz;
914	>	c0 = Kcy * c * c - Kcx - Kcy - Kcz;
915	>	a001 = px * px * Khx;
916	>	a111 = pz * pz * Khz;
917	>	b01 = -2.0 * Khx * px * (1.0 + px);
918	>	b11 = -2.0 * Khz * pz * (1.0 + pz);
919	>	c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
920	>	+ Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
921	>	break;
922	>	case rnemdPzScale ://we don't really do this, do we?
923	>	c = 1 - targetFlux_ / Pcz;
924	>	a000 = Kcx;
925	>	a110 = Kcy;
926	>	c0 = Kcz * c * c - Kcx - Kcy - Kcz;
927	>	a001 = px * px * Khx;
928	>	a111 = py * py * Khy;
929	>	b01 = -2.0 * Khx * px * (1.0 + px);
930	>	b11 = -2.0 * Khy * py * (1.0 + py);
931	>	c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
932	>	+ Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
933	>	break;
934		default :
935	<	break;
935	>	break;
936		}
937	<	vector<StuntDouble*>::iterator sdi;
938	<	Vector3d vel;
939	<	for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
940	<	vel = (*sdi)->getVel();
941	<	vel.x() *= x;
942	<	vel.y() *= y;
943	<	vel.z() *= z;
944	<	(*sdi)->setVel(vel);
937	>
938	>	RealType v1 = a000 * a111 - a001 * a110;
939	>	RealType v2 = a000 * b01;
940	>	RealType v3 = a000 * b11;
941	>	RealType v4 = a000 * c1 - a001 * c0;
942	>	RealType v8 = a110 * b01;
943	>	RealType v10 = - b01 * c0;
944	>
945	>	RealType u0 = v2 * v10 - v4 * v4;
946	>	RealType u1 = -2.0 * v3 * v4;
947	>	RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
948	>	RealType u3 = -2.0 * v1 * v3;
949	>	RealType u4 = - v1 * v1;
950	>	//rescale coefficients
951	>	RealType maxAbs = fabs(u0);
952	>	if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
953	>	if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
954	>	if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
955	>	if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
956	>	u0 /= maxAbs;
957	>	u1 /= maxAbs;
958	>	u2 /= maxAbs;
959	>	u3 /= maxAbs;
960	>	u4 /= maxAbs;
961	>	//max_element(start, end) is also available.
962	>	Polynomial<RealType> poly; //same as DoublePolynomial poly;
963	>	poly.setCoefficient(4, u4);
964	>	poly.setCoefficient(3, u3);
965	>	poly.setCoefficient(2, u2);
966	>	poly.setCoefficient(1, u1);
967	>	poly.setCoefficient(0, u0);
968	>	vector<RealType> realRoots = poly.FindRealRoots();
969	>
970	>	vector<RealType>::iterator ri;
971	>	RealType r1, r2, alpha0;
972	>	vector<pair<RealType,RealType> > rps;
973	>	for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) {
974	>	r2 = *ri;
975	>	//check if FindRealRoots() give the right answer
976	>	if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
977	>	sprintf(painCave.errMsg,
978	>	"RNEMD Warning: polynomial solve seems to have an error!");
979	>	painCave.isFatal = 0;
980	>	simError();
981	>	failRootCount_++;
982	>	}
983	>	//might not be useful w/o rescaling coefficients
984	>	alpha0 = -c0 - a110 * r2 * r2;
985	>	if (alpha0 >= 0.0) {
986	>	r1 = sqrt(alpha0 / a000);
987	>	if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111))
988	>	< 1e-6)
989	>	{ rps.push_back(make_pair(r1, r2)); }
990	>	if (r1 > 1e-6) { //r1 non-negative
991	>	r1 = -r1;
992	>	if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111))
993	>	< 1e-6)
994	>	{ rps.push_back(make_pair(r1, r2)); }
995	>	}
996	>	}
997		}
998	<	//convert to hotBin coefficient
999	<	x = 1.0 + px * (1.0 - x);
1000	<	y = 1.0 + py * (1.0 - y);
1001	<	z = 1.0 + pz * (1.0 - z);
1002	<	for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1003	<	vel = (*sdi)->getVel();
1004	<	vel.x() *= x;
1005	<	vel.y() *= y;
1006	<	vel.z() *= z;
1007	<	(*sdi)->setVel(vel);
998	>	// Consider combining together the solving pair part w/ the searching
999	>	// best solution part so that we don't need the pairs vector
1000	>	if (!rps.empty()) {
1001	>	RealType smallestDiff = HONKING_LARGE_VALUE;
1002	>	RealType diff;
1003	>	pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
1004	>	vector<pair<RealType,RealType> >::iterator rpi;
1005	>	for (rpi = rps.begin(); rpi != rps.end(); rpi++) {
1006	>	r1 = (*rpi).first;
1007	>	r2 = (*rpi).second;
1008	>	switch(rnemdType_) {
1009	>	case rnemdKineticScale :
1010	>	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1011	>	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2)
1012	>	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
1013	>	break;
1014	>	case rnemdPxScale :
1015	>	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1016	>	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
1017	>	break;
1018	>	case rnemdPyScale :
1019	>	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1020	>	+ fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
1021	>	break;
1022	>	case rnemdPzScale :
1023	>	diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1024	>	+ fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
1025	>	default :
1026	>	break;
1027	>	}
1028	>	if (diff < smallestDiff) {
1029	>	smallestDiff = diff;
1030	>	bestPair = *rpi;
1031	>	}
1032	>	}
1033	>	#ifdef IS_MPI
1034	>	if (worldRank == 0) {
1035	>	#endif
1036	>	sprintf(painCave.errMsg,
1037	>	"RNEMD: roots r1= %lf\tr2 = %lf\n",
1038	>	bestPair.first, bestPair.second);
1039	>	painCave.isFatal = 0;
1040	>	painCave.severity = OPENMD_INFO;
1041	>	simError();
1042	>	#ifdef IS_MPI
1043	>	}
1044	>	#endif
1045	>
1046	>	switch(rnemdType_) {
1047	>	case rnemdKineticScale :
1048	>	x = bestPair.first;
1049	>	y = bestPair.first;
1050	>	z = bestPair.second;
1051	>	break;
1052	>	case rnemdPxScale :
1053	>	x = c;
1054	>	y = bestPair.first;
1055	>	z = bestPair.second;
1056	>	break;
1057	>	case rnemdPyScale :
1058	>	x = bestPair.first;
1059	>	y = c;
1060	>	z = bestPair.second;
1061	>	break;
1062	>	case rnemdPzScale :
1063	>	x = bestPair.first;
1064	>	y = bestPair.second;
1065	>	z = c;
1066	>	break;
1067	>	default :
1068	>	break;
1069	>	}
1070	>	vector<StuntDouble*>::iterator sdi;
1071	>	Vector3d vel;
1072	>	for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1073	>	vel = (*sdi)->getVel();
1074	>	vel.x() *= x;
1075	>	vel.y() *= y;
1076	>	vel.z() *= z;
1077	>	(*sdi)->setVel(vel);
1078	>	}
1079	>	//convert to hotBin coefficient
1080	>	x = 1.0 + px * (1.0 - x);
1081	>	y = 1.0 + py * (1.0 - y);
1082	>	z = 1.0 + pz * (1.0 - z);
1083	>	for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1084	>	vel = (*sdi)->getVel();
1085	>	vel.x() *= x;
1086	>	vel.y() *= y;
1087	>	vel.z() *= z;
1088	>	(*sdi)->setVel(vel);
1089	>	}
1090	>	successfulScale = true;
1091	>	exchangeSum_ += targetFlux_;
1092		}
1093	<	exchangeSum_ += targetFlux_;
1094	<	//we may want to check whether the exchange has been successful
906	<	} else {
1093	>	}
1094	>	if (successfulScale != true) {
1095		sprintf(painCave.errMsg,
1096	<	"RNEMD: exchange NOT performed!\n");
1096	>	"RNEMD: exchange NOT performed!\n");
1097		painCave.isFatal = 0;
1098		painCave.severity = OPENMD_INFO;
1099		simError();
1100		failTrialCount_++;
1101		}
1102	+	}
1103
1104	+	void RNEMD::doShiftScale() {
1105	+
1106	+	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1107	+	Mat3x3d hmat = currentSnap_->getHmat();
1108	+
1109	+	seleMan_.setSelectionSet(evaluator_.evaluate());
1110	+
1111	+	int selei;
1112	+	StuntDouble* sd;
1113	+	int idx;
1114	+
1115	+	vector<StuntDouble*> hotBin, coldBin;
1116	+
1117	+	Vector3d Ph(V3Zero);
1118	+	RealType Mh = 0.0;
1119	+	RealType Kh = 0.0;
1120	+	Vector3d Pc(V3Zero);
1121	+	RealType Mc = 0.0;
1122	+	RealType Kc = 0.0;
1123	+
1124	+	for (sd = seleMan_.beginSelected(selei); sd != NULL;
1125	+	sd = seleMan_.nextSelected(selei)) {
1126	+
1127	+	idx = sd->getLocalIndex();
1128	+
1129	+	Vector3d pos = sd->getPos();
1130	+
1131	+	// wrap the stuntdouble's position back into the box:
1132	+
1133	+	if (usePeriodicBoundaryConditions_)
1134	+	currentSnap_->wrapVector(pos);
1135	+
1136	+	// which bin is this stuntdouble in?
1137	+	// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1138	+
1139	+	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
1140	+
1141	+	// if we're in bin 0 or the middleBin
1142	+	if (binNo == 0 || binNo == midBin_) {
1143	+
1144	+	RealType mass = sd->getMass();
1145	+	Vector3d vel = sd->getVel();
1146	+
1147	+	if (binNo == 0) {
1148	+	hotBin.push_back(sd);
1149	+	//std::cerr << "before, velocity = " << vel << endl;
1150	+	Ph += mass * vel;
1151	+	//std::cerr << "after, velocity = " << vel << endl;
1152	+	Mh += mass;
1153	+	Kh += mass * vel.lengthSquare();
1154	+	if (rnemdType_ == rnemdShiftScaleVAM) {
1155	+	if (sd->isDirectional()) {
1156	+	Vector3d angMom = sd->getJ();
1157	+	Mat3x3d I = sd->getI();
1158	+	if (sd->isLinear()) {
1159	+	int i = sd->linearAxis();
1160	+	int j = (i + 1) % 3;
1161	+	int k = (i + 2) % 3;
1162	+	Kh += angMom[j] * angMom[j] / I(j, j) +
1163	+	angMom[k] * angMom[k] / I(k, k);
1164	+	} else {
1165	+	Kh += angMom[0] * angMom[0] / I(0, 0) +
1166	+	angMom[1] * angMom[1] / I(1, 1) +
1167	+	angMom[2] * angMom[2] / I(2, 2);
1168	+	}
1169	+	}
1170	+	}
1171	+	} else { //midBin_
1172	+	coldBin.push_back(sd);
1173	+	Pc += mass * vel;
1174	+	Mc += mass;
1175	+	Kc += mass * vel.lengthSquare();
1176	+	if (rnemdType_ == rnemdShiftScaleVAM) {
1177	+	if (sd->isDirectional()) {
1178	+	Vector3d angMom = sd->getJ();
1179	+	Mat3x3d I = sd->getI();
1180	+	if (sd->isLinear()) {
1181	+	int i = sd->linearAxis();
1182	+	int j = (i + 1) % 3;
1183	+	int k = (i + 2) % 3;
1184	+	Kc += angMom[j] * angMom[j] / I(j, j) +
1185	+	angMom[k] * angMom[k] / I(k, k);
1186	+	} else {
1187	+	Kc += angMom[0] * angMom[0] / I(0, 0) +
1188	+	angMom[1] * angMom[1] / I(1, 1) +
1189	+	angMom[2] * angMom[2] / I(2, 2);
1190	+	}
1191	+	}
1192	+	}
1193	+	}
1194	+	}
1195	+	}
1196	+
1197	+	Kh *= 0.5;
1198	+	Kc *= 0.5;
1199	+
1200	+	std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1201	+	<< "\tKc= " << Kc << endl;
1202	+	std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1203	+
1204	+	#ifdef IS_MPI
1205	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1206	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
1207	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM);
1208	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM);
1209	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM);
1210	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM);
1211	+	#endif
1212	+
1213	+	bool successfulExchange = false;
1214	+	if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1215	+	Vector3d vc = Pc / Mc;
1216	+	Vector3d ac = njzp_ / Mc + vc;
1217	+	RealType cNumerator = Kc - targetJzKE_ - 0.5 * Mc * ac.lengthSquare();
1218	+	if (cNumerator > 0.0) {
1219	+	RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare();
1220	+	if (cDenominator > 0.0) {
1221	+	RealType c = sqrt(cNumerator / cDenominator);
1222	+	if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients
1223	+	Vector3d vh = Ph / Mh;
1224	+	Vector3d ah = jzp_ / Mh + vh;
1225	+	RealType hNumerator = Kh + targetJzKE_
1226	+	- 0.5 * Mh * ah.lengthSquare();
1227	+	if (hNumerator > 0.0) {
1228	+	RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare();
1229	+	if (hDenominator > 0.0) {
1230	+	RealType h = sqrt(hNumerator / hDenominator);
1231	+	if ((h > 0.9) && (h < 1.1)) {
1232	+	std::cerr << "cold slab scaling coefficient: " << c << "\n";
1233	+	std::cerr << "hot slab scaling coefficient: " << h << "\n";
1234	+	vector<StuntDouble*>::iterator sdi;
1235	+	Vector3d vel;
1236	+	for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1237	+	//vel = (*sdi)->getVel();
1238	+	vel = ((*sdi)->getVel() - vc) * c + ac;
1239	+	(*sdi)->setVel(vel);
1240	+	if (rnemdType_ == rnemdShiftScaleVAM) {
1241	+	if ((*sdi)->isDirectional()) {
1242	+	Vector3d angMom = (*sdi)->getJ() * c;
1243	+	(*sdi)->setJ(angMom);
1244	+	}
1245	+	}
1246	+	}
1247	+	for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1248	+	//vel = (*sdi)->getVel();
1249	+	vel = ((*sdi)->getVel() - vh) * h + ah;
1250	+	(*sdi)->setVel(vel);
1251	+	if (rnemdType_ == rnemdShiftScaleVAM) {
1252	+	if ((*sdi)->isDirectional()) {
1253	+	Vector3d angMom = (*sdi)->getJ() * h;
1254	+	(*sdi)->setJ(angMom);
1255	+	}
1256	+	}
1257	+	}
1258	+	successfulExchange = true;
1259	+	exchangeSum_ += targetFlux_;
1260	+	// this is a redundant variable for doShiftScale() so that
1261	+	// RNEMD can output one exchange quantity needed in a job.
1262	+	// need a better way to do this.
1263	+	}
1264	+	}
1265	+	}
1266	+	}
1267	+	}
1268	+	}
1269	+	}
1270	+	if (successfulExchange != true) {
1271	+	sprintf(painCave.errMsg,
1272	+	"RNEMD: exchange NOT performed!\n");
1273	+	painCave.isFatal = 0;
1274	+	painCave.severity = OPENMD_INFO;
1275	+	simError();
1276	+	failTrialCount_++;
1277	+	}
1278		}
1279
1280		void RNEMD::doRNEMD() {
1281
1282		switch(rnemdType_) {
1283		case rnemdKineticScale :
1284	+	case rnemdKineticScaleVAM :
1285	+	case rnemdKineticScaleAM :
1286		case rnemdPxScale :
1287		case rnemdPyScale :
1288		case rnemdPzScale :
#	Line 929 | Line 1294 | namespace OpenMD {
1294		case rnemdPz :
1295		doSwap();
1296		break;
1297	+	case rnemdShiftScaleV :
1298	+	case rnemdShiftScaleVAM :
1299	+	doShiftScale();
1300	+	break;
1301		case rnemdUnknown :
1302		default :
1303		break;
#	Line 977 | Line 1346 | namespace OpenMD {
1346
1347		int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
1348		rnemdLogWidth_;
1349	<	// no symmetrization allowed due to arbitary rnemdLogWidth_ value
1349	>	// no symmetrization allowed due to arbitary rnemdLogWidth_
1350		/*
1351		if (rnemdLogWidth_ == midBin_ + 1)
1352		if (binNo > midBin_)
1353		binNo = nBins_ - binNo;
1354		*/
1355		RealType mass = sd->getMass();
1356	+	mHist_[binNo] += mass;
1357		Vector3d vel = sd->getVel();
1358		RealType value;
1359	<	RealType xVal, yVal, zVal;
1359	>	//RealType xVal, yVal, zVal;
1360
1361	<	switch(rnemdType_) {
1362	<	case rnemdKineticSwap :
1363	<	case rnemdKineticScale :
994	<
995	<	value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
996	<
997	<	valueCount_[binNo] += 3;
1361	>	if (outputTemp_) {
1362	>	value = mass * vel.lengthSquare();
1363	>	tempCount_[binNo] += 3;
1364		if (sd->isDirectional()) {
1365		Vector3d angMom = sd->getJ();
1366		Mat3x3d I = sd->getI();
1001	–
1367		if (sd->isLinear()) {
1368		int i = sd->linearAxis();
1369		int j = (i + 1) % 3;
1370		int k = (i + 2) % 3;
1371		value += angMom[j] * angMom[j] / I(j, j) +
1372		angMom[k] * angMom[k] / I(k, k);
1373	<
1009	<	valueCount_[binNo] +=2;
1010	<
1373	>	tempCount_[binNo] +=2;
1374		} else {
1375	<	value += angMom[0]*angMom[0]/I(0, 0)
1376	<	+ angMom[1]*angMom[1]/I(1, 1)
1377	<	+ angMom[2]*angMom[2]/I(2, 2);
1378	<	valueCount_[binNo] +=3;
1375	>	value += angMom[0] * angMom[0] / I(0, 0) +
1376	>	angMom[1]*angMom[1]/I(1, 1) +
1377	>	angMom[2]*angMom[2]/I(2, 2);
1378	>	tempCount_[binNo] +=3;
1379		}
1380		}
1381	<	value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb;
1382	<
1383	<	break;
1384	<	case rnemdPx :
1385	<	case rnemdPxScale :
1381	>	value = value / PhysicalConstants::energyConvert
1382	>	/ PhysicalConstants::kb;//may move to getStatus()
1383	>	tempHist_[binNo] += value;
1384	>	}
1385	>	if (outputVx_) {
1386		value = mass * vel[0];
1387	<	valueCount_[binNo]++;
1388	<	break;
1389	<	case rnemdPy :
1390	<	case rnemdPyScale :
1387	>	//vxzCount_[binNo]++;
1388	>	pxzHist_[binNo] += value;
1389	>	}
1390	>	if (outputVy_) {
1391		value = mass * vel[1];
1392	<	valueCount_[binNo]++;
1393	<	break;
1031	<	case rnemdPz :
1032	<	case rnemdPzScale :
1033	<	value = pos.z(); //temporarily for homogeneous systems ONLY
1034	<	valueCount_[binNo]++;
1035	<	break;
1036	<	case rnemdUnknown :
1037	<	default :
1038	<	value = 1.0;
1039	<	valueCount_[binNo]++;
1040	<	break;
1392	>	//vyzCount_[binNo]++;
1393	>	pyzHist_[binNo] += value;
1394		}
1042	–	valueHist_[binNo] += value;
1395
1396		if (output3DTemp_) {
1397	<	xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert
1397	>	value = mass * vel.x() * vel.x();
1398	>	xTempHist_[binNo] += value;
1399	>	value = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1400		/ PhysicalConstants::kb;
1401	<	yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1401	>	yTempHist_[binNo] += value;
1402	>	value = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1403		/ PhysicalConstants::kb;
1404	<	zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1050	<	/ PhysicalConstants::kb;
1051	<	xTempHist_[binNo] += xVal;
1052	<	yTempHist_[binNo] += yVal;
1053	<	zTempHist_[binNo] += zVal;
1404	>	zTempHist_[binNo] += value;
1405		xyzTempCount_[binNo]++;
1406	+	}
1407	+	if (outputRotTemp_) {
1408	+	if (sd->isDirectional()) {
1409	+	Vector3d angMom = sd->getJ();
1410	+	Mat3x3d I = sd->getI();
1411	+	if (sd->isLinear()) {
1412	+	int i = sd->linearAxis();
1413	+	int j = (i + 1) % 3;
1414	+	int k = (i + 2) % 3;
1415	+	value = angMom[j] * angMom[j] / I(j, j) +
1416	+	angMom[k] * angMom[k] / I(k, k);
1417	+	rotTempCount_[binNo] +=2;
1418	+	} else {
1419	+	value = angMom[0] * angMom[0] / I(0, 0) +
1420	+	angMom[1] * angMom[1] / I(1, 1) +
1421	+	angMom[2] * angMom[2] / I(2, 2);
1422	+	rotTempCount_[binNo] +=3;
1423	+	}
1424	+	}
1425	+	value = value / PhysicalConstants::energyConvert
1426	+	/ PhysicalConstants::kb;//may move to getStatus()
1427	+	rotTempHist_[binNo] += value;
1428		}
1429	+
1430		}
1431		}
1432
1433		void RNEMD::getStarted() {
1434		collectData();
1435	<	/* now should be able to output profile in step 0, but might not be useful
1436	<	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1437	<	Stats& stat = currentSnap_->statData;
1438	<	stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1435	>	/*now can output profile in step 0, but might not be useful;
1436	>	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1437	>	Stats& stat = currentSnap_->statData;
1438	>	stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1439		*/
1440	+	//may output a header for the log file here
1441		getStatus();
1442		}
1443
#	Line 1081 | Line 1456 | namespace OpenMD {
1456		// all processors have the same number of bins, and STL vectors pack their
1457		// arrays, so in theory, this should be safe:
1458
1459	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueHist_[0],
1460	<	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1461	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[0],
1462	<	rnemdLogWidth_, MPI::INT, MPI::SUM);
1459	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &mHist_[0],
1460	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1461	>	if (outputTemp_) {
1462	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempHist_[0],
1463	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1464	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempCount_[0],
1465	>	rnemdLogWidth_, MPI::INT, MPI::SUM);
1466	>	}
1467	>	if (outputVx_) {
1468	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pxzHist_[0],
1469	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1470	>	//MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &vxzCount_[0],
1471	>	//                        rnemdLogWidth_, MPI::INT, MPI::SUM);
1472	>	}
1473	>	if (outputVy_) {
1474	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pyzHist_[0],
1475	>	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1476	>	//MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &vyzCount_[0],
1477	>	//                        rnemdLogWidth_, MPI::INT, MPI::SUM);
1478	>	}
1479		if (output3DTemp_) {
1480		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0],
1481		rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
#	Line 1095 | Line 1486 | namespace OpenMD {
1486		MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0],
1487		rnemdLogWidth_, MPI::INT, MPI::SUM);
1488		}
1489	+	if (outputRotTemp_) {
1490	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempHist_[0],
1491	+	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1492	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempCount_[0],
1493	+	rnemdLogWidth_, MPI::INT, MPI::SUM);
1494	+	}
1495	+
1496		// If we're the root node, should we print out the results
1497		int worldRank = MPI::COMM_WORLD.Get_rank();
1498		if (worldRank == 0) {
1499		#endif
1500	<	rnemdLog_ << time;
1501	<	for (j = 0; j < rnemdLogWidth_; j++) {
1502	<	rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j];
1500	>
1501	>	if (outputTemp_) {
1502	>	tempLog_ << time;
1503	>	for (j = 0; j < rnemdLogWidth_; j++) {
1504	>	tempLog_ << "\t" << tempHist_[j] / (RealType)tempCount_[j];
1505	>	}
1506	>	tempLog_ << endl;
1507		}
1508	<	rnemdLog_ << "\n";
1508	>	if (outputVx_) {
1509	>	vxzLog_ << time;
1510	>	for (j = 0; j < rnemdLogWidth_; j++) {
1511	>	vxzLog_ << "\t" << pxzHist_[j] / mHist_[j];
1512	>	}
1513	>	vxzLog_ << endl;
1514	>	}
1515	>	if (outputVy_) {
1516	>	vyzLog_ << time;
1517	>	for (j = 0; j < rnemdLogWidth_; j++) {
1518	>	vyzLog_ << "\t" << pyzHist_[j] / mHist_[j];
1519	>	}
1520	>	vyzLog_ << endl;
1521	>	}
1522	>
1523		if (output3DTemp_) {
1524	<	xTempLog_ << time;
1524	>	RealType temp;
1525	>	xTempLog_ << time;
1526		for (j = 0; j < rnemdLogWidth_; j++) {
1527	<	xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j];
1527	>	if (outputVx_)
1528	>	xTempHist_[j] -= pxzHist_[j] * pxzHist_[j] / mHist_[j];
1529	>	temp = xTempHist_[j] / (RealType)xyzTempCount_[j]
1530	>	/ PhysicalConstants::energyConvert / PhysicalConstants::kb;
1531	>	xTempLog_ << "\t" << temp;
1532		}
1533	<	xTempLog_ << "\n";
1533	>	xTempLog_ << endl;
1534		yTempLog_ << time;
1535		for (j = 0; j < rnemdLogWidth_; j++) {
1536		yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1537		}
1538	<	yTempLog_ << "\n";
1538	>	yTempLog_ << endl;
1539		zTempLog_ << time;
1540		for (j = 0; j < rnemdLogWidth_; j++) {
1541		zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1542		}
1543	<	zTempLog_ << "\n";
1543	>	zTempLog_ << endl;
1544		}
1545	+	if (outputRotTemp_) {
1546	+	rotTempLog_ << time;
1547	+	for (j = 0; j < rnemdLogWidth_; j++) {
1548	+	rotTempLog_ << "\t" << rotTempHist_[j] / (RealType)rotTempCount_[j];
1549	+	}
1550	+	rotTempLog_ << endl;
1551	+	}
1552	+
1553		#ifdef IS_MPI
1554		}
1555		#endif
1556	+
1557		for (j = 0; j < rnemdLogWidth_; j++) {
1558	<	valueCount_[j] = 0;
1129	<	valueHist_[j] = 0.0;
1558	>	mHist_[j] = 0.0;
1559		}
1560	+	if (outputTemp_)
1561	+	for (j = 0; j < rnemdLogWidth_; j++) {
1562	+	tempCount_[j] = 0;
1563	+	tempHist_[j] = 0.0;
1564	+	}
1565	+	if (outputVx_)
1566	+	for (j = 0; j < rnemdLogWidth_; j++) {
1567	+	//pxzCount_[j] = 0;
1568	+	pxzHist_[j] = 0.0;
1569	+	}
1570	+	if (outputVy_)
1571	+	for (j = 0; j < rnemdLogWidth_; j++) {
1572	+	//pyzCount_[j] = 0;
1573	+	pyzHist_[j] = 0.0;
1574	+	}
1575	+
1576		if (output3DTemp_)
1577		for (j = 0; j < rnemdLogWidth_; j++) {
1578		xTempHist_[j] = 0.0;
#	Line 1135 | Line 1580 | namespace OpenMD {
1580		zTempHist_[j] = 0.0;
1581		xyzTempCount_[j] = 0;
1582		}
1583	+	if (outputRotTemp_)
1584	+	for (j = 0; j < rnemdLogWidth_; j++) {
1585	+	rotTempCount_[j] = 0;
1586	+	rotTempHist_[j] = 0.0;
1587	+	}
1588		}
1589		}
1590	+

Diff Legend

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