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

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1946 by gezelter, Tue Nov 12 02:18:35 2013 UTC vs.
Revision 2068 by gezelter, Thu Mar 5 15:40:58 2015 UTC

#	Line 291 \| Line 291 \| namespace OpenMD {
291		kineticFlux_ = 0.0;
292		}
293		if (hasMomentumFluxVector) {
294	<	momentumFluxVector_ = rnemdParams->getMomentumFluxVector();
294	>	std::vector<RealType> mf = rnemdParams->getMomentumFluxVector();
295	>	if (mf.size() != 3) {
296	>	sprintf(painCave.errMsg,
297	>	"RNEMD: Incorrect number of parameters specified for momentumFluxVector.\n"
298	>	"\tthere should be 3 parameters, but %lu were specified.\n",
299	>	mf.size());
300	>	painCave.isFatal = 1;
301	>	simError();
302	>	}
303	>	momentumFluxVector_.x() = mf[0];
304	>	momentumFluxVector_.y() = mf[1];
305	>	momentumFluxVector_.z() = mf[2];
306		} else {
307		momentumFluxVector_ = V3Zero;
308		if (hasMomentumFlux) {
#	Line 316 \| Line 327 \| namespace OpenMD {
327		break;
328		}
329		}
330	<	if (hasAngularMomentumFluxVector) {
331	<	angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
332	<	} else {
333	<	angularMomentumFluxVector_ = V3Zero;
334	<	if (hasAngularMomentumFlux) {
335	<	RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
336	<	switch (rnemdFluxType_) {
337	<	case rnemdLx:
338	<	angularMomentumFluxVector_.x() = angularMomentumFlux;
339	<	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	<	}
348	<	}
349	<
350	<	if (hasCoordinateOrigin) {
351	<	coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
352	<	} else {
353	<	coordinateOrigin_ = V3Zero;
330	>	}
331	>	if (hasAngularMomentumFluxVector) {
332	>	std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
333	>	if (amf.size() != 3) {
334	>	sprintf(painCave.errMsg,
335	>	"RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
336	>	"\tthere should be 3 parameters, but %lu were specified.\n",
337	>	amf.size());
338	>	painCave.isFatal = 1;
339	>	simError();
340		}
341	<
342	<	// do some sanity checking
343	<
344	<	int selectionCount = seleMan_.getSelectionCount();
345	<
346	<	int nIntegrable = info->getNGlobalIntegrableObjects();
347	<
362	<	if (selectionCount > nIntegrable) {
363	<	sprintf(painCave.errMsg,
364	<	"RNEMD: The current objectSelection,\n"
365	<	"\t\t%s\n"
366	<	"\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);
341	>	angularMomentumFluxVector_.x() = amf[0];
342	>	angularMomentumFluxVector_.y() = amf[1];
343	>	angularMomentumFluxVector_.z() = amf[2];
344	>	} else {
345	>	angularMomentumFluxVector_ = V3Zero;
346	>	if (hasAngularMomentumFlux) {
347	>	RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
348		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;
349		case rnemdLx:
350	+	angularMomentumFluxVector_.x() = angularMomentumFlux;
351	+	break;
352		case rnemdLy:
353	+	angularMomentumFluxVector_.y() = angularMomentumFlux;
354	+	break;
355		case rnemdLz:
356	<	case rnemdLvector:
470	<	outputMask_.set(ANGULARVELOCITY);
356	>	angularMomentumFluxVector_.z() = angularMomentumFlux;
357		break;
358		case rnemdKeLx:
359	+	angularMomentumFluxVector_.x() = angularMomentumFlux;
360	+	break;
361		case rnemdKeLy:
362	+	angularMomentumFluxVector_.y() = angularMomentumFlux;
363	+	break;
364		case rnemdKeLz:
365	<	case rnemdKeLvector:
476	<	outputMask_.set(TEMPERATURE);
477	<	outputMask_.set(ANGULARVELOCITY);
365	>	angularMomentumFluxVector_.z() = angularMomentumFlux;
366		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;
367		default:
368		break;
369		}
370	<	}
371	<
494	<	if (hasOutputFileName) {
495	<	rnemdFileName_ = rnemdParams->getOutputFileName();
496	<	} else {
497	<	rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
498	<	}
370	>	}
371	>	}
372
373	<	exchangeTime_ = rnemdParams->getExchangeTime();
374	<
375	<	Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
376	<	// total exchange sums are zeroed out at the beginning:
377	<
378	<	kineticExchange_ = 0.0;
379	<	momentumExchange_ = V3Zero;
380	<	angularMomentumExchange_ = V3Zero;
381	<
382	<	std::ostringstream selectionAstream;
383	<	std::ostringstream selectionBstream;
373	>	if (hasCoordinateOrigin) {
374	>	std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
375	>	if (co.size() != 3) {
376	>	sprintf(painCave.errMsg,
377	>	"RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
378	>	"\tthere should be 3 parameters, but %lu were specified.\n",
379	>	co.size());
380	>	painCave.isFatal = 1;
381	>	simError();
382	>	}
383	>	coordinateOrigin_.x() = co[0];
384	>	coordinateOrigin_.y() = co[1];
385	>	coordinateOrigin_.z() = co[2];
386	>	} else {
387	>	coordinateOrigin_ = V3Zero;
388	>	}
389
390	<	if (hasSelectionA_) {
391	<	selectionA_ = rnemdParams->getSelectionA();
392	<	} else {
393	<	if (usePeriodicBoundaryConditions_) {
394	<	Mat3x3d hmat = currentSnap_->getHmat();
395	<
396	<	if (hasSlabWidth)
397	<	slabWidth_ = rnemdParams->getSlabWidth();
398	<	else
399	<	slabWidth_ = hmat(2,2) / 10.0;
400	<
401	<	if (hasSlabACenter)
402	<	slabACenter_ = rnemdParams->getSlabACenter();
403	<	else
404	<	slabACenter_ = 0.0;
405	<
406	<	selectionAstream << "select wrappedz > "
407	<	<< slabACenter_ - 0.5*slabWidth_
408	<	<< " && wrappedz < "
409	<	<< slabACenter_ + 0.5*slabWidth_;
410	<	selectionA_ = selectionAstream.str();
411	<	} else {
412	<	if (hasSphereARadius)
413	<	sphereARadius_ = rnemdParams->getSphereARadius();
414	<	else {
415	<	// use an initial guess to the size of the inner slab to be 1/10 the
416	<	// radius of an approximately spherical hull:
417	<	Thermo thermo(info);
418	<	RealType hVol = thermo.getHullVolume();
419	<	sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
420	<	}
421	<	selectionAstream << "select r < " << sphereARadius_;
422	<	selectionA_ = selectionAstream.str();
423	<	}
390	>	// do some sanity checking
391	>
392	>	int selectionCount = seleMan_.getSelectionCount();
393	>	int nIntegrable = info->getNGlobalIntegrableObjects();
394	>	if (selectionCount > nIntegrable) {
395	>	sprintf(painCave.errMsg,
396	>	"RNEMD: The current objectSelection,\n"
397	>	"\t\t%s\n"
398	>	"\thas resulted in %d selected objects. However,\n"
399	>	"\tthe total number of integrable objects in the system\n"
400	>	"\tis only %d. This is almost certainly not what you want\n"
401	>	"\tto do. A likely cause of this is forgetting the _RB_0\n"
402	>	"\tselector in the selection script!\n",
403	>	rnemdObjectSelection_.c_str(),
404	>	selectionCount, nIntegrable);
405	>	painCave.isFatal = 0;
406	>	painCave.severity = OPENMD_WARNING;
407	>	simError();
408	>	}
409	>
410	>	areaAccumulator_ = new Accumulator();
411	>
412	>	nBins_ = rnemdParams->getOutputBins();
413	>	binWidth_ = rnemdParams->getOutputBinWidth();
414	>
415	>	data_.resize(RNEMD::ENDINDEX);
416	>	OutputData z;
417	>	z.units = "Angstroms";
418	>	z.title = "Z";
419	>	z.dataType = "RealType";
420	>	z.accumulator.reserve(nBins_);
421	>	for (int i = 0; i < nBins_; i++)
422	>	z.accumulator.push_back( new Accumulator() );
423	>	data_[Z] = z;
424	>	outputMap_["Z"] = Z;
425	>
426	>	OutputData r;
427	>	r.units = "Angstroms";
428	>	r.title = "R";
429	>	r.dataType = "RealType";
430	>	r.accumulator.reserve(nBins_);
431	>	for (int i = 0; i < nBins_; i++)
432	>	r.accumulator.push_back( new Accumulator() );
433	>	data_[R] = r;
434	>	outputMap_["R"] = R;
435	>
436	>	OutputData temperature;
437	>	temperature.units = "K";
438	>	temperature.title = "Temperature";
439	>	temperature.dataType = "RealType";
440	>	temperature.accumulator.reserve(nBins_);
441	>	for (int i = 0; i < nBins_; i++)
442	>	temperature.accumulator.push_back( new Accumulator() );
443	>	data_[TEMPERATURE] = temperature;
444	>	outputMap_["TEMPERATURE"] = TEMPERATURE;
445	>
446	>	OutputData velocity;
447	>	velocity.units = "angstroms/fs";
448	>	velocity.title = "Velocity";
449	>	velocity.dataType = "Vector3d";
450	>	velocity.accumulator.reserve(nBins_);
451	>	for (int i = 0; i < nBins_; i++)
452	>	velocity.accumulator.push_back( new VectorAccumulator() );
453	>	data_[VELOCITY] = velocity;
454	>	outputMap_["VELOCITY"] = VELOCITY;
455	>
456	>	OutputData angularVelocity;
457	>	angularVelocity.units = "angstroms^2/fs";
458	>	angularVelocity.title = "AngularVelocity";
459	>	angularVelocity.dataType = "Vector3d";
460	>	angularVelocity.accumulator.reserve(nBins_);
461	>	for (int i = 0; i < nBins_; i++)
462	>	angularVelocity.accumulator.push_back( new VectorAccumulator() );
463	>	data_[ANGULARVELOCITY] = angularVelocity;
464	>	outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
465	>
466	>	OutputData density;
467	>	density.units = "g cm^-3";
468	>	density.title = "Density";
469	>	density.dataType = "RealType";
470	>	density.accumulator.reserve(nBins_);
471	>	for (int i = 0; i < nBins_; i++)
472	>	density.accumulator.push_back( new Accumulator() );
473	>	data_[DENSITY] = density;
474	>	outputMap_["DENSITY"] = DENSITY;
475	>
476	>	if (hasOutputFields) {
477	>	parseOutputFileFormat(rnemdParams->getOutputFields());
478	>	} else {
479	>	if (usePeriodicBoundaryConditions_)
480	>	outputMask_.set(Z);
481	>	else
482	>	outputMask_.set(R);
483	>	switch (rnemdFluxType_) {
484	>	case rnemdKE:
485	>	case rnemdRotKE:
486	>	case rnemdFullKE:
487	>	outputMask_.set(TEMPERATURE);
488	>	break;
489	>	case rnemdPx:
490	>	case rnemdPy:
491	>	outputMask_.set(VELOCITY);
492	>	break;
493	>	case rnemdPz:
494	>	case rnemdPvector:
495	>	outputMask_.set(VELOCITY);
496	>	outputMask_.set(DENSITY);
497	>	break;
498	>	case rnemdLx:
499	>	case rnemdLy:
500	>	case rnemdLz:
501	>	case rnemdLvector:
502	>	outputMask_.set(ANGULARVELOCITY);
503	>	break;
504	>	case rnemdKeLx:
505	>	case rnemdKeLy:
506	>	case rnemdKeLz:
507	>	case rnemdKeLvector:
508	>	outputMask_.set(TEMPERATURE);
509	>	outputMask_.set(ANGULARVELOCITY);
510	>	break;
511	>	case rnemdKePx:
512	>	case rnemdKePy:
513	>	outputMask_.set(TEMPERATURE);
514	>	outputMask_.set(VELOCITY);
515	>	break;
516	>	case rnemdKePvector:
517	>	outputMask_.set(TEMPERATURE);
518	>	outputMask_.set(VELOCITY);
519	>	outputMask_.set(DENSITY);
520	>	break;
521	>	default:
522	>	break;
523		}
524	+	}
525
526	<	if (hasSelectionB_) {
527	<	selectionB_ = rnemdParams->getSelectionB();
528	<
529	<	} else {
530	<	if (usePeriodicBoundaryConditions_) {
531	<	Mat3x3d hmat = currentSnap_->getHmat();
526	>	if (hasOutputFileName) {
527	>	rnemdFileName_ = rnemdParams->getOutputFileName();
528	>	} else {
529	>	rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
530	>	}
531	>
532	>	exchangeTime_ = rnemdParams->getExchangeTime();
533	>
534	>	Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
535	>	// total exchange sums are zeroed out at the beginning:
536	>
537	>	kineticExchange_ = 0.0;
538	>	momentumExchange_ = V3Zero;
539	>	angularMomentumExchange_ = V3Zero;
540	>
541	>	std::ostringstream selectionAstream;
542	>	std::ostringstream selectionBstream;
543	>
544	>	if (hasSelectionA_) {
545	>	selectionA_ = rnemdParams->getSelectionA();
546	>	} else {
547	>	if (usePeriodicBoundaryConditions_) {
548	>	Mat3x3d hmat = currentSnap_->getHmat();
549
550	<	if (hasSlabWidth)
551	<	slabWidth_ = rnemdParams->getSlabWidth();
552	<	else
553	<	slabWidth_ = hmat(2,2) / 10.0;
550	>	if (hasSlabWidth)
551	>	slabWidth_ = rnemdParams->getSlabWidth();
552	>	else
553	>	slabWidth_ = hmat(2,2) / 10.0;
554
555	<	if (hasSlabBCenter)
556	<	slabBCenter_ = rnemdParams->getSlabBCenter();
557	<	else
558	<	slabBCenter_ = hmat(2,2) / 2.0;
555	>	if (hasSlabACenter)
556	>	slabACenter_ = rnemdParams->getSlabACenter();
557	>	else
558	>	slabACenter_ = 0.0;
559
560	<	selectionBstream << "select wrappedz > "
561	<	<< slabBCenter_ - 0.5*slabWidth_
562	<	<< " && wrappedz < "
563	<	<< slabBCenter_ + 0.5*slabWidth_;
560	>	selectionAstream << "select wrappedz > "
561	>	<< slabACenter_ - 0.5*slabWidth_
562	>	<< " && wrappedz < "
563	>	<< slabACenter_ + 0.5*slabWidth_;
564	>	selectionA_ = selectionAstream.str();
565	>	} else {
566	>	if (hasSphereARadius)
567	>	sphereARadius_ = rnemdParams->getSphereARadius();
568	>	else {
569	>	// use an initial guess to the size of the inner slab to be 1/10 the
570	>	// radius of an approximately spherical hull:
571	>	Thermo thermo(info);
572	>	RealType hVol = thermo.getHullVolume();
573	>	sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
574	>	}
575	>	selectionAstream << "select r < " << sphereARadius_;
576	>	selectionA_ = selectionAstream.str();
577	>	}
578	>	}
579	>
580	>	if (hasSelectionB_) {
581	>	selectionB_ = rnemdParams->getSelectionB();
582	>
583	>	} else {
584	>	if (usePeriodicBoundaryConditions_) {
585	>	Mat3x3d hmat = currentSnap_->getHmat();
586	>
587	>	if (hasSlabWidth)
588	>	slabWidth_ = rnemdParams->getSlabWidth();
589	>	else
590	>	slabWidth_ = hmat(2,2) / 10.0;
591	>
592	>	if (hasSlabBCenter)
593	>	slabBCenter_ = rnemdParams->getSlabBCenter();
594	>	else
595	>	slabBCenter_ = hmat(2,2) / 2.0;
596	>
597	>	selectionBstream << "select wrappedz > "
598	>	<< slabBCenter_ - 0.5*slabWidth_
599	>	<< " && wrappedz < "
600	>	<< slabBCenter_ + 0.5*slabWidth_;
601	>	selectionB_ = selectionBstream.str();
602	>	} else {
603	>	if (hasSphereBRadius_) {
604	>	sphereBRadius_ = rnemdParams->getSphereBRadius();
605	>	selectionBstream << "select r > " << sphereBRadius_;
606		selectionB_ = selectionBstream.str();
607		} else {
608	<	if (hasSphereBRadius_) {
609	<	sphereBRadius_ = rnemdParams->getSphereBRadius();
610	<	selectionBstream << "select r > " << sphereBRadius_;
574	<	selectionB_ = selectionBstream.str();
575	<	} else {
576	<	selectionB_ = "select hull";
577	<	BisHull_ = true;
578	<	hasSelectionB_ = true;
579	<	}
608	>	selectionB_ = "select hull";
609	>	BisHull_ = true;
610	>	hasSelectionB_ = true;
611		}
612		}
613		}
614	<
614	>
615	>
616		// object evaluator:
617		evaluator_.loadScriptString(rnemdObjectSelection_);
618		seleMan_.setSelectionSet(evaluator_.evaluate());
#	Line 622 \| Line 654 \| namespace OpenMD {
654
655		StuntDouble* sd;
656
657	<	RealType min_val;
658	<	bool min_found = false;
659	<	StuntDouble* min_sd;
657	>	RealType min_val(0.0);
658	>	int min_found = 0;
659	>	StuntDouble* min_sd = NULL;
660
661	<	RealType max_val;
662	<	bool max_found = false;
663	<	StuntDouble* max_sd;
661	>	RealType max_val(0.0);
662	>	int max_found = 0;
663	>	StuntDouble* max_sd = NULL;
664
665		for (sd = seleManA_.beginSelected(selei); sd != NULL;
666		sd = seleManA_.nextSelected(selei)) {
#	Line 682 \| Line 714 \| namespace OpenMD {
714		if (!max_found) {
715		max_val = value;
716		max_sd = sd;
717	<	max_found = true;
717	>	max_found = 1;
718		} else {
719		if (max_val < value) {
720		max_val = value;
#	Line 744 \| Line 776 \| namespace OpenMD {
776		if (!min_found) {
777		min_val = value;
778		min_sd = sd;
779	<	min_found = true;
779	>	min_found = 1;
780		} else {
781		if (min_val > value) {
782		min_val = value;
#	Line 754 \| Line 786 \| namespace OpenMD {
786		}
787
788		#ifdef IS_MPI
789	<	int worldRank = MPI::COMM_WORLD.Get_rank();
790	<
791	<	bool my_min_found = min_found;
792	<	bool my_max_found = max_found;
789	>	int worldRank;
790	>	MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
791	>
792	>	int my_min_found = min_found;
793	>	int my_max_found = max_found;
794
795		// Even if we didn't find a minimum, did someone else?
796	<	MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
796	>	MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
797	>	MPI_COMM_WORLD);
798		// Even if we didn't find a maximum, did someone else?
799	<	MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
799	>	MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
800	>	MPI_COMM_WORLD);
801		#endif
802
803		if (max_found && min_found) {
#	Line 781 \| Line 816 \| namespace OpenMD {
816		min_vals.rank = worldRank;
817
818		// Who had the minimum?
819	<	MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
820	<	1, MPI::REALTYPE_INT, MPI::MINLOC);
819	>	MPI_Allreduce(&min_vals, &min_vals,
820	>	1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD);
821		min_val = min_vals.val;
822
823		if (my_max_found) {
#	Line 793 \| Line 828 \| namespace OpenMD {
828		max_vals.rank = worldRank;
829
830		// Who had the maximum?
831	<	MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
832	<	1, MPI::REALTYPE_INT, MPI::MAXLOC);
831	>	MPI_Allreduce(&max_vals, &max_vals,
832	>	1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
833		max_val = max_vals.val;
834		#endif
835
#	Line 854 \| Line 889 \| namespace OpenMD {
889
890		Vector3d min_vel;
891		Vector3d max_vel = max_sd->getVel();
892	<	MPI::Status status;
892	>	MPI_Status status;
893
894		// point-to-point swap of the velocity vector
895	<	MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
896	<	min_vals.rank, 0,
897	<	min_vel.getArrayPointer(), 3, MPI::REALTYPE,
898	<	min_vals.rank, 0, status);
895	>	MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
896	>	min_vals.rank, 0,
897	>	min_vel.getArrayPointer(), 3, MPI_REALTYPE,
898	>	min_vals.rank, 0, MPI_COMM_WORLD, &status);
899
900		switch(rnemdFluxType_) {
901		case rnemdKE :
#	Line 871 \| Line 906 \| namespace OpenMD {
906		Vector3d max_angMom = max_sd->getJ();
907
908		// point-to-point swap of the angular momentum vector
909	<	MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
910	<	MPI::REALTYPE, min_vals.rank, 1,
911	<	min_angMom.getArrayPointer(), 3,
912	<	MPI::REALTYPE, min_vals.rank, 1,
913	<	status);
909	>	MPI_Sendrecv(max_angMom.getArrayPointer(), 3,
910	>	MPI_REALTYPE, min_vals.rank, 1,
911	>	min_angMom.getArrayPointer(), 3,
912	>	MPI_REALTYPE, min_vals.rank, 1,
913	>	MPI_COMM_WORLD, &status);
914
915		max_sd->setJ(min_angMom);
916		}
#	Line 900 \| Line 935 \| namespace OpenMD {
935
936		Vector3d max_vel;
937		Vector3d min_vel = min_sd->getVel();
938	<	MPI::Status status;
938	>	MPI_Status status;
939
940		// point-to-point swap of the velocity vector
941	<	MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
942	<	max_vals.rank, 0,
943	<	max_vel.getArrayPointer(), 3, MPI::REALTYPE,
944	<	max_vals.rank, 0, status);
941	>	MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
942	>	max_vals.rank, 0,
943	>	max_vel.getArrayPointer(), 3, MPI_REALTYPE,
944	>	max_vals.rank, 0, MPI_COMM_WORLD, &status);
945
946		switch(rnemdFluxType_) {
947		case rnemdKE :
#	Line 917 \| Line 952 \| namespace OpenMD {
952		Vector3d max_angMom;
953
954		// point-to-point swap of the angular momentum vector
955	<	MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
956	<	MPI::REALTYPE, max_vals.rank, 1,
957	<	max_angMom.getArrayPointer(), 3,
958	<	MPI::REALTYPE, max_vals.rank, 1,
959	<	status);
955	>	MPI_Sendrecv(min_angMom.getArrayPointer(), 3,
956	>	MPI_REALTYPE, max_vals.rank, 1,
957	>	max_angMom.getArrayPointer(), 3,
958	>	MPI_REALTYPE, max_vals.rank, 1,
959	>	MPI_COMM_WORLD, &status);
960
961		min_sd->setJ(max_angMom);
962		}
#	Line 985 \| Line 1020 \| namespace OpenMD {
1020		int selej;
1021
1022		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
988	–	RealType time = currentSnap_->getTime();
1023		Mat3x3d hmat = currentSnap_->getHmat();
1024
1025		StuntDouble* sd;
#	Line 1090 \| Line 1124 \| namespace OpenMD {
1124		Kcw *= 0.5;
1125
1126		#ifdef IS_MPI
1127	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
1128	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
1129	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
1130	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
1131	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
1132	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
1127	>	MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1128	>	MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1129	>	MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1130	>	MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1131	>	MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1132	>	MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1133
1134	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
1135	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
1136	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
1137	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM);
1134	>	MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1135	>	MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1136	>	MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1137	>	MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1138
1139	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
1140	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
1141	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
1142	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM);
1139	>	MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1140	>	MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1141	>	MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1142	>	MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1143		#endif
1144
1145		//solve coldBin coeff's first
#	Line 1433 \| Line 1467 \| namespace OpenMD {
1467		int selej;
1468
1469		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1436	–	RealType time = currentSnap_->getTime();
1470		Mat3x3d hmat = currentSnap_->getHmat();
1471
1472		StuntDouble* sd;
#	Line 1582 \| Line 1615 \| namespace OpenMD {
1615		Kc *= 0.5;
1616
1617		#ifdef IS_MPI
1618	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1619	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
1620	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM);
1621	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM);
1622	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM);
1623	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM);
1624	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM);
1625	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM);
1626	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9,
1627	<	MPI::REALTYPE, MPI::SUM);
1628	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9,
1629	<	MPI::REALTYPE, MPI::SUM);
1618	>	MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,
1619	>	MPI_COMM_WORLD);
1620	>	MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,
1621	>	MPI_COMM_WORLD);
1622	>	MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,
1623	>	MPI_COMM_WORLD);
1624	>	MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,
1625	>	MPI_COMM_WORLD);
1626	>	MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1627	>	MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1628	>	MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1629	>	MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1630	>	MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,
1631	>	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1632	>	MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,
1633	>	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1634		#endif
1635
1636
#	Line 1736 \| Line 1773 \| namespace OpenMD {
1773		areaA = evaluatorA_.getSurfaceArea();
1774		else {
1775
1739	–	cerr << "selection A did not have surface area, recomputing\n";
1776		int isd;
1777		StuntDouble* sd;
1778		vector<StuntDouble*> aSites;
#	Line 1774 \| Line 1810 \| namespace OpenMD {
1810		}
1811
1812		if (hasSelectionB_) {
1813	<	if (evaluatorB_.hasSurfaceArea())
1813	>	if (evaluatorB_.hasSurfaceArea()) {
1814		areaB = evaluatorB_.getSurfaceArea();
1815	<	else {
1780	<	cerr << "selection B did not have surface area, recomputing\n";
1815	>	} else {
1816
1817		int isd;
1818		StuntDouble* sd;
#	Line 1868 \| Line 1903 \| namespace OpenMD {
1903		void RNEMD::collectData() {
1904		if (!doRNEMD_) return;
1905		Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1906	<
1906	>
1907		// collectData can be called more frequently than the doRNEMD, so use the
1908		// computed area from the last exchange time:
1909		RealType area = getDividingArea();
#	Line 1983 \| Line 2018 \| namespace OpenMD {
2018		}
2019		}
2020		}
2021	<
2021	>
2022		#ifdef IS_MPI
2023
2024		for (int i = 0; i < nBins_; i++) {
2025	<
2026	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[i],
2027	<	1, MPI::INT, MPI::SUM);
2028	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[i],
2029	<	1, MPI::REALTYPE, MPI::SUM);
2030	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binP[i],
2031	<	3, MPI::REALTYPE, MPI::SUM);
2032	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binL[i],
2033	<	3, MPI::REALTYPE, MPI::SUM);
2034	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binI[i],
2035	<	9, MPI::REALTYPE, MPI::SUM);
2036	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[i],
2037	<	1, MPI::REALTYPE, MPI::SUM);
2038	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[i],
2039	<	1, MPI::INT, MPI::SUM);
2040	<	//MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmega[i],
2041	<	// 1, MPI::REALTYPE, MPI::SUM);
2025	>
2026	>	MPI_Allreduce(MPI_IN_PLACE, &binCount[i],
2027	>	1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2028	>	MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2029	>	1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2030	>	MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(),
2031	>	3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2032	>	MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(),
2033	>	3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2034	>	MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(),
2035	>	9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2036	>	MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2037	>	1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2038	>	MPI_Allreduce(MPI_IN_PLACE, &binDOF[i],
2039	>	1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2040	>	//MPI_Allreduce(MPI_IN_PLACE, &binOmega[i],
2041	>	// 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2042		}
2043
2044		#endif
#	Line 2090 \| Line 2125 \| namespace OpenMD {
2125		painCave.severity = OPENMD_ERROR;
2126		simError();
2127		}
2128	<	}
2128	>	}
2129		}
2130
2131		void RNEMD::writeOutputFile() {
#	Line 2099 \| Line 2134 \| namespace OpenMD {
2134
2135		#ifdef IS_MPI
2136		// If we're the root node, should we print out the results
2137	<	int worldRank = MPI::COMM_WORLD.Get_rank();
2137	>	int worldRank;
2138	>	MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
2139	>
2140		if (worldRank == 0) {
2141		#endif
2142		rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out \| std::ios::trunc );
#	Line 2230 \| Line 2267 \| namespace OpenMD {
2267		}
2268
2269		rnemdFile_ << "#######################################################\n";
2270	<	rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
2270	>	rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
2271		rnemdFile_ << "#######################################################\n";
2272
2273
#	Line 2239 \| Line 2276 \| namespace OpenMD {
2276		for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2277		if (outputMask_[i]) {
2278		if (data_[i].dataType == "RealType")
2279	<	writeRealStdDev(i,j);
2279	>	writeRealErrorBars(i,j);
2280		else if (data_[i].dataType == "Vector3d")
2281	<	writeVectorStdDev(i,j);
2281	>	writeVectorErrorBars(i,j);
2282		else {
2283		sprintf( painCave.errMsg,
2284		"RNEMD found an unknown data type for: %s ",
#	Line 2312 \| Line 2349 \| namespace OpenMD {
2349		}
2350		}
2351
2352	<	void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2352	>	void RNEMD::writeRealErrorBars(int index, unsigned int bin) {
2353		if (!doRNEMD_) return;
2354		assert(index >=0 && index < ENDINDEX);
2355		assert(int(bin) < nBins_);
#	Line 2322 \| Line 2359 \| namespace OpenMD {
2359		count = data_[index].accumulator[bin]->count();
2360		if (count == 0) return;
2361
2362	<	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2362	>	dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2363
2364		if (! isinf(s) && ! isnan(s)) {
2365		rnemdFile_ << "\t" << s;
#	Line 2335 \| Line 2372 \| namespace OpenMD {
2372		}
2373		}
2374
2375	<	void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2375	>	void RNEMD::writeVectorErrorBars(int index, unsigned int bin) {
2376		if (!doRNEMD_) return;
2377		assert(index >=0 && index < ENDINDEX);
2378		assert(int(bin) < nBins_);
#	Line 2345 \| Line 2382 \| namespace OpenMD {
2382		count = data_[index].accumulator[bin]->count();
2383		if (count == 0) return;
2384
2385	<	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2385	>	dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2386		if (isinf(s[0]) \|\| isnan(s[0]) \|\|
2387		isinf(s[1]) \|\| isnan(s[1]) \|\|
2388		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 1946 by gezelter, Tue Nov 12 02:18:35 2013 UTC vs. Revision 2068 by gezelter, Thu Mar 5 15:40:58 2015 UTC

Diff Legend

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1946 by gezelter, Tue Nov 12 02:18:35 2013 UTC vs.
Revision 2068 by gezelter, Thu Mar 5 15:40:58 2015 UTC