45#include "applications/staticProps/RNEMDStats.hpp"
56#include "applications/staticProps/SpatialStatistics.hpp"
57#include "brains/DataStorage.hpp"
59#include "io/Globals.hpp"
66#include "rnemd/RNEMDParameters.hpp"
67#include "types/AtomType.hpp"
68#include "types/FixedChargeAdapter.hpp"
69#include "types/FluctuatingChargeAdapter.hpp"
70#include "utils/Accumulator.hpp"
71#include "utils/AccumulatorView.hpp"
72#include "utils/BaseAccumulator.hpp"
73#include "utils/Constants.hpp"
76using namespace OpenMD::Utils;
80 RNEMDZ::RNEMDZ(
SimInfo* info,
const std::string& filename,
81 const std::string& sele,
int nzbins,
int axis) :
83 setOutputName(
getPrefix(filename) +
".rnemdZ");
85 evaluator_.loadScriptString(sele);
86 seleMan_.setSelectionSet(evaluator_.evaluate());
87 AtomTypeSet osTypes = seleMan_.getSelectedAtomTypes();
88 std::copy(osTypes.begin(), osTypes.end(), std::back_inserter(outputTypes_));
89 bool usePeriodicBoundaryConditions_ =
90 info_->getSimParams()->getUsePeriodicBoundaryConditions();
92 data_.resize(RNEMDZ::ENDINDEX);
95 z.units =
"Angstroms";
97 z.dataHandling = DataHandling::Average;
98 for (
unsigned int i = 0; i < nBins_; i++)
99 z.accumulator.push_back(
100 std::make_unique<AccumulatorView<RealAccumulator>>());
101 data_[Z] = std::move(z);
103 OutputData temperature;
104 temperature.units =
"K";
105 temperature.title =
"Temperature";
106 temperature.dataHandling = DataHandling::Average;
107 for (
unsigned int i = 0; i < nBins_; i++)
108 temperature.accumulator.push_back(
109 std::make_unique<AccumulatorView<RealAccumulator>>());
110 data_[TEMPERATURE] = std::move(temperature);
113 velocity.units =
"angstroms/fs";
114 velocity.title =
"Velocity";
115 velocity.dataHandling = DataHandling::Average;
116 for (
unsigned int i = 0; i < nBins_; i++)
117 velocity.accumulator.push_back(
118 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
119 data_[VELOCITY] = std::move(velocity);
122 density.units =
"g cm^-3";
123 density.title =
"Density";
124 density.dataHandling = DataHandling::Average;
125 for (
unsigned int i = 0; i < nBins_; i++)
126 density.accumulator.push_back(
127 std::make_unique<AccumulatorView<RealAccumulator>>());
128 data_[DENSITY] = std::move(density);
131 activity.units =
"unitless";
132 activity.title =
"Activity";
133 activity.dataHandling = DataHandling::Average;
134 unsigned int nTypes = outputTypes_.size();
137 for (
unsigned int i = 0; i < nBins_; i++)
138 activity.accumulator.push_back(
139 std::make_unique<AccumulatorView<StdVectorAccumulator>>());
140 data_[ACTIVITY] = std::move(activity);
144 eField.units =
"kcal/mol/angstroms/e";
145 eField.title =
"Electric Field";
146 eField.dataHandling = DataHandling::Average;
147 for (
unsigned int i = 0; i < nBins_; i++)
148 eField.accumulator.push_back(
149 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
152 ePot.units =
"kcal/mol/e";
153 ePot.title =
"Electrostatic Potential";
154 ePot.dataHandling = DataHandling::Average;
155 for (
unsigned int i = 0; i < nBins_; i++)
156 ePot.accumulator.push_back(
157 std::make_unique<AccumulatorView<RealAccumulator>>());
161 charge.title =
"Charge";
162 charge.dataHandling = DataHandling::Average;
163 for (
unsigned int i = 0; i < nBins_; i++)
164 charge.accumulator.push_back(
165 std::make_unique<AccumulatorView<RealAccumulator>>());
167 OutputData chargeVelocity;
168 chargeVelocity.units =
"e/fs";
169 chargeVelocity.title =
"Charge_Velocity";
170 chargeVelocity.dataHandling = DataHandling::Average;
171 for (
unsigned int i = 0; i < nBins_; i++)
172 chargeVelocity.accumulator.push_back(
173 std::make_unique<AccumulatorView<RealAccumulator>>());
176 outputMask_.set(TEMPERATURE);
177 outputMask_.set(VELOCITY);
178 outputMask_.set(DENSITY);
179 outputMask_.set(ACTIVITY);
181 int atomStorageLayout = info_->getAtomStorageLayout();
182 int rigidBodyStorageLayout = info->getRigidBodyStorageLayout();
183 int cutoffGroupStorageLayout = info->getCutoffGroupStorageLayout();
185 if (atomStorageLayout & DataStorage::dslElectricField) {
186 outputMask_.set(ELECTRICFIELD);
187 outputMask_.set(ELECTROSTATICPOTENTIAL);
189 data_[ELECTRICFIELD] = std::move(eField);
190 data_[ELECTROSTATICPOTENTIAL] = std::move(ePot);
193 if (info_->usesElectrostaticAtoms() ||
194 atomStorageLayout & DataStorage::dslFlucQPosition) {
195 outputMask_.set(CHARGE);
197 data_[CHARGE] = std::move(charge);
200 if (atomStorageLayout & DataStorage::dslFlucQVelocity) {
201 outputMask_.set(CHARGEVELOCITY);
203 data_[CHARGEVELOCITY] = std::move(chargeVelocity);
207 void RNEMDZ::processFrame(
int istep) {
208 SlabStatistics::processFrame(istep);
211 seleMan_.setSelectionSet(evaluator_.evaluate());
224 std::vector<RealType> binMass(nBins_, 0.0);
225 std::vector<Vector3d> binP(nBins_, V3Zero);
226 std::vector<RealType> binCharge(nBins_, 0.0);
227 std::vector<RealType> binChargeVelocity(nBins_, 0.0);
228 std::vector<RealType> binKE(nBins_, 0.0);
229 std::vector<Vector3d> binEField(nBins_, V3Zero);
230 std::vector<int> binDOF(nBins_, 0);
231 std::vector<int> binCount(nBins_, 0);
232 std::vector<std::vector<int>> binTypeCounts;
233 std::vector<int> binEFieldCount(nBins_, 0);
235 if (outputMask_[ACTIVITY]) {
236 binTypeCounts.resize(nBins_);
237 for (
unsigned int i = 0; i < nBins_; i++) {
238 binTypeCounts[i].resize(outputTypes_.size(), 0);
242 SimInfo::MoleculeIterator miter;
243 std::vector<StuntDouble*>::iterator iiter;
244 std::vector<AtomType*>::iterator at;
249 Molecule::ConstraintPairIterator cpi;
253 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
254 sd = mol->nextIntegrableObject(iiter)) {
255 if (seleMan_.isSelected(sd)) {
271 KE += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
272 angMom[k] * angMom[k] / Ia(k, k));
275 KE += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
276 angMom[1] * angMom[1] / Ia(1, 1) +
277 angMom[2] * angMom[2] / Ia(2, 2));
282 if (outputMask_[ACTIVITY]) {
285 atype =
static_cast<Atom*
>(sd)->getAtomType();
286 at = std::find(outputTypes_.begin(), outputTypes_.end(), atype);
287 if (at != outputTypes_.end()) {
288 typeIndex = std::distance(outputTypes_.begin(), at);
293 if (binNo >= 0 && binNo <
int(nBins_)) {
295 binMass[binNo] += mass;
296 binP[binNo] += mass * vel;
298 binDOF[binNo] += DOF;
300 if (outputMask_[ACTIVITY] && typeIndex != -1)
301 binTypeCounts[binNo][typeIndex]++;
303 if (outputMask_[CHARGE] || outputMask_[CHARGEVELOCITY]) {
305 AtomType* atomType =
static_cast<Atom*
>(sd)->getAtomType();
307 if (fca.isFixedCharge()) { q = fca.getCharge(); }
310 if (fqa.isFluctuatingCharge()) {
315 if (outputMask_[CHARGE]) binCharge[binNo] += q;
316 if (outputMask_[CHARGEVELOCITY]) binChargeVelocity[binNo] += w;
319 std::vector<Atom*>::iterator ai;
321 for (atom = rb->beginAtom(ai); atom != NULL;
322 atom = rb->nextAtom(ai)) {
323 binNo = getBin(atom->getPos());
324 AtomType* atomType = atom->getAtomType();
326 if (fca.isFixedCharge()) { q = fca.getCharge(); }
330 if (fqa.isFluctuatingCharge()) {
335 if (outputMask_[CHARGE]) binCharge[binNo] += q;
336 if (outputMask_[CHARGEVELOCITY])
337 binChargeVelocity[binNo] += w;
346 if (outputMask_[ELECTRICFIELD]) {
349 std::vector<Atom*>::iterator ai;
351 for (atom = rb->beginAtom(ai); atom != NULL;
352 atom = rb->nextAtom(ai)) {
353 binNo = getBin(atom->getPos());
354 eField = atom->getElectricField();
355 binEFieldCount[binNo]++;
356 binEField[binNo] += eField;
360 binNo = getBin(sd->
getPos());
362 binEFieldCount[binNo]++;
363 binEField[binNo] += eField;
367 if (seleMan_.isSelected(mol)) {
368 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
369 consPair = mol->nextConstraintPair(cpi)) {
370 Vector3d posA = consPair->getConsElem1()->getPos();
371 Vector3d posB = consPair->getConsElem2()->getPos();
373 if (usePeriodicBoundaryConditions_) {
379 int binCons = getBin(coc);
380 binDOF[binCons] -= 1;
385 for (
unsigned int i = 0; i < nBins_; i++) {
386 RealType temp(0.0), ePot(0.0);
388 RealType z, den(0.0), binVolume(0.0), dz(0.0);
389 std::vector<RealType> nden(outputTypes_.size(), 0.0);
391 z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat_(axis_, axis_);
392 data_[Z].accumulator[i]->add(z);
394 binVolume = volume_ / nBins_;
395 dz = hmat_(axis_, axis_) / (RealType)nBins_;
399 if (outputMask_[ELECTRICFIELD] && binEFieldCount[i] > 0) {
400 eField = binEField[i] / RealType(binEFieldCount[i]);
401 data_[ELECTRICFIELD].accumulator[i]->add(eField);
404 if (outputMask_[ELECTROSTATICPOTENTIAL] && binEFieldCount[i] > 0) {
405 ePot += eField[axis_] * dz;
406 data_[ELECTROSTATICPOTENTIAL].accumulator[i]->add(ePot);
411 if (outputMask_[DENSITY]) {
412 den = binMass[i] * Constants::densityConvert / binVolume;
413 data_[DENSITY].accumulator[i]->add(den);
416 if (outputMask_[ACTIVITY]) {
417 for (
unsigned int j = 0; j < outputTypes_.size(); j++) {
418 nden[j] = (binTypeCounts[i][j] / binVolume) *
419 Constants::concentrationConvert;
421 data_[ACTIVITY].accumulator[i]->add(nden);
424 if (binCount[i] > 0) {
427 if (outputMask_[VELOCITY]) {
428 vel = binP[i] / binMass[i];
429 data_[VELOCITY].accumulator[i]->add(vel);
432 if (outputMask_[TEMPERATURE]) {
434 temp = 2.0 * binKE[i] /
435 (binDOF[i] * Constants::kb * Constants::energyConvert);
436 data_[TEMPERATURE].accumulator[i]->add(temp);
438 std::cerr <<
"No degrees of freedom in this bin?\n";
442 if (outputMask_[CHARGE])
443 data_[CHARGE].accumulator[i]->add(binCharge[i]);
445 if (outputMask_[CHARGEVELOCITY])
446 data_[CHARGEVELOCITY].accumulator[i]->add(binChargeVelocity[i]);
451 RNEMDR::RNEMDR(
SimInfo* info,
const std::string& filename,
452 const std::string& sele,
const std::string& comsele,
453 int nrbins, RealType binWidth) :
455 setOutputName(
getPrefix(filename) +
".rnemdR");
458 data_.resize(RNEMDR::ENDINDEX);
461 r.units =
"Angstroms";
463 r.dataHandling = DataHandling::Average;
464 for (
int i = 0; i < nBins_; i++)
465 r.accumulator.push_back(
466 std::make_unique<AccumulatorView<RealAccumulator>>());
467 data_[R] = std::move(r);
469 OutputData temperature;
470 temperature.units =
"K";
471 temperature.title =
"Temperature";
472 temperature.dataHandling = DataHandling::Average;
473 for (
unsigned int i = 0; i < nBins_; i++)
474 temperature.accumulator.push_back(
475 std::make_unique<AccumulatorView<RealAccumulator>>());
476 data_[TEMPERATURE] = std::move(temperature);
478 OutputData angularVelocity;
479 angularVelocity.units =
"angstroms/fs";
480 angularVelocity.title =
"Velocity";
481 angularVelocity.dataHandling = DataHandling::Average;
482 for (
unsigned int i = 0; i < nBins_; i++)
483 angularVelocity.accumulator.push_back(
484 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
485 data_[ANGULARVELOCITY] = std::move(angularVelocity);
488 density.units =
"g cm^-3";
489 density.title =
"Density";
490 density.dataHandling = DataHandling::Average;
491 for (
unsigned int i = 0; i < nBins_; i++)
492 density.accumulator.push_back(
493 std::make_unique<AccumulatorView<RealAccumulator>>());
494 data_[DENSITY] = std::move(density);
497 outputMask_.set(TEMPERATURE);
498 outputMask_.set(ANGULARVELOCITY);
499 outputMask_.set(DENSITY);
502 void RNEMDR::processFrame(
int istep) {
503 ShellStatistics::processFrame(istep);
506 seleMan_.setSelectionSet(evaluator_.evaluate());
519 std::vector<int> binCount(nBins_, 0);
520 std::vector<RealType> binMass(nBins_, 0.0);
521 std::vector<Vector3d> binP(nBins_, V3Zero);
522 std::vector<RealType> binOmega(nBins_, 0.0);
523 std::vector<Vector3d> binL(nBins_, V3Zero);
524 std::vector<Mat3x3d> binI(nBins_);
525 std::vector<RealType> binKE(nBins_, 0.0);
526 std::vector<int> binDOF(nBins_, 0);
528 SimInfo::MoleculeIterator miter;
529 std::vector<StuntDouble*>::iterator iiter;
530 std::vector<AtomType*>::iterator at;
534 Molecule::ConstraintPairIterator cpi;
539 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
540 sd = mol->nextIntegrableObject(iiter)) {
541 if (seleMan_.isSelected(sd)) {
543 binNo = getBin(sd->
getPos());
545 if (binNo >= 0 && binNo <
int(nBins_)) {
548 rPos = sd->
getPos() - coordinateOrigin_;
559 KE += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
560 angMom[k] * angMom[k] / Ia(k, k));
563 KE += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
564 angMom[1] * angMom[1] / Ia(1, 1) +
565 angMom[2] * angMom[2] / Ia(2, 2));
570 L = mass *
cross(rPos, vel);
571 I = outProduct(rPos, rPos) * mass;
572 r2 = rPos.lengthSquare();
573 I(0, 0) += mass * r2;
574 I(1, 1) += mass * r2;
575 I(2, 2) += mass * r2;
578 binMass[binNo] += mass;
579 binP[binNo] += mass * vel;
583 binDOF[binNo] += DOF;
587 if (seleMan_.isSelected(mol)) {
588 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
589 consPair = mol->nextConstraintPair(cpi)) {
590 Vector3d posA = consPair->getConsElem1()->getPos();
591 Vector3d posB = consPair->getConsElem2()->getPos();
594 int binCons = getBin(coc);
595 if (binCons >= 0 && binCons <
int(nBins_)) { binDOF[binCons] -= 1; }
600 for (
unsigned int i = 0; i < nBins_; i++) {
601 RealType r, rinner, router, den(0.0), binVolume(0.0), temp(0.0);
604 r = (((RealType)i + 0.5) * binWidth_);
605 rinner = (RealType)i * binWidth_;
606 router = (RealType)(i + 1) * binWidth_;
608 (4.0 * Constants::PI * (pow(router, 3) - pow(rinner, 3))) / 3.0;
610 data_[R].accumulator[i]->add(r);
614 den = binMass[i] * Constants::densityConvert / binVolume;
615 data_[DENSITY].accumulator[i]->add(den);
620 omega = binI[i].inverse() * binL[i];
621 data_[ANGULARVELOCITY].accumulator[i]->add(omega);
623 temp = 2.0 * binKE[i] /
624 (binDOF[i] * Constants::kb * Constants::energyConvert);
625 data_[TEMPERATURE].accumulator[i]->add(temp);
630 RNEMDRTheta::RNEMDRTheta(
SimInfo* info,
const std::string& filename,
631 const std::string& sele,
const std::string& comsele,
632 int nrbins, RealType binWidth,
int nangleBins) :
634 nAngleBins_(nangleBins) {
635 Globals* simParams = info->getSimParams();
637 bool hasAngularMomentumFluxVector =
638 rnemdParams->haveAngularMomentumFluxVector();
640 if (hasAngularMomentumFluxVector) {
641 std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
643 if (amf.size() != 3) {
644 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
645 "RNEMDRTheta: Incorrect number of parameters specified for "
646 "angularMomentumFluxVector.\n"
647 "\tthere should be 3 parameters, but %lu were specified.\n",
649 painCave.isFatal = 1;
652 fluxVector_.x() = amf[0];
653 fluxVector_.y() = amf[1];
654 fluxVector_.z() = amf[2];
656 std::string fluxStr = rnemdParams->getFluxType();
658 if (fluxStr.find(
"Lx") != std::string::npos) {
660 }
else if (fluxStr.find(
"Ly") != std::string::npos) {
667 fluxVector_.normalize();
669 setOutputName(
getPrefix(filename) +
".rnemdRTheta");
672 r_.units =
"Angstroms";
674 for (
int i = 0; i < nBins_; i++)
675 r_.accumulator.push_back(
676 std::make_unique<AccumulatorView<RealAccumulator>>());
678 angularVelocity_.units =
"1/fs";
679 angularVelocity_.title =
"Projected Angular Velocity";
680 for (
unsigned int i = 0; i < nBins_; i++) {
681 angularVelocity_.accumulator.push_back(
682 std::make_unique<AccumulatorView<StdVectorAccumulator>>());
686 std::pair<int, int> RNEMDRTheta::getBins(
Vector3d pos) {
687 std::pair<int, int> result;
689 Vector3d rPos = pos - coordinateOrigin_;
690 RealType cosAngle =
dot(rPos, fluxVector_) / rPos.length();
692 result.first = int(rPos.length() / binWidth_);
693 result.second = int((nAngleBins_)*0.5 * (cosAngle + 1.0));
697 void RNEMDRTheta::processFrame(
int istep) {
698 ShellStatistics::processFrame(istep);
701 seleMan_.setSelectionSet(evaluator_.evaluate());
707 std::vector<std::vector<int>> binCount(nBins_);
708 std::vector<std::vector<Mat3x3d>> binI(nBins_);
709 std::vector<std::vector<Vector3d>> binL(nBins_);
711 for (std::size_t i {}; i < nBins_; ++i) {
712 binCount[i].resize(nAngleBins_);
713 binI[i].resize(nAngleBins_);
714 binL[i].resize(nAngleBins_);
721 std::pair<int, int> bins = getBins(sd->
getPos());
723 if (bins.first >= 0 && bins.first <
int(nBins_)) {
724 if (bins.second >= 0 && bins.second < nAngleBins_) {
730 I = outProduct(rPos, rPos) * m;
731 RealType r2 = rPos.lengthSquare();
736 binI[bins.first][bins.second] += I;
737 binL[bins.first][bins.second] += L;
738 binCount[bins.first][bins.second]++;
743 for (
unsigned int i = 0; i < nBins_; i++) {
744 RealType r = (((RealType)i + 0.5) * binWidth_);
745 r_.accumulator[i]->add(r);
747 std::vector<RealType> projections(nAngleBins_);
749 for (
int j = 0; j < nAngleBins_; j++) {
752 if (binCount[i][j] > 0) { omega = binI[i][j].inverse() * binL[i][j]; }
755 projections[j] =
dot(omega, fluxVector_);
758 angularVelocity_.accumulator[i]->add(projections);
762 void RNEMDRTheta::writeOutput() {
763 std::ofstream outStream(outputFilename_.c_str());
765 if (outStream.is_open()) {
767 outStream <<
"# SPATIAL STATISTICS\n";
768 outStream <<
"#nBins = " << nBins_ <<
"\t binWidth = " << binWidth_
769 <<
" maxR = " << nBins_ * binWidth_ <<
"\n";
770 outStream <<
"#fluxVector = " << fluxVector_ <<
"\tBins = " << nAngleBins_
772 outStream <<
"#\t" << angularVelocity_.title <<
"("
773 << angularVelocity_.units <<
")\t\t";
775 outStream << std::endl;
777 outStream.precision(8);
779 for (
unsigned int i = 0; i < nBins_; i++) {
780 std::size_t n {r_.accumulator[i]->getCount()};
783 std::string message =
784 "StaticAnalyser detected a numerical error writing: " +
785 angularVelocity_.title +
" for bin " + std::to_string(i);
787 angularVelocity_.accumulator[i]->writeData(outStream, message);
790 outStream << std::endl;
AtomType is what OpenMD looks to for unchanging data about an atom.
bool isDynamic()
Tests if the result from evaluation of script is dynamic.
StuntDouble * nextSelected(int &i)
Finds the next selected StuntDouble in the selection.
StuntDouble * beginSelected(int &i)
Finds the first selected StuntDouble in the selection.
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
void wrapVector(Vector3d &v)
Wrapping the vector according to periodic boundary condition.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
Vector3d getVel()
Returns the current velocity of this stuntDouble.
int linearAxis()
Returns the linear axis of the rigidbody, atom and directional atom will always return -1.
RealType getMass()
Returns the mass of this stuntDouble.
virtual Mat3x3d getI()=0
Returns the inertia tensor of this stuntDouble.
bool isLinear()
Tests the if this stuntDouble is a linear rigidbody.
Vector3d getPos()
Returns the current position of this stuntDouble.
RealType getFlucQPos()
Returns the current fluctuating charge of this stuntDouble.
Vector3d getElectricField()
Returns the current electric field of this stuntDouble.
bool isRigidBody()
Tests if this stuntDouble is a rigid body.
Vector3d getJ()
Returns the current angular momentum of this stuntDouble (body -fixed).
bool isAtom()
Tests if this stuntDouble is an atom.
bool isDirectional()
Tests if this stuntDouble is a directional one.
RealType getFlucQVel()
Returns the current charge velocity of this stuntDouble.
Real lengthSquare()
Returns the squared length of this vector.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
std::string getPrefix(const std::string &str)