45#include "applications/dynamicProps/VelAngularVelOutProdCorrFunc.hpp"
50 VelAngularVelOutProdCorrFunc::VelAngularVelOutProdCorrFunc(
51 SimInfo* info,
const std::string& filename,
const std::string& sele1,
52 const std::string& sele2) :
53 ObjectCCF<Mat3x3d>(info, filename, sele1, sele2) {
55 "Velocity - Angular Velocity Outer Product Correlation Function");
56 setOutputName(
getPrefix(dumpFilename_) +
".vwOutProdcorr");
58 velocity_.resize(nFrames_);
59 angularVelocity_.resize(nFrames_);
61 sumVelocity_ = V3Zero;
62 sumAngularVelocity_ = V3Zero;
65 angularVelocityCount_ = 0;
67 propertyTemp = V3Zero;
70 int VelAngularVelOutProdCorrFunc::computeProperty1(
int frame,
72 Vector3d vel = sd->getVel();
75 velocity_[frame].push_back(propertyTemp);
76 sumVelocity_ += propertyTemp;
78 return velocity_[frame].size() - 1;
81 int VelAngularVelOutProdCorrFunc::computeProperty2(
int frame,
83 if (sd->isDirectional()) {
84 Mat3x3d momentInertia = sd->getI();
85 Vector3d angMom = sd->getJ();
86 Vector3d omega = momentInertia.inverse() * angMom;
90 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
91 "The selection contains non-directional entities. Your selection should include\
92 Directional atoms and/or Rigid Bodies.\n");
97 angularVelocity_[frame].push_back(propertyTemp);
98 sumAngularVelocity_ += propertyTemp;
99 angularVelocityCount_++;
100 return angularVelocity_[frame].size() - 1;
103 Mat3x3d VelAngularVelOutProdCorrFunc::calcCorrVal(
int frame1,
int frame2,
107 outProduct(velocity_[frame1][id1], angularVelocity_[frame2][id2]);
110 outProduct(velocity_[frame2][id2], angularVelocity_[frame1][id1]);
112 tmpMat_3 = 0.5 * (tmpMat_1 + tmpMat_2);
116 void VelAngularVelOutProdCorrFunc::postCorrelate() {
118 sumVelocity_ /= RealType(velocityCount_);
121 sumAngularVelocity_ /= RealType(angularVelocityCount_);
123 Mat3x3d correlationOfAverages_ =
124 outProduct(sumVelocity_, sumAngularVelocity_);
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)