ViewVC Help

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

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

Comparing:
trunk/src/integrators/RNEMD.cpp (file contents), Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/integrators/RNEMD.cpp (file contents), Revision 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC

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

Comparing:
trunk/src/integrators/RNEMD.cpp (property svn:keywords), Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/integrators/RNEMD.cpp (property svn:keywords), Revision 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC

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

Diff Legend

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