45#include "applications/dynamicProps/AngularVelVelOutProdCorrFunc.hpp"
50 AngularVelVelOutProdCorrFunc::AngularVelVelOutProdCorrFunc(
51 SimInfo* info,
const std::string& filename,
const std::string& sele1,
52 const std::string& sele2) :
53 ObjectCCF<Mat3x3d>(info, filename, sele1, sele2) {
55 "Angular Velocity - Velocity Outer Product Correlation Function");
56 setOutputName(
getPrefix(dumpFilename_) +
".wvOutProdcorr");
58 velocity_.resize(nFrames_);
59 angularVelocity_.resize(nFrames_);
61 sumVelocity_ = V3Zero;
62 sumAngularVelocity_ = V3Zero;
65 angularVelocityCount_ = 0;
67 propertyTemp = V3Zero;
70 int AngularVelVelOutProdCorrFunc::computeProperty1(
int frame,
72 if (sd->isDirectional()) {
73 Mat3x3d momentInertia = sd->getI();
74 Vector3d angMom = sd->getJ();
75 Vector3d omega = momentInertia.inverse() * angMom;
79 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
80 "The selection contains non-directional entities. Your selection should include\
81 Directional atoms and/or Rigid Bodies.\n");
86 angularVelocity_[frame].push_back(propertyTemp);
87 sumAngularVelocity_ += propertyTemp;
88 angularVelocityCount_++;
89 return angularVelocity_[frame].size() - 1;
92 int AngularVelVelOutProdCorrFunc::computeProperty2(
int frame,
94 Vector3d vel = sd->getVel();
97 velocity_[frame].push_back(propertyTemp);
98 sumVelocity_ += propertyTemp;
100 return velocity_[frame].size() - 1;
103 Mat3x3d AngularVelVelOutProdCorrFunc::calcCorrVal(
int frame1,
int frame2,
107 outProduct(angularVelocity_[frame1][id1], velocity_[frame2][id2]);
110 outProduct(angularVelocity_[frame2][id2], velocity_[frame1][id1]);
112 tmpMat_3 = 0.5 * (tmpMat_1 + tmpMat_2);
116 void AngularVelVelOutProdCorrFunc::postCorrelate() {
118 sumAngularVelocity_ /= RealType(angularVelocityCount_);
121 sumVelocity_ /= RealType(velocityCount_);
123 Mat3x3d correlationOfAverages_ =
124 outProduct(sumAngularVelocity_, sumVelocity_);
126 for (
unsigned int i = 0; i < nTimeBins_; ++i) {
128 histogram_[i] /= RealType(count_[i]);
132 histogram_[i] -= correlationOfAverages_;
134 histogram_[i] = M3Zero;
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)