47#include "applications/staticProps/P2R.hpp"
54#include "brains/Thermo.hpp"
57#include "utils/Revision.hpp"
58#include "utils/simError.h"
62 P2R::P2R(SimInfo* info,
const std::string& filename,
const std::string& sele1,
64 StaticAnalyser(info, filename, nbins),
65 doVect_(true), doOffset_(false), selectionScript1_(sele1),
66 seleMan1_(info), seleMan2_(info), evaluator1_(info), evaluator2_(info) {
67 evaluator1_.loadScriptString(sele1);
68 if (!evaluator1_.isDynamic()) {
69 seleMan1_.setSelectionSet(evaluator1_.evaluate());
73 "2nd order Legendre Polynomial Correlation using r as reference axis");
74 setOutputName(
getPrefix(filename) +
".P2R");
75 std::stringstream params;
76 const std::string paramString = params.str();
77 setParameterString(paramString);
80 P2R::P2R(SimInfo* info,
const std::string& filename,
const std::string& sele1,
81 const std::string& sele2,
unsigned int nbins) :
82 StaticAnalyser(info, filename, nbins),
83 doVect_(false), doOffset_(false), selectionScript1_(sele1),
84 selectionScript2_(sele2), seleMan1_(info), seleMan2_(info),
85 evaluator1_(info), evaluator2_(info) {
86 setOutputName(
getPrefix(filename) +
".P2R");
88 evaluator1_.loadScriptString(sele1);
89 if (!evaluator1_.isDynamic()) {
90 seleMan1_.setSelectionSet(evaluator1_.evaluate());
93 evaluator2_.loadScriptString(sele2);
94 if (!evaluator2_.isDynamic()) {
95 seleMan2_.setSelectionSet(evaluator2_.evaluate());
99 "2nd order Legendre Polynomial Correlation using r as reference axis");
100 setOutputName(
getPrefix(filename) +
".P2R");
101 std::stringstream params;
102 const std::string paramString = params.str();
103 setParameterString(paramString);
106 P2R::P2R(SimInfo* info,
const std::string& filename,
const std::string& sele1,
107 int seleOffset,
unsigned int nbins) :
108 StaticAnalyser(info, filename, nbins),
109 doVect_(false), doOffset_(true), doOffset2_(false),
110 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
111 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset) {
112 setOutputName(
getPrefix(filename) +
".P2R");
114 evaluator1_.loadScriptString(sele1);
115 if (!evaluator1_.isDynamic()) {
116 seleMan1_.setSelectionSet(evaluator1_.evaluate());
120 "2nd order Legendre Polynomial Correlation using r as reference axis");
121 setOutputName(
getPrefix(filename) +
".P2R");
122 std::stringstream params;
123 const std::string paramString = params.str();
124 setParameterString(paramString);
127 P2R::P2R(SimInfo* info,
const std::string& filename,
const std::string& sele1,
128 int seleOffset,
int seleOffset2,
unsigned int nbins) :
129 StaticAnalyser(info, filename, nbins),
130 doVect_(false), doOffset_(true), doOffset2_(true),
131 selectionScript1_(sele1), seleMan1_(info), seleMan2_(info),
132 evaluator1_(info), evaluator2_(info), seleOffset_(seleOffset),
133 seleOffset2_(seleOffset2) {
134 setOutputName(
getPrefix(filename) +
".P2R");
136 evaluator1_.loadScriptString(sele1);
137 if (!evaluator1_.isDynamic()) {
138 seleMan1_.setSelectionSet(evaluator1_.evaluate());
142 "2nd order Legendre Polynomial Correlation using r as reference axis");
143 setOutputName(
getPrefix(filename) +
".P2R");
144 std::stringstream params;
145 const std::string paramString = params.str();
146 setParameterString(paramString);
149 void P2R::process() {
154 bool usePeriodicBoundaryConditions_ =
155 info_->getSimParams()->getUsePeriodicBoundaryConditions();
157 Thermo thermo(info_);
158 DumpReader reader(info_, dumpFilename_);
159 int nFrames = reader.getNFrames();
161 nProcessed_ = nFrames / step_;
163 for (
int istep = 0; istep < nFrames; istep += step_) {
164 reader.readFrame(istep);
165 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
167 Vector3d CenterOfMass = thermo.getCom();
169 if (evaluator1_.isDynamic()) {
170 seleMan1_.setSelectionSet(evaluator1_.evaluate());
174 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
175 sd1 = seleMan1_.nextSelected(ii)) {
176 Vector3d pos = sd1->getPos();
178 Vector3d r1 = CenterOfMass - pos;
181 if (sd1->isDirectional()) {
182 Vector3d vec = sd1->getA().transpose() * V3Z;
185 RealType cosangle =
dot(r1, vec);
187 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
193 for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
194 sd1 = seleMan1_.nextSelected(ii)) {
201 int sd1Aind = sd1->getGlobalIndex() + seleOffset2_;
202 StuntDouble* sd1A = info_->getIOIndexToIntegrableObject(sd1Aind);
203 r1 = CenterOfMass - sd1A->getPos();
205 r1 = CenterOfMass - sd1->getPos();
208 if (usePeriodicBoundaryConditions_)
209 currentSnapshot_->wrapVector(r1);
211 int sd2Index = sd1->getGlobalIndex() + seleOffset_;
212 sd2 = info_->getIOIndexToIntegrableObject(sd2Index);
214 Vector3d r2 = CenterOfMass - sd2->getPos();
215 if (usePeriodicBoundaryConditions_)
216 currentSnapshot_->wrapVector(r1);
218 Vector3d rc = 0.5 * (r1 + r2);
219 if (usePeriodicBoundaryConditions_)
220 currentSnapshot_->wrapVector(rc);
222 Vector3d vec = r1 - r2;
223 if (usePeriodicBoundaryConditions_)
224 currentSnapshot_->wrapVector(vec);
228 RealType cosangle =
dot(rc, vec);
229 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
233 if (evaluator2_.isDynamic()) {
234 seleMan2_.setSelectionSet(evaluator2_.evaluate());
237 if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount()) {
238 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
239 "In frame %d, the number of selected StuntDoubles are\n"
240 "\tnot the same in --sele1 and sele2\n",
242 painCave.severity = OPENMD_INFO;
243 painCave.isFatal = 0;
247 for (sd1 = seleMan1_.beginSelected(ii),
248 sd2 = seleMan2_.beginSelected(jj);
249 sd1 != NULL && sd2 != NULL; sd1 = seleMan1_.nextSelected(ii),
250 sd2 = seleMan2_.nextSelected(jj)) {
251 Vector3d r1 = CenterOfMass - sd1->getPos();
252 if (usePeriodicBoundaryConditions_)
253 currentSnapshot_->wrapVector(r1);
255 Vector3d r2 = CenterOfMass - sd2->getPos();
256 if (usePeriodicBoundaryConditions_)
257 currentSnapshot_->wrapVector(r1);
259 Vector3d rc = 0.5 * (r1 + r2);
260 if (usePeriodicBoundaryConditions_)
261 currentSnapshot_->wrapVector(rc);
263 Vector3d vec = r1 - r2;
264 if (usePeriodicBoundaryConditions_)
265 currentSnapshot_->wrapVector(vec);
269 RealType cosangle =
dot(rc, vec);
270 P2_ += 0.5 * (3.0 * cosangle * cosangle - 1.0);
280 void P2R::processHistogram() {
281 if (count_ > 0) P2_ /= count_;
284 void P2R::writeP2R() {
285 std::ofstream ofs(outputFilename_.c_str());
289 ofs <<
"# " << getAnalysisType() <<
"\n";
290 ofs <<
"# OpenMD " << rev.getFullRevision() <<
"\n";
291 ofs <<
"# " << rev.getBuildDate() <<
"\n";
292 ofs <<
"#nFrames:\t" << nProcessed_ <<
"\n";
293 ofs <<
"#selection1: (" << selectionScript1_ <<
")";
294 if (!doVect_) { ofs <<
"\tselection2: (" << selectionScript2_ <<
")"; }
296 if (!paramString_.empty())
297 ofs <<
"# parameters: " << paramString_ <<
"\n";
299 ofs <<
"#P2 correlation:\n";
303 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
304 "P2R: unable to open %s\n", outputFilename_.c_str());
305 painCave.isFatal = 1;
void normalize()
Normalizes this vector in place.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
std::string getPrefix(const std::string &str)