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 1339 by gezelter, Thu Apr 23 18:31:05 2009 UTC vs.
branches/development/src/integrators/RNEMD.cpp (file contents), Revision 1722 by gezelter, Thu May 24 14:23:40 2012 UTC

#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	+	* [4]  Vardeman & Gezelter, in progress (2009).
40		*/
41
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"
49		#include "primitives/StuntDouble.hpp"
50	<	#include "utils/OOPSEConstant.hpp"
50	>	#include "utils/PhysicalConstants.hpp"
51		#include "utils/Tuple.hpp"
52
53		#ifndef IS_MPI
#	Line 53 | Line 56
56		#include "math/ParallelRandNumGen.hpp"
57		#endif
58
59	<	/* Remove me after testing*/
57	<	/*
58	<	#include <cstdio>
59	<	#include <iostream>
60	<	*/
61	<	/*End remove me*/
59	>	#define HONKING_LARGE_VALUE 1.0e10
60
61	<	namespace oopse {
61	>	using namespace std;
62	>	namespace OpenMD {
63
64	<	RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
65	<
64	>	RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
65	>	usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
66	>
67	>	failTrialCount_ = 0;
68	>	failRootCount_ = 0;
69	>
70		int seedValue;
71		Globals * simParams = info->getSimParams();
72
73	<	stringToEnumMap_["Kinetic"] = rnemdKinetic;
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();
77	–
78	–	std::cerr << "calling  evaluator with " << rnemdObjectSelection_ << "\n";
88		evaluator_.loadScriptString(rnemdObjectSelection_);
89	<	std::cerr << "done\n";
89	>	seleMan_.setSelectionSet(evaluator_.evaluate());
90	>
91	>	// do some sanity checking
92	>
93	>	int selectionCount = seleMan_.getSelectionCount();
94	>	int nIntegrable = info->getNGlobalIntegrableObjects();
95	>
96	>	if (selectionCount > nIntegrable) {
97	>	sprintf(painCave.errMsg,
98	>	"RNEMD: The current RNEMD_objectSelection,\n"
99	>	"\t\t%s\n"
100	>	"\thas resulted in %d selected objects.  However,\n"
101	>	"\tthe total number of integrable objects in the system\n"
102	>	"\tis only %d.  This is almost certainly not what you want\n"
103	>	"\tto do.  A likely cause of this is forgetting the _RB_0\n"
104	>	"\tselector in the selection script!\n",
105	>	rnemdObjectSelection_.c_str(),
106	>	selectionCount, nIntegrable);
107	>	painCave.isFatal = 0;
108	>	painCave.severity = OPENMD_WARNING;
109	>	simError();
110	>	}
111
112	<	const std::string st = simParams->getRNEMD_swapType();
112	>	const string st = simParams->getRNEMD_exchangeType();
113
114	<	std::map<std::string, RNEMDTypeEnum>::iterator i;
114	>	map<string, RNEMDTypeEnum>::iterator i;
115		i = stringToEnumMap_.find(st);
116	<	rnemdType_  = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
116	>	rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second;
117	>	if (rnemdType_ == rnemdUnknown) {
118	>	sprintf(painCave.errMsg,
119	>	"RNEMD: The current RNEMD_exchangeType,\n"
120	>	"\t\t%s\n"
121	>	"\tis not one of the recognized exchange types.\n",
122	>	st.c_str());
123	>	painCave.isFatal = 1;
124	>	painCave.severity = OPENMD_ERROR;
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_outputXyzTemperature()) {
151	>	output3DTemp_ = simParams->getRNEMD_outputXyzTemperature();
152	>	}
153	>	outputRotTemp_ = false;
154	>	if (simParams->haveRNEMD_outputRotTemperature()) {
155	>	outputRotTemp_ = simParams->getRNEMD_outputRotTemperature();
156	>	}
157
158	<	set_RNEMD_swapTime(simParams->getRNEMD_swapTime());
158	>	#ifdef IS_MPI
159	>	if (worldRank == 0) {
160	>	#endif
161	>
162	>	//may have rnemdWriter separately
163	>	string rnemdFileName;
164	>
165	>	if (outputTemp_) {
166	>	rnemdFileName = "temperature.log";
167	>	tempLog_.open(rnemdFileName.c_str());
168	>	}
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	>
178	>	if (output3DTemp_) {
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	>	}
193	>	#endif
194	>
195	>	set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
196		set_RNEMD_nBins(simParams->getRNEMD_nBins());
197	<	exchangeSum_ = 0.0;
198	<	counter_ = 0; //added by shenyu
199	<	//profile_.open("profile", std::ios::out);
197	>	midBin_ = nBins_ / 2;
198	>	if (simParams->haveRNEMD_binShift()) {
199	>	if (simParams->getRNEMD_binShift()) {
200	>	zShift_ = 0.5 / (RealType)(nBins_);
201	>	} else {
202	>	zShift_ = 0.0;
203	>	}
204	>	} else {
205	>	zShift_ = 0.0;
206	>	}
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) {
214	>	cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
215	>	cerr << "Automaically set back to default.\n";
216	>	rnemdLogWidth_ = nBins_;
217	>	}*/
218	>	} else {
219	>	set_RNEMD_logWidth(nBins_);
220	>	}
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()) {
238	+	set_RNEMD_target_flux(simParams->getRNEMD_targetFlux());
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()) {
271		seedValue = simParams->getSeed();
#	Line 110 | Line 285 | namespace oopse {
285
286		RNEMD::~RNEMD() {
287		delete randNumGen_;
288	<	//profile_.close();
288	>
289	>	#ifdef IS_MPI
290	>	if (worldRank == 0) {
291	>	#endif
292	>
293	>	sprintf(painCave.errMsg,
294	>	"RNEMD: total failed trials: %d\n",
295	>	failTrialCount_);
296	>	painCave.isFatal = 0;
297	>	painCave.severity = OPENMD_INFO;
298	>	simError();
299	>
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_);
309	>	painCave.isFatal = 0;
310	>	painCave.severity = OPENMD_INFO;
311	>	simError();
312	>	}
313	>	if (output3DTemp_) {
314	>	xTempLog_.close();
315	>	yTempLog_.close();
316	>	zTempLog_.close();
317	>	}
318	>	if (outputRotTemp_) rotTempLog_.close();
319	>
320	>	#ifdef IS_MPI
321	>	}
322	>	#endif
323		}
324
325		void RNEMD::doSwap() {
117	–	//std::cerr << "in RNEMD!\n";
118	–	//std::cerr << "nBins = " << nBins_ << "\n";
119	–	int midBin = nBins_ / 2;
120	–	//std::cerr << "midBin = " << midBin << "\n";
121	–	//std::cerr << "swapTime = " << swapTime_ << "\n";
122	–	//std::cerr << "swapType = " << rnemdType_ << "\n";
123	–	//std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
326
327		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
328		Mat3x3d hmat = currentSnap_->getHmat();
329
128	–	//std::cerr << "hmat = " << hmat << "\n";
129	–
330		seleMan_.setSelectionSet(evaluator_.evaluate());
331
132	–	//std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
133	–
332		int selei;
333		StuntDouble* sd;
334		int idx;
#	Line 158 | Line 356 | namespace oopse {
356		// which bin is this stuntdouble in?
357		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
358
359	<	int binNo = int((nBins_-1) * (pos.z() + 0.5*hmat(2,2)) / hmat(2,2));
359	>	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
360
163	–	//std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
361
362		// if we're in bin 0 or the middleBin
363	<	if (binNo == 0 || binNo == midBin) {
363	>	if (binNo == 0 || binNo == midBin_) {
364
365		RealType mass = sd->getMass();
366		Vector3d vel = sd->getVel();
367		RealType value;
368
369		switch(rnemdType_) {
370	<	case rnemdKinetic :
370	>	case rnemdKineticSwap :
371
372	<	value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
373	<	vel[2]*vel[2]);
374	<
178	<	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);
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);
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	<	}
390	<	value = value * 0.5 / OOPSEConstant::energyConvert;
389	>	} //angular momenta exchange enabled
390	>	//energyConvert temporarily disabled
391	>	//make exchangeSum_ comparable between swap & scale
392	>	//value = value * 0.5 / PhysicalConstants::energyConvert;
393	>	value *= 0.5;
394		break;
395		case rnemdPx :
396		value = mass * vel[0];
#	Line 202 | Line 401 | namespace oopse {
401		case rnemdPz :
402		value = mass * vel[2];
403		break;
205	–	case rnemdUnknown :
404		default :
405		break;
406		}
#	Line 218 | Line 416 | namespace oopse {
416		min_sd = sd;
417		}
418		}
419	<	} else {
419	>	} else { //midBin_
420		if (!max_found) {
421		max_val = value;
422		max_sd = sd;
#	Line 232 | Line 430 | namespace oopse {
430		}
431		}
432		}
235	–	//std::cerr << "smallest value = " << min_val  << "\n";
236	–	//std::cerr << "largest value = " << max_val << "\n";
237	–
238	–	// missing:  swap information in parallel
433
434	+	#ifdef IS_MPI
435	+	int nProc, worldRank;
436	+
437	+	nProc = MPI::COMM_WORLD.Get_size();
438	+	worldRank = MPI::COMM_WORLD.Get_rank();
439	+
440	+	bool my_min_found = min_found;
441	+	bool my_max_found = max_found;
442	+
443	+	// Even if we didn't find a minimum, did someone else?
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	+	#endif
448	+
449		if (max_found && min_found) {
450	<	if (min_val< max_val) {
451	<	Vector3d min_vel = min_sd->getVel();
452	<	Vector3d max_vel = max_sd->getVel();
453	<	RealType temp_vel;
454	<	switch(rnemdType_) {
455	<	case rnemdKinetic :
456	<	min_sd->setVel(max_vel);
457	<	max_sd->setVel(min_vel);
458	<	if (min_sd->isDirectional() && max_sd->isDirectional()) {
250	<	Vector3d min_angMom = min_sd->getJ();
251	<	Vector3d max_angMom = max_sd->getJ();
252	<	min_sd->setJ(max_angMom);
253	<	max_sd->setJ(min_angMom);
254	<	}
255	<	break;
256	<	case rnemdPx :
257	<	temp_vel = min_vel.x();
258	<	min_vel.x() = max_vel.x();
259	<	max_vel.x() = temp_vel;
260	<	min_sd->setVel(min_vel);
261	<	max_sd->setVel(max_vel);
262	<	break;
263	<	case rnemdPy :
264	<	temp_vel = min_vel.y();
265	<	min_vel.y() = max_vel.y();
266	<	max_vel.y() = temp_vel;
267	<	min_sd->setVel(min_vel);
268	<	max_sd->setVel(max_vel);
269	<	break;
270	<	case rnemdPz :
271	<	temp_vel = min_vel.z();
272	<	min_vel.z() = max_vel.z();
273	<	max_vel.z() = temp_vel;
274	<	min_sd->setVel(min_vel);
275	<	max_sd->setVel(max_vel);
276	<	break;
277	<	case rnemdUnknown :
278	<	default :
279	<	break;
280	<	}
281	<	exchangeSum_ += max_val - min_val;
450	>
451	>	#ifdef IS_MPI
452	>	struct {
453	>	RealType val;
454	>	int rank;
455	>	} max_vals, min_vals;
456	>
457	>	if (my_min_found) {
458	>	min_vals.val = min_val;
459		} else {
460	<	std::cerr << "exchange NOT performed.\nmin_val > max_val.\n";
460	>	min_vals.val = HONKING_LARGE_VALUE;
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;
468	>
469	>	if (my_max_found) {
470	>	max_vals.val = max_val;
471	>	} else {
472	>	max_vals.val = -HONKING_LARGE_VALUE;
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;
480	>	#endif
481	>
482	>	if (min_val < max_val) {
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	>
490	>	Vector3d min_vel = min_sd->getVel();
491	>	Vector3d max_vel = max_sd->getVel();
492	>	RealType temp_vel;
493	>
494	>	switch(rnemdType_) {
495	>	case rnemdKineticSwap :
496	>	min_sd->setVel(max_vel);
497	>	max_sd->setVel(min_vel);
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	>	}//angular momenta exchange enabled
504	>	//assumes same rigid body identity
505	>	break;
506	>	case rnemdPx :
507	>	temp_vel = min_vel.x();
508	>	min_vel.x() = max_vel.x();
509	>	max_vel.x() = temp_vel;
510	>	min_sd->setVel(min_vel);
511	>	max_sd->setVel(max_vel);
512	>	break;
513	>	case rnemdPy :
514	>	temp_vel = min_vel.y();
515	>	min_vel.y() = max_vel.y();
516	>	max_vel.y() = temp_vel;
517	>	min_sd->setVel(min_vel);
518	>	max_sd->setVel(max_vel);
519	>	break;
520	>	case rnemdPz :
521	>	temp_vel = min_vel.z();
522	>	min_vel.z() = max_vel.z();
523	>	max_vel.z() = temp_vel;
524	>	min_sd->setVel(min_vel);
525	>	max_sd->setVel(max_vel);
526	>	break;
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) {
534	>	// I had the max, but not the minimum
535	>
536	>	Vector3d min_vel;
537	>	Vector3d max_vel = max_sd->getVel();
538	>	MPI::Status status;
539	>
540	>	// point-to-point swap of the velocity vector
541	>	MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
542	>	min_vals.rank, 0,
543	>	min_vel.getArrayPointer(), 3, MPI::REALTYPE,
544	>	min_vals.rank, 0, status);
545	>
546	>	switch(rnemdType_) {
547	>	case rnemdKineticSwap :
548	>	max_sd->setVel(min_vel);
549	>	//angular momenta exchange enabled
550	>	if (max_sd->isDirectional()) {
551	>	Vector3d min_angMom;
552	>	Vector3d max_angMom = max_sd->getJ();
553	>
554	>	// point-to-point swap of the angular momentum vector
555	>	MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
556	>	MPI::REALTYPE, min_vals.rank, 1,
557	>	min_angMom.getArrayPointer(), 3,
558	>	MPI::REALTYPE, min_vals.rank, 1,
559	>	status);
560	>
561	>	max_sd->setJ(min_angMom);
562	>	}
563	>	break;
564	>	case rnemdPx :
565	>	max_vel.x() = min_vel.x();
566	>	max_sd->setVel(max_vel);
567	>	break;
568	>	case rnemdPy :
569	>	max_vel.y() = min_vel.y();
570	>	max_sd->setVel(max_vel);
571	>	break;
572	>	case rnemdPz :
573	>	max_vel.z() = min_vel.z();
574	>	max_sd->setVel(max_vel);
575	>	break;
576	>	default :
577	>	break;
578	>	}
579	>	} else if (min_vals.rank == worldRank) {
580	>	// I had the minimum but not the maximum:
581	>
582	>	Vector3d max_vel;
583	>	Vector3d min_vel = min_sd->getVel();
584	>	MPI::Status status;
585	>
586	>	// point-to-point swap of the velocity vector
587	>	MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
588	>	max_vals.rank, 0,
589	>	max_vel.getArrayPointer(), 3, MPI::REALTYPE,
590	>	max_vals.rank, 0, status);
591	>
592	>	switch(rnemdType_) {
593	>	case rnemdKineticSwap :
594	>	min_sd->setVel(max_vel);
595	>	//angular momenta exchange enabled
596	>	if (min_sd->isDirectional()) {
597	>	Vector3d min_angMom = min_sd->getJ();
598	>	Vector3d max_angMom;
599	>
600	>	// point-to-point swap of the angular momentum vector
601	>	MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
602	>	MPI::REALTYPE, max_vals.rank, 1,
603	>	max_angMom.getArrayPointer(), 3,
604	>	MPI::REALTYPE, max_vals.rank, 1,
605	>	status);
606	>
607	>	min_sd->setJ(max_angMom);
608	>	}
609	>	break;
610	>	case rnemdPx :
611	>	min_vel.x() = max_vel.x();
612	>	min_sd->setVel(min_vel);
613	>	break;
614	>	case rnemdPy :
615	>	min_vel.y() = max_vel.y();
616	>	min_sd->setVel(min_vel);
617	>	break;
618	>	case rnemdPz :
619	>	min_vel.z() = max_vel.z();
620	>	min_sd->setVel(min_vel);
621	>	break;
622	>	default :
623	>	break;
624	>	}
625	>	}
626	>	#endif
627	>	exchangeSum_ += max_val - min_val;
628	>	} else {
629	>	sprintf(painCave.errMsg,
630	>	"RNEMD: exchange NOT performed because min_val > max_val\n");
631	>	painCave.isFatal = 0;
632	>	painCave.severity = OPENMD_INFO;
633	>	simError();
634	>	failTrialCount_++;
635	>	}
636	>	} else {
637	>	sprintf(painCave.errMsg,
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();
643	>	failTrialCount_++;
644	>	}
645	>
646	>	}
647	>
648	>	void RNEMD::doScale() {
649	>
650	>	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
651	>	Mat3x3d hmat = currentSnap_->getHmat();
652	>
653	>	seleMan_.setSelectionSet(evaluator_.evaluate());
654	>
655	>	int selei;
656	>	StuntDouble* sd;
657	>	int idx;
658	>
659	>	vector<StuntDouble*> hotBin, coldBin;
660	>
661	>	RealType Phx = 0.0;
662	>	RealType Phy = 0.0;
663	>	RealType Phz = 0.0;
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)) {
678	>
679	>	idx = sd->getLocalIndex();
680	>
681	>	Vector3d pos = sd->getPos();
682	>
683	>	// wrap the stuntdouble's position back into the box:
684	>
685	>	if (usePeriodicBoundaryConditions_)
686	>	currentSnap_->wrapVector(pos);
687	>
688	>	// which bin is this stuntdouble in?
689	>	// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
690	>
691	>	int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_;
692	>
693	>	// if we're in bin 0 or the middleBin
694	>	if (binNo == 0 || binNo == midBin_) {
695	>
696	>	RealType mass = sd->getMass();
697	>	Vector3d vel = sd->getVel();
698	>
699	>	if (binNo == 0) {
700	>	hotBin.push_back(sd);
701	>	Phx += mass * vel.x();
702	>	Phy += mass * vel.y();
703	>	Phz += mass * vel.z();
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();
727	>	Pcy += mass * vel.y();
728	>	Pcz += mass * vel.z();
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	>
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);
771	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
772	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
773	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
774	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
775	>
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	>	//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	>	if (rnemdType_ == rnemdKineticScaleVAM) {
798	>	c = 1.0 - targetFlux_ / (Kcx + Kcy + Kcz + Kcw);
799	>	} else {
800	>	c = 1.0 - targetFlux_ / Kcw;
801	>	}
802	>
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		} else {
872	<	std::cerr << "exchange NOT performed.\none of the two slabs empty.\n";
872	>	RealType a000, a110, c0, a001, a111, b01, b11, c1;
873	>	switch(rnemdType_) {
874	>	case rnemdKineticScale :
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	>	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	>	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;
936	>	}
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	>	// 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		}
1094	<	std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
1094	>	if (successfulScale != true) {
1095	>	sprintf(painCave.errMsg,
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::getStatus() {
292	<	//std::cerr << "in RNEMD!\n";
293	<	//std::cerr << "nBins = " << nBins_ << "\n";
294	<	int midBin = nBins_ / 2;
295	<	//std::cerr << "midBin = " << midBin << "\n";
296	<	//std::cerr << "swapTime = " << swapTime_ << "\n";
297	<	//std::cerr << "exchangeSum = " << exchangeSum_ << "\n";
298	<	//std::cerr << "swapType = " << rnemdType_ << "\n";
299	<	//std::cerr << "selection = " << rnemdObjectSelection_ << "\n";
1104	>	void RNEMD::doShiftScale() {
1105
1106		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1107		Mat3x3d hmat = currentSnap_->getHmat();
1108
304	–	//std::cerr << "hmat = " << hmat << "\n";
305	–
1109		seleMan_.setSelectionSet(evaluator_.evaluate());
1110
1111	<	//std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n";
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 :
1289	+	doScale();
1290	+	break;
1291	+	case rnemdKineticSwap :
1292	+	case rnemdPx :
1293	+	case rnemdPy :
1294	+	case rnemdPz :
1295	+	doSwap();
1296	+	break;
1297	+	case rnemdShiftScaleV :
1298	+	case rnemdShiftScaleVAM :
1299	+	doShiftScale();
1300	+	break;
1301	+	case rnemdUnknown :
1302	+	default :
1303	+	break;
1304	+	}
1305	+	}
1306	+
1307	+	void RNEMD::collectData() {
1308	+
1309	+	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1310	+	Mat3x3d hmat = currentSnap_->getHmat();
1311	+
1312	+	seleMan_.setSelectionSet(evaluator_.evaluate());
1313	+
1314		int selei;
1315		StuntDouble* sd;
1316		int idx;
313	–	/*
314	–	RealType min_val;
315	–	bool min_found = false;
316	–	StuntDouble* min_sd;
1317
1318	<	RealType max_val;
1319	<	bool max_found = false;
1320	<	StuntDouble* max_sd;
1318	>	// alternative approach, track all molecules instead of only those
1319	>	// selected for scaling/swapping:
1320	>	/*
1321	>	SimInfo::MoleculeIterator miter;
1322	>	vector<StuntDouble*>::iterator iiter;
1323	>	Molecule* mol;
1324	>	StuntDouble* integrableObject;
1325	>	for (mol = info_->beginMolecule(miter); mol != NULL;
1326	>	mol = info_->nextMolecule(miter))
1327	>	integrableObject is essentially sd
1328	>	for (integrableObject = mol->beginIntegrableObject(iiter);
1329	>	integrableObject != NULL;
1330	>	integrableObject = mol->nextIntegrableObject(iiter))
1331		*/
322	–	std::vector<RealType> valueHist;  // keeps track of what's being averaged
323	–	std::vector<int> valueCount; // keeps track of the number of degrees of
324	–	// freedom being averaged
325	–	valueHist.resize(nBins_);
326	–	valueCount.resize(nBins_);
327	–	//do they initialize themselves to zero automatically?
1332		for (sd = seleMan_.beginSelected(selei); sd != NULL;
1333		sd = seleMan_.nextSelected(selei)) {
1334
#	Line 340 | Line 1344 | namespace oopse {
1344		// which bin is this stuntdouble in?
1345		// wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1346
1347	<	int binNo = int((nBins_-1) * (pos.z()+0.5*hmat(2,2)) / hmat(2,2));
1348	<
1349	<	//std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n";
1350	<
1347	>	int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) %
1348	>	rnemdLogWidth_;
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();
349	–	//std::cerr << "mass = " << mass << " vel = " << vel << "\n";
1358		RealType value;
1359	+	//RealType xVal, yVal, zVal;
1360
1361	<	switch(rnemdType_) {
1362	<	case rnemdKinetic :
1363	<
355	<	value = mass * (vel[0]*vel[0] + vel[1]*vel[1] +
356	<	vel[2]*vel[2]);
357	<
358	<	valueCount[binNo] += 3;
359	<
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();
363	–
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	<
1374	<	valueCount[binNo] +=2;
1375	<
1376	<	} else {
1377	<	value += angMom[0]*angMom[0]/I(0, 0)
1378	<	+ angMom[1]*angMom[1]/I(1, 1)
376	<	+ angMom[2]*angMom[2]/I(2, 2);
377	<	valueCount[binNo] +=3;
378	<
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	>	tempCount_[binNo] +=3;
1379		}
1380		}
1381	<	//std::cerr <<"this value = " << value << "\n";
1382	<	value *= 0.5 / OOPSEConstant::energyConvert;  // get it in kcal / mol
1383	<	value *= 2.0 / OOPSEConstant::kb;             // convert to temperature
1384	<	//std::cerr <<"this value = " << value << "\n";
1385	<	break;
386	<	case rnemdPx :
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 :
1387	>	//vxzCount_[binNo]++;
1388	>	pxzHist_[binNo] += value;
1389	>	}
1390	>	if (outputVy_) {
1391		value = mass * vel[1];
1392	<	valueCount[binNo]++;
1393	<	break;
394	<	case rnemdPz :
395	<	value = mass * vel[2];
396	<	valueCount[binNo]++;
397	<	break;
398	<	case rnemdUnknown :
399	<	default :
400	<	break;
1392	>	//vyzCount_[binNo]++;
1393	>	pyzHist_[binNo] += value;
1394		}
1395	<	//std::cerr << "bin = " << binNo << " value = " << value ;
1396	<	valueHist[binNo] += value;
1397	<	//std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n";
1395	>
1396	>	if (output3DTemp_) {
1397	>	value = mass * vel.x() * vel.x();
1398	>	xTempHist_[binNo] += value;
1399	>	value = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert
1400	>	/ PhysicalConstants::kb;
1401	>	yTempHist_[binNo] += value;
1402	>	value = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert
1403	>	/ PhysicalConstants::kb;
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		}
406	–
407	–	std::cout << counter_++;
408	–	for (int j = 0; j < nBins_; j++)
409	–	std::cout << "\t" << valueHist[j] / (RealType)valueCount[j];
410	–	std::cout << "\n";
1431		}
1432	+
1433	+	void RNEMD::getStarted() {
1434	+	collectData();
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	+
1444	+	void RNEMD::getStatus() {
1445	+
1446	+	Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1447	+	Stats& stat = currentSnap_->statData;
1448	+	RealType time = currentSnap_->getTime();
1449	+
1450	+	stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1451	+	//or to be more meaningful, define another item as exchangeSum_ / time
1452	+	int j;
1453	+
1454	+	#ifdef IS_MPI
1455	+
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, &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);
1482	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0],
1483	+	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1484	+	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0],
1485	+	rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
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	+
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	+	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	+	RealType temp;
1525	+	xTempLog_ << time;
1526	+	for (j = 0; j < rnemdLogWidth_; 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_ << endl;
1534	+	yTempLog_ << time;
1535	+	for (j = 0; j < rnemdLogWidth_; j++) {
1536	+	yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j];
1537	+	}
1538	+	yTempLog_ << endl;
1539	+	zTempLog_ << time;
1540	+	for (j = 0; j < rnemdLogWidth_; j++) {
1541	+	zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j];
1542	+	}
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	+	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;
1579	+	yTempHist_[j] = 0.0;
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	+

Comparing:
trunk/src/integrators/RNEMD.cpp (property svn:keywords), Revision 1339 by gezelter, Thu Apr 23 18:31:05 2009 UTC vs.
branches/development/src/integrators/RNEMD.cpp (property svn:keywords), Revision 1722 by gezelter, Thu May 24 14:23:40 2012 UTC

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

Diff Legend

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