[ViewVC] Diff of: OpenMD/trunk/src/rnemd/RNEMD.cpp

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

#	Line 68 \| Line 68 \| namespace OpenMD {
68		using namespace std;
69		namespace OpenMD {
70
71	<	RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
72	<	evaluatorA_(info), seleManA_(info),
73	<	commonA_(info), evaluatorB_(info),
74	<	seleManB_(info), commonB_(info),
75	<	hasData_(false), hasDividingArea_(false),
76	<	usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
71	>	RNEMD::RNEMD(SimInfo* info) : info_(info),
72	>	evaluator_(info_), seleMan_(info_),
73	>	evaluatorA_(info_), seleManA_(info_),
74	>	evaluatorB_(info_), seleManB_(info_),
75	>	commonA_(info_), commonB_(info_),
76	>	usePeriodicBoundaryConditions_(info_->getSimParams()->getUsePeriodicBoundaryConditions()),
77	>	hasDividingArea_(false),
78	>	hasData_(false) {
79
80		trialCount_ = 0;
81		failTrialCount_ = 0;
#	Line 291 \| Line 293 \| namespace OpenMD {
293		kineticFlux_ = 0.0;
294		}
295		if (hasMomentumFluxVector) {
296	<	momentumFluxVector_ = rnemdParams->getMomentumFluxVector();
296	>	std::vector<RealType> mf = rnemdParams->getMomentumFluxVector();
297	>	if (mf.size() != 3) {
298	>	sprintf(painCave.errMsg,
299	>	"RNEMD: Incorrect number of parameters specified for momentumFluxVector.\n"
300	>	"\tthere should be 3 parameters, but %lu were specified.\n",
301	>	mf.size());
302	>	painCave.isFatal = 1;
303	>	simError();
304	>	}
305	>	momentumFluxVector_.x() = mf[0];
306	>	momentumFluxVector_.y() = mf[1];
307	>	momentumFluxVector_.z() = mf[2];
308		} else {
309		momentumFluxVector_ = V3Zero;
310		if (hasMomentumFlux) {
#	Line 315 \| Line 328 \| namespace OpenMD {
328		default:
329		break;
330		}
318	–	}
319	–	if (hasAngularMomentumFluxVector) {
320	–	angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
321	–	} else {
322	–	angularMomentumFluxVector_ = V3Zero;
323	–	if (hasAngularMomentumFlux) {
324	–	RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
325	–	switch (rnemdFluxType_) {
326	–	case rnemdLx:
327	–	angularMomentumFluxVector_.x() = angularMomentumFlux;
328	–	break;
329	–	case rnemdLy:
330	–	angularMomentumFluxVector_.y() = angularMomentumFlux;
331	–	break;
332	–	case rnemdLz:
333	–	angularMomentumFluxVector_.z() = angularMomentumFlux;
334	–	break;
335	–	case rnemdKeLx:
336	–	angularMomentumFluxVector_.x() = angularMomentumFlux;
337	–	break;
338	–	case rnemdKeLy:
339	–	angularMomentumFluxVector_.y() = angularMomentumFlux;
340	–	break;
341	–	case rnemdKeLz:
342	–	angularMomentumFluxVector_.z() = angularMomentumFlux;
343	–	break;
344	–	default:
345	–	break;
346	–	}
347	–	}
331		}
332	<
333	<	if (hasCoordinateOrigin) {
334	<	coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
335	<	} else {
336	<	coordinateOrigin_ = V3Zero;
337	<	}
338	<
339	<	// do some sanity checking
340	<
341	<	int selectionCount = seleMan_.getSelectionCount();
342	<
343	<	int nIntegrable = info->getNGlobalIntegrableObjects();
344	<
345	<	if (selectionCount > nIntegrable) {
346	<	sprintf(painCave.errMsg,
347	<	"RNEMD: The current objectSelection,\n"
348	<	"\t\t%s\n"
349	<	"\thas resulted in %d selected objects. However,\n"
367	<	"\tthe total number of integrable objects in the system\n"
368	<	"\tis only %d. This is almost certainly not what you want\n"
369	<	"\tto do. A likely cause of this is forgetting the _RB_0\n"
370	<	"\tselector in the selection script!\n",
371	<	rnemdObjectSelection_.c_str(),
372	<	selectionCount, nIntegrable);
373	<	painCave.isFatal = 0;
374	<	painCave.severity = OPENMD_WARNING;
375	<	simError();
376	<	}
377	<
378	<	areaAccumulator_ = new Accumulator();
379	<
380	<	nBins_ = rnemdParams->getOutputBins();
381	<	binWidth_ = rnemdParams->getOutputBinWidth();
382	<
383	<	data_.resize(RNEMD::ENDINDEX);
384	<	OutputData z;
385	<	z.units = "Angstroms";
386	<	z.title = "Z";
387	<	z.dataType = "RealType";
388	<	z.accumulator.reserve(nBins_);
389	<	for (int i = 0; i < nBins_; i++)
390	<	z.accumulator.push_back( new Accumulator() );
391	<	data_[Z] = z;
392	<	outputMap_["Z"] = Z;
393	<
394	<	OutputData r;
395	<	r.units = "Angstroms";
396	<	r.title = "R";
397	<	r.dataType = "RealType";
398	<	r.accumulator.reserve(nBins_);
399	<	for (int i = 0; i < nBins_; i++)
400	<	r.accumulator.push_back( new Accumulator() );
401	<	data_[R] = r;
402	<	outputMap_["R"] = R;
403	<
404	<	OutputData temperature;
405	<	temperature.units = "K";
406	<	temperature.title = "Temperature";
407	<	temperature.dataType = "RealType";
408	<	temperature.accumulator.reserve(nBins_);
409	<	for (int i = 0; i < nBins_; i++)
410	<	temperature.accumulator.push_back( new Accumulator() );
411	<	data_[TEMPERATURE] = temperature;
412	<	outputMap_["TEMPERATURE"] = TEMPERATURE;
413	<
414	<	OutputData velocity;
415	<	velocity.units = "angstroms/fs";
416	<	velocity.title = "Velocity";
417	<	velocity.dataType = "Vector3d";
418	<	velocity.accumulator.reserve(nBins_);
419	<	for (int i = 0; i < nBins_; i++)
420	<	velocity.accumulator.push_back( new VectorAccumulator() );
421	<	data_[VELOCITY] = velocity;
422	<	outputMap_["VELOCITY"] = VELOCITY;
423	<
424	<	OutputData angularVelocity;
425	<	angularVelocity.units = "angstroms^2/fs";
426	<	angularVelocity.title = "AngularVelocity";
427	<	angularVelocity.dataType = "Vector3d";
428	<	angularVelocity.accumulator.reserve(nBins_);
429	<	for (int i = 0; i < nBins_; i++)
430	<	angularVelocity.accumulator.push_back( new VectorAccumulator() );
431	<	data_[ANGULARVELOCITY] = angularVelocity;
432	<	outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
433	<
434	<	OutputData density;
435	<	density.units = "g cm^-3";
436	<	density.title = "Density";
437	<	density.dataType = "RealType";
438	<	density.accumulator.reserve(nBins_);
439	<	for (int i = 0; i < nBins_; i++)
440	<	density.accumulator.push_back( new Accumulator() );
441	<	data_[DENSITY] = density;
442	<	outputMap_["DENSITY"] = DENSITY;
443	<
444	<	if (hasOutputFields) {
445	<	parseOutputFileFormat(rnemdParams->getOutputFields());
446	<	} else {
447	<	if (usePeriodicBoundaryConditions_)
448	<	outputMask_.set(Z);
449	<	else
450	<	outputMask_.set(R);
332	>	}
333	>	if (hasAngularMomentumFluxVector) {
334	>	std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
335	>	if (amf.size() != 3) {
336	>	sprintf(painCave.errMsg,
337	>	"RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
338	>	"\tthere should be 3 parameters, but %lu were specified.\n",
339	>	amf.size());
340	>	painCave.isFatal = 1;
341	>	simError();
342	>	}
343	>	angularMomentumFluxVector_.x() = amf[0];
344	>	angularMomentumFluxVector_.y() = amf[1];
345	>	angularMomentumFluxVector_.z() = amf[2];
346	>	} else {
347	>	angularMomentumFluxVector_ = V3Zero;
348	>	if (hasAngularMomentumFlux) {
349	>	RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
350		switch (rnemdFluxType_) {
452	–	case rnemdKE:
453	–	case rnemdRotKE:
454	–	case rnemdFullKE:
455	–	outputMask_.set(TEMPERATURE);
456	–	break;
457	–	case rnemdPx:
458	–	case rnemdPy:
459	–	outputMask_.set(VELOCITY);
460	–	break;
461	–	case rnemdPz:
462	–	case rnemdPvector:
463	–	outputMask_.set(VELOCITY);
464	–	outputMask_.set(DENSITY);
465	–	break;
351		case rnemdLx:
352	+	angularMomentumFluxVector_.x() = angularMomentumFlux;
353	+	break;
354		case rnemdLy:
355	+	angularMomentumFluxVector_.y() = angularMomentumFlux;
356	+	break;
357		case rnemdLz:
358	<	case rnemdLvector:
470	<	outputMask_.set(ANGULARVELOCITY);
358	>	angularMomentumFluxVector_.z() = angularMomentumFlux;
359		break;
360		case rnemdKeLx:
361	+	angularMomentumFluxVector_.x() = angularMomentumFlux;
362	+	break;
363		case rnemdKeLy:
364	+	angularMomentumFluxVector_.y() = angularMomentumFlux;
365	+	break;
366		case rnemdKeLz:
367	<	case rnemdKeLvector:
476	<	outputMask_.set(TEMPERATURE);
477	<	outputMask_.set(ANGULARVELOCITY);
367	>	angularMomentumFluxVector_.z() = angularMomentumFlux;
368		break;
479	–	case rnemdKePx:
480	–	case rnemdKePy:
481	–	outputMask_.set(TEMPERATURE);
482	–	outputMask_.set(VELOCITY);
483	–	break;
484	–	case rnemdKePvector:
485	–	outputMask_.set(TEMPERATURE);
486	–	outputMask_.set(VELOCITY);
487	–	outputMask_.set(DENSITY);
488	–	break;
369		default:
370		break;
371		}
372	<	}
373	<
494	<	if (hasOutputFileName) {
495	<	rnemdFileName_ = rnemdParams->getOutputFileName();
496	<	} else {
497	<	rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
498	<	}
372	>	}
373	>	}
374
375	<	exchangeTime_ = rnemdParams->getExchangeTime();
376	<
377	<	Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
378	<	// total exchange sums are zeroed out at the beginning:
379	<
380	<	kineticExchange_ = 0.0;
381	<	momentumExchange_ = V3Zero;
382	<	angularMomentumExchange_ = V3Zero;
383	<
384	<	std::ostringstream selectionAstream;
385	<	std::ostringstream selectionBstream;
375	>	if (hasCoordinateOrigin) {
376	>	std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
377	>	if (co.size() != 3) {
378	>	sprintf(painCave.errMsg,
379	>	"RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
380	>	"\tthere should be 3 parameters, but %lu were specified.\n",
381	>	co.size());
382	>	painCave.isFatal = 1;
383	>	simError();
384	>	}
385	>	coordinateOrigin_.x() = co[0];
386	>	coordinateOrigin_.y() = co[1];
387	>	coordinateOrigin_.z() = co[2];
388	>	} else {
389	>	coordinateOrigin_ = V3Zero;
390	>	}
391
392	<	if (hasSelectionA_) {
393	<	selectionA_ = rnemdParams->getSelectionA();
394	<	} else {
395	<	if (usePeriodicBoundaryConditions_) {
396	<	Mat3x3d hmat = currentSnap_->getHmat();
397	<
398	<	if (hasSlabWidth)
399	<	slabWidth_ = rnemdParams->getSlabWidth();
400	<	else
401	<	slabWidth_ = hmat(2,2) / 10.0;
402	<
403	<	if (hasSlabACenter)
404	<	slabACenter_ = rnemdParams->getSlabACenter();
405	<	else
406	<	slabACenter_ = 0.0;
407	<
408	<	selectionAstream << "select wrappedz > "
409	<	<< slabACenter_ - 0.5*slabWidth_
410	<	<< " && wrappedz < "
411	<	<< slabACenter_ + 0.5*slabWidth_;
412	<	selectionA_ = selectionAstream.str();
413	<	} else {
414	<	if (hasSphereARadius)
415	<	sphereARadius_ = rnemdParams->getSphereARadius();
416	<	else {
417	<	// use an initial guess to the size of the inner slab to be 1/10 the
418	<	// radius of an approximately spherical hull:
419	<	Thermo thermo(info);
420	<	RealType hVol = thermo.getHullVolume();
421	<	sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
422	<	}
423	<	selectionAstream << "select r < " << sphereARadius_;
424	<	selectionA_ = selectionAstream.str();
425	<	}
392	>	// do some sanity checking
393	>
394	>	int selectionCount = seleMan_.getSelectionCount();
395	>	int nIntegrable = info->getNGlobalIntegrableObjects();
396	>	if (selectionCount > nIntegrable) {
397	>	sprintf(painCave.errMsg,
398	>	"RNEMD: The current objectSelection,\n"
399	>	"\t\t%s\n"
400	>	"\thas resulted in %d selected objects. However,\n"
401	>	"\tthe total number of integrable objects in the system\n"
402	>	"\tis only %d. This is almost certainly not what you want\n"
403	>	"\tto do. A likely cause of this is forgetting the _RB_0\n"
404	>	"\tselector in the selection script!\n",
405	>	rnemdObjectSelection_.c_str(),
406	>	selectionCount, nIntegrable);
407	>	painCave.isFatal = 0;
408	>	painCave.severity = OPENMD_WARNING;
409	>	simError();
410	>	}
411	>
412	>	areaAccumulator_ = new Accumulator();
413	>
414	>	nBins_ = rnemdParams->getOutputBins();
415	>	binWidth_ = rnemdParams->getOutputBinWidth();
416	>
417	>	data_.resize(RNEMD::ENDINDEX);
418	>	OutputData z;
419	>	z.units = "Angstroms";
420	>	z.title = "Z";
421	>	z.dataType = "RealType";
422	>	z.accumulator.reserve(nBins_);
423	>	for (int i = 0; i < nBins_; i++)
424	>	z.accumulator.push_back( new Accumulator() );
425	>	data_[Z] = z;
426	>	outputMap_["Z"] = Z;
427	>
428	>	OutputData r;
429	>	r.units = "Angstroms";
430	>	r.title = "R";
431	>	r.dataType = "RealType";
432	>	r.accumulator.reserve(nBins_);
433	>	for (int i = 0; i < nBins_; i++)
434	>	r.accumulator.push_back( new Accumulator() );
435	>	data_[R] = r;
436	>	outputMap_["R"] = R;
437	>
438	>	OutputData temperature;
439	>	temperature.units = "K";
440	>	temperature.title = "Temperature";
441	>	temperature.dataType = "RealType";
442	>	temperature.accumulator.reserve(nBins_);
443	>	for (int i = 0; i < nBins_; i++)
444	>	temperature.accumulator.push_back( new Accumulator() );
445	>	data_[TEMPERATURE] = temperature;
446	>	outputMap_["TEMPERATURE"] = TEMPERATURE;
447	>
448	>	OutputData velocity;
449	>	velocity.units = "angstroms/fs";
450	>	velocity.title = "Velocity";
451	>	velocity.dataType = "Vector3d";
452	>	velocity.accumulator.reserve(nBins_);
453	>	for (int i = 0; i < nBins_; i++)
454	>	velocity.accumulator.push_back( new VectorAccumulator() );
455	>	data_[VELOCITY] = velocity;
456	>	outputMap_["VELOCITY"] = VELOCITY;
457	>
458	>	OutputData angularVelocity;
459	>	angularVelocity.units = "angstroms^2/fs";
460	>	angularVelocity.title = "AngularVelocity";
461	>	angularVelocity.dataType = "Vector3d";
462	>	angularVelocity.accumulator.reserve(nBins_);
463	>	for (int i = 0; i < nBins_; i++)
464	>	angularVelocity.accumulator.push_back( new VectorAccumulator() );
465	>	data_[ANGULARVELOCITY] = angularVelocity;
466	>	outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
467	>
468	>	OutputData density;
469	>	density.units = "g cm^-3";
470	>	density.title = "Density";
471	>	density.dataType = "RealType";
472	>	density.accumulator.reserve(nBins_);
473	>	for (int i = 0; i < nBins_; i++)
474	>	density.accumulator.push_back( new Accumulator() );
475	>	data_[DENSITY] = density;
476	>	outputMap_["DENSITY"] = DENSITY;
477	>
478	>	if (hasOutputFields) {
479	>	parseOutputFileFormat(rnemdParams->getOutputFields());
480	>	} else {
481	>	if (usePeriodicBoundaryConditions_)
482	>	outputMask_.set(Z);
483	>	else
484	>	outputMask_.set(R);
485	>	switch (rnemdFluxType_) {
486	>	case rnemdKE:
487	>	case rnemdRotKE:
488	>	case rnemdFullKE:
489	>	outputMask_.set(TEMPERATURE);
490	>	break;
491	>	case rnemdPx:
492	>	case rnemdPy:
493	>	outputMask_.set(VELOCITY);
494	>	break;
495	>	case rnemdPz:
496	>	case rnemdPvector:
497	>	outputMask_.set(VELOCITY);
498	>	outputMask_.set(DENSITY);
499	>	break;
500	>	case rnemdLx:
501	>	case rnemdLy:
502	>	case rnemdLz:
503	>	case rnemdLvector:
504	>	outputMask_.set(ANGULARVELOCITY);
505	>	break;
506	>	case rnemdKeLx:
507	>	case rnemdKeLy:
508	>	case rnemdKeLz:
509	>	case rnemdKeLvector:
510	>	outputMask_.set(TEMPERATURE);
511	>	outputMask_.set(ANGULARVELOCITY);
512	>	break;
513	>	case rnemdKePx:
514	>	case rnemdKePy:
515	>	outputMask_.set(TEMPERATURE);
516	>	outputMask_.set(VELOCITY);
517	>	break;
518	>	case rnemdKePvector:
519	>	outputMask_.set(TEMPERATURE);
520	>	outputMask_.set(VELOCITY);
521	>	outputMask_.set(DENSITY);
522	>	break;
523	>	default:
524	>	break;
525		}
526	+	}
527
528	<	if (hasSelectionB_) {
529	<	selectionB_ = rnemdParams->getSelectionB();
530	<
531	<	} else {
532	<	if (usePeriodicBoundaryConditions_) {
533	<	Mat3x3d hmat = currentSnap_->getHmat();
528	>	if (hasOutputFileName) {
529	>	rnemdFileName_ = rnemdParams->getOutputFileName();
530	>	} else {
531	>	rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
532	>	}
533	>
534	>	exchangeTime_ = rnemdParams->getExchangeTime();
535	>
536	>	Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
537	>	// total exchange sums are zeroed out at the beginning:
538	>
539	>	kineticExchange_ = 0.0;
540	>	momentumExchange_ = V3Zero;
541	>	angularMomentumExchange_ = V3Zero;
542	>
543	>	std::ostringstream selectionAstream;
544	>	std::ostringstream selectionBstream;
545	>
546	>	if (hasSelectionA_) {
547	>	selectionA_ = rnemdParams->getSelectionA();
548	>	} else {
549	>	if (usePeriodicBoundaryConditions_) {
550	>	Mat3x3d hmat = currentSnap_->getHmat();
551
552	<	if (hasSlabWidth)
553	<	slabWidth_ = rnemdParams->getSlabWidth();
554	<	else
555	<	slabWidth_ = hmat(2,2) / 10.0;
552	>	if (hasSlabWidth)
553	>	slabWidth_ = rnemdParams->getSlabWidth();
554	>	else
555	>	slabWidth_ = hmat(2,2) / 10.0;
556
557	<	if (hasSlabBCenter)
558	<	slabBCenter_ = rnemdParams->getSlabBCenter();
559	<	else
560	<	slabBCenter_ = hmat(2,2) / 2.0;
557	>	if (hasSlabACenter)
558	>	slabACenter_ = rnemdParams->getSlabACenter();
559	>	else
560	>	slabACenter_ = 0.0;
561
562	<	selectionBstream << "select wrappedz > "
563	<	<< slabBCenter_ - 0.5*slabWidth_
564	<	<< " && wrappedz < "
565	<	<< slabBCenter_ + 0.5*slabWidth_;
562	>	selectionAstream << "select wrappedz > "
563	>	<< slabACenter_ - 0.5*slabWidth_
564	>	<< " && wrappedz < "
565	>	<< slabACenter_ + 0.5*slabWidth_;
566	>	selectionA_ = selectionAstream.str();
567	>	} else {
568	>	if (hasSphereARadius)
569	>	sphereARadius_ = rnemdParams->getSphereARadius();
570	>	else {
571	>	// use an initial guess to the size of the inner slab to be 1/10 the
572	>	// radius of an approximately spherical hull:
573	>	Thermo thermo(info);
574	>	RealType hVol = thermo.getHullVolume();
575	>	sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
576	>	}
577	>	selectionAstream << "select r < " << sphereARadius_;
578	>	selectionA_ = selectionAstream.str();
579	>	}
580	>	}
581	>
582	>	if (hasSelectionB_) {
583	>	selectionB_ = rnemdParams->getSelectionB();
584	>
585	>	} else {
586	>	if (usePeriodicBoundaryConditions_) {
587	>	Mat3x3d hmat = currentSnap_->getHmat();
588	>
589	>	if (hasSlabWidth)
590	>	slabWidth_ = rnemdParams->getSlabWidth();
591	>	else
592	>	slabWidth_ = hmat(2,2) / 10.0;
593	>
594	>	if (hasSlabBCenter)
595	>	slabBCenter_ = rnemdParams->getSlabBCenter();
596	>	else
597	>	slabBCenter_ = hmat(2,2) / 2.0;
598	>
599	>	selectionBstream << "select wrappedz > "
600	>	<< slabBCenter_ - 0.5*slabWidth_
601	>	<< " && wrappedz < "
602	>	<< slabBCenter_ + 0.5*slabWidth_;
603	>	selectionB_ = selectionBstream.str();
604	>	} else {
605	>	if (hasSphereBRadius_) {
606	>	sphereBRadius_ = rnemdParams->getSphereBRadius();
607	>	selectionBstream << "select r > " << sphereBRadius_;
608		selectionB_ = selectionBstream.str();
609		} else {
610	<	if (hasSphereBRadius_) {
611	<	sphereBRadius_ = rnemdParams->getSphereBRadius();
612	<	selectionBstream << "select r > " << sphereBRadius_;
574	<	selectionB_ = selectionBstream.str();
575	<	} else {
576	<	selectionB_ = "select hull";
577	<	BisHull_ = true;
578	<	hasSelectionB_ = true;
579	<	}
610	>	selectionB_ = "select hull";
611	>	BisHull_ = true;
612	>	hasSelectionB_ = true;
613		}
614		}
615		}
616	<
616	>
617	>
618		// object evaluator:
619		evaluator_.loadScriptString(rnemdObjectSelection_);
620		seleMan_.setSelectionSet(evaluator_.evaluate());
#	Line 622 \| Line 656 \| namespace OpenMD {
656
657		StuntDouble* sd;
658
659	<	RealType min_val;
660	<	bool min_found = false;
661	<	StuntDouble* min_sd;
659	>	RealType min_val(0.0);
660	>	int min_found = 0;
661	>	StuntDouble* min_sd = NULL;
662
663	<	RealType max_val;
664	<	bool max_found = false;
665	<	StuntDouble* max_sd;
663	>	RealType max_val(0.0);
664	>	int max_found = 0;
665	>	StuntDouble* max_sd = NULL;
666
667		for (sd = seleManA_.beginSelected(selei); sd != NULL;
668		sd = seleManA_.nextSelected(selei)) {
#	Line 642 \| Line 676 \| namespace OpenMD {
676
677		RealType mass = sd->getMass();
678		Vector3d vel = sd->getVel();
679	<	RealType value;
679	>	RealType value(0.0);
680
681		switch(rnemdFluxType_) {
682		case rnemdKE :
#	Line 682 \| Line 716 \| namespace OpenMD {
716		if (!max_found) {
717		max_val = value;
718		max_sd = sd;
719	<	max_found = true;
719	>	max_found = 1;
720		} else {
721		if (max_val < value) {
722		max_val = value;
#	Line 703 \| Line 737 \| namespace OpenMD {
737
738		RealType mass = sd->getMass();
739		Vector3d vel = sd->getVel();
740	<	RealType value;
740	>	RealType value(0.0);
741
742		switch(rnemdFluxType_) {
743		case rnemdKE :
#	Line 744 \| Line 778 \| namespace OpenMD {
778		if (!min_found) {
779		min_val = value;
780		min_sd = sd;
781	<	min_found = true;
781	>	min_found = 1;
782		} else {
783		if (min_val > value) {
784		min_val = value;
#	Line 857 \| Line 891 \| namespace OpenMD {
891
892		Vector3d min_vel;
893		Vector3d max_vel = max_sd->getVel();
894	<	MPI_Status* status;
894	>	MPI_Status status;
895
896		// point-to-point swap of the velocity vector
897		MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
898		min_vals.rank, 0,
899		min_vel.getArrayPointer(), 3, MPI_REALTYPE,
900	<	min_vals.rank, 0, MPI_COMM_WORLD, status);
900	>	min_vals.rank, 0, MPI_COMM_WORLD, &status);
901
902		switch(rnemdFluxType_) {
903		case rnemdKE :
#	Line 878 \| Line 912 \| namespace OpenMD {
912		MPI_REALTYPE, min_vals.rank, 1,
913		min_angMom.getArrayPointer(), 3,
914		MPI_REALTYPE, min_vals.rank, 1,
915	<	MPI_COMM_WORLD, status);
915	>	MPI_COMM_WORLD, &status);
916
917		max_sd->setJ(min_angMom);
918		}
#	Line 903 \| Line 937 \| namespace OpenMD {
937
938		Vector3d max_vel;
939		Vector3d min_vel = min_sd->getVel();
940	<	MPI_Status* status;
940	>	MPI_Status status;
941
942		// point-to-point swap of the velocity vector
943		MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
944		max_vals.rank, 0,
945		max_vel.getArrayPointer(), 3, MPI_REALTYPE,
946	<	max_vals.rank, 0, MPI_COMM_WORLD, status);
946	>	max_vals.rank, 0, MPI_COMM_WORLD, &status);
947
948		switch(rnemdFluxType_) {
949		case rnemdKE :
#	Line 924 \| Line 958 \| namespace OpenMD {
958		MPI_REALTYPE, max_vals.rank, 1,
959		max_angMom.getArrayPointer(), 3,
960		MPI_REALTYPE, max_vals.rank, 1,
961	<	MPI_COMM_WORLD, status);
961	>	MPI_COMM_WORLD, &status);
962
963		min_sd->setJ(max_angMom);
964		}
#	Line 988 \| Line 1022 \| namespace OpenMD {
1022		int selej;
1023
1024		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
991	–	RealType time = currentSnap_->getTime();
1025		Mat3x3d hmat = currentSnap_->getHmat();
1026
1027		StuntDouble* sd;
#	Line 1115 \| Line 1148 \| namespace OpenMD {
1148		RealType px = Pcx / Phx;
1149		RealType py = Pcy / Phy;
1150		RealType pz = Pcz / Phz;
1151	<	RealType c, x, y, z;
1151	>	RealType c(0.0), x(0.0), y(0.0), z(0.0);
1152		bool successfulScale = false;
1153		if ((rnemdFluxType_ == rnemdFullKE) \|\|
1154		(rnemdFluxType_ == rnemdRotKE)) {
#	Line 1185 \| Line 1218 \| namespace OpenMD {
1218		}
1219		}
1220		} else {
1221	<	RealType a000, a110, c0, a001, a111, b01, b11, c1;
1221	>	RealType a000(0.0), a110(0.0), c0(0.0);
1222	>	RealType a001(0.0), a111(0.0), b01(0.0), b11(0.0), c1(0.0);
1223		switch(rnemdFluxType_) {
1224		case rnemdKE :
1225		/* used hotBin coeff's & only scale x & y dimensions
#	Line 1288 \| Line 1322 \| namespace OpenMD {
1322		vector<pair<RealType,RealType> > rps;
1323		for (ri = realRoots.begin(); ri !=realRoots.end(); ++ri) {
1324		r2 = *ri;
1325	<	//check if FindRealRoots() give the right answer
1325	>	// Check to see if FindRealRoots() gave the right answer:
1326		if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
1327		sprintf(painCave.errMsg,
1328		"RNEMD Warning: polynomial solve seems to have an error!");
#	Line 1296 \| Line 1330 \| namespace OpenMD {
1330		simError();
1331		failRootCount_++;
1332		}
1333	<	//might not be useful w/o rescaling coefficients
1333	>	// Might not be useful w/o rescaling coefficients
1334		alpha0 = -c0 - a110 * r2 * r2;
1335		if (alpha0 >= 0.0) {
1336		r1 = sqrt(alpha0 / a000);
#	Line 1311 \| Line 1345 \| namespace OpenMD {
1345		}
1346		}
1347		}
1348	<	// Consider combining together the solving pair part w/ the searching
1349	<	// best solution part so that we don't need the pairs vector
1348	>	// Consider combining together the part for solving for the pair
1349	>	// w/ the searching for the best solution part so that we don't
1350	>	// need the pairs vector:
1351		if (!rps.empty()) {
1352		RealType smallestDiff = HONKING_LARGE_VALUE;
1353	<	RealType diff;
1354	<	pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
1355	<	vector<pair<RealType,RealType> >::iterator rpi;
1353	>	RealType diff(0.0);
1354	>	std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0);
1355	>	std::vector<std::pair<RealType,RealType> >::iterator rpi;
1356		for (rpi = rps.begin(); rpi != rps.end(); ++rpi) {
1357		r1 = (*rpi).first;
1358		r2 = (*rpi).second;
#	Line 1436 \| Line 1471 \| namespace OpenMD {
1471		int selej;
1472
1473		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1439	–	RealType time = currentSnap_->getTime();
1474		Mat3x3d hmat = currentSnap_->getHmat();
1475
1476		StuntDouble* sd;
#	Line 1743 \| Line 1777 \| namespace OpenMD {
1777		areaA = evaluatorA_.getSurfaceArea();
1778		else {
1779
1746	–	cerr << "selection A did not have surface area, recomputing\n";
1780		int isd;
1781		StuntDouble* sd;
1782		vector<StuntDouble*> aSites;
#	Line 1781 \| Line 1814 \| namespace OpenMD {
1814		}
1815
1816		if (hasSelectionB_) {
1817	<	if (evaluatorB_.hasSurfaceArea())
1817	>	if (evaluatorB_.hasSurfaceArea()) {
1818		areaB = evaluatorB_.getSurfaceArea();
1819	<	else {
1787	<	cerr << "selection B did not have surface area, recomputing\n";
1819	>	} else {
1820
1821		int isd;
1822		StuntDouble* sd;
#	Line 1875 \| Line 1907 \| namespace OpenMD {
1907		void RNEMD::collectData() {
1908		if (!doRNEMD_) return;
1909		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1910	<
1910	>
1911		// collectData can be called more frequently than the doRNEMD, so use the
1912		// computed area from the last exchange time:
1913		RealType area = getDividingArea();
#	Line 1990 \| Line 2022 \| namespace OpenMD {
2022		}
2023		}
2024		}
2025	<
2025	>
2026		#ifdef IS_MPI
2027
2028		for (int i = 0; i < nBins_; i++) {
#	Line 1999 \| Line 2031 \| namespace OpenMD {
2031		1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2032		MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2033		1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2034	<	MPI_Allreduce(MPI_IN_PLACE, &binP[i],
2034	>	MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(),
2035		3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2036	<	MPI_Allreduce(MPI_IN_PLACE, &binL[i],
2036	>	MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(),
2037		3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2038	<	MPI_Allreduce(MPI_IN_PLACE, &binI[i],
2038	>	MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(),
2039		9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2040		MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2041		1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
#	Line 2097 \| Line 2129 \| namespace OpenMD {
2129		painCave.severity = OPENMD_ERROR;
2130		simError();
2131		}
2132	<	}
2132	>	}
2133		}
2134
2135		void RNEMD::writeOutputFile() {
#	Line 2239 \| Line 2271 \| namespace OpenMD {
2271		}
2272
2273		rnemdFile_ << "#######################################################\n";
2274	<	rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
2274	>	rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
2275		rnemdFile_ << "#######################################################\n";
2276
2277
#	Line 2248 \| Line 2280 \| namespace OpenMD {
2280		for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2281		if (outputMask_[i]) {
2282		if (data_[i].dataType == "RealType")
2283	<	writeRealStdDev(i,j);
2283	>	writeRealErrorBars(i,j);
2284		else if (data_[i].dataType == "Vector3d")
2285	<	writeVectorStdDev(i,j);
2285	>	writeVectorErrorBars(i,j);
2286		else {
2287		sprintf( painCave.errMsg,
2288		"RNEMD found an unknown data type for: %s ",
#	Line 2321 \| Line 2353 \| namespace OpenMD {
2353		}
2354		}
2355
2356	<	void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2356	>	void RNEMD::writeRealErrorBars(int index, unsigned int bin) {
2357		if (!doRNEMD_) return;
2358		assert(index >=0 && index < ENDINDEX);
2359		assert(int(bin) < nBins_);
#	Line 2331 \| Line 2363 \| namespace OpenMD {
2363		count = data_[index].accumulator[bin]->count();
2364		if (count == 0) return;
2365
2366	<	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2366	>	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2367
2368		if (! isinf(s) && ! isnan(s)) {
2369		rnemdFile_ << "\t" << s;
#	Line 2344 \| Line 2376 \| namespace OpenMD {
2376		}
2377		}
2378
2379	<	void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2379	>	void RNEMD::writeVectorErrorBars(int index, unsigned int bin) {
2380		if (!doRNEMD_) return;
2381		assert(index >=0 && index < ENDINDEX);
2382		assert(int(bin) < nBins_);
#	Line 2354 \| Line 2386 \| namespace OpenMD {
2386		count = data_[index].accumulator[bin]->count();
2387		if (count == 0) return;
2388
2389	<	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2389	>	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2390		if (isinf(s[0]) \|\| isnan(s[0]) \|\|
2391		isinf(s[1]) \|\| isnan(s[1]) \|\|
2392		isinf(s[2]) \|\| isnan(s[2]) ) {

Diff Legend

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

Comparing trunk/src/rnemd/RNEMD.cpp (file contents): Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC vs. Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

Diff Legend

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC