45#include "applications/staticProps/TranslationalOrderParamR.hpp"
53#include "utils/Revision.hpp"
54#include "utils/simError.h"
56#define HONKING_LARGE_VALUE 1.0e10
60 TranslationalOrderParamR::TranslationalOrderParamR(
61 SimInfo* info,
const std::string& filename,
const std::string& sele1,
62 const std::string& sele2,
const std::string& sele3, RealType rCut,
63 RealType len,
int nrbins) :
64 StaticAnalyser(info, filename, nrbins),
65 selectionScript1_(sele1), selectionScript2_(sele2),
66 selectionScript3_(sele3), seleMan1_(info), seleMan2_(info),
67 seleMan3_(info), evaluator1_(info), evaluator2_(info), evaluator3_(info),
68 len_(len), nBins_(nrbins) {
69 setAnalysisType(
"Translational Order Parameter(r)");
71 evaluator1_.loadScriptString(sele1);
72 if (!evaluator1_.isDynamic()) {
73 seleMan1_.setSelectionSet(evaluator1_.evaluate());
74 validateSelection1(seleMan1_);
77 evaluator2_.loadScriptString(sele2);
78 if (!evaluator2_.isDynamic()) {
79 seleMan2_.setSelectionSet(evaluator2_.evaluate());
80 validateSelection2(seleMan2_);
83 if (!evaluator1_.isDynamic() && !evaluator2_.isDynamic()) {
86 common_ = seleMan1_ & seleMan2_;
87 sele1_minus_common_ = seleMan1_ - common_;
88 sele2_minus_common_ = seleMan2_ - common_;
90 nSelected1_ = seleMan1_.getSelectionCount();
91 nSelected2_ = seleMan2_.getSelectionCount();
92 int nIntersect = common_.getSelectionCount();
94 nPairs_ = nSelected1_ * nSelected2_ - (nIntersect + 1) * nIntersect / 2;
97 evaluator3_.loadScriptString(sele3);
98 if (!evaluator3_.isDynamic()) {
99 seleMan3_.setSelectionSet(evaluator3_.evaluate());
102 deltaR_ = len_ / nBins_;
105 sliceT_.resize(nBins_);
106 sliceCount_.resize(nBins_);
107 std::fill(sliceT_.begin(), sliceQ_.end(), 0.0);
108 std::fill(sliceCount_.begin(), sliceCount_.end(), 0);
110 setOutputName(
getPrefix(filename) +
".Tr");
113 void TranslationalOrderParamR::process() {
121 Vector3d ri, rj, rk, rik, rkj;
125 std::vector<std::pair<RealType, StuntDouble*>> myNeighbors;
129 bool usePeriodicBoundaryConditions_ =
130 info_->getSimParams()->getUsePeriodicBoundaryConditions();
132 DumpReader reader(info_, dumpFilename_);
133 int nFrames = reader.getNFrames();
135 for (
int istep = 0; istep < nFrames; istep += step_) {
136 reader.readFrame(istep);
137 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
139 if (evaluator1_.isDynamic()) {
140 seleMan1_.setSelectionSet(evaluator1_.evaluate());
141 validateSelection1(seleMan1_);
144 if (evaluator2_.isDynamic()) {
145 seleMan2_.setSelectionSet(evaluator2_.evaluate());
146 validateSelection2(seleMan2_);
149 if (evaluator3_.isDynamic()) {
150 seleMan3_.setSelectionSet(evaluator3_.evaluate());
153 if (evaluator1_.isDynamic() || evaluator2_.isDynamic()) {
154 common_ = seleMan1_ & seleMan2_;
155 sele1_minus_common_ = seleMan1_ - common_;
156 sele2_minus_common_ = seleMan2_ - common_;
157 nSelected1_ = seleMan1_.getSelectionCount();
158 nSelected2_ = seleMan2_.getSelectionCount();
159 int nIntersect = common_.getSelectionCount();
161 nPairs_ = nSelected1_ * nSelected2_ - (nIntersect + 1) * nIntersect / 2;
164 processNonOverlapping(sele1_minus_common_, seleMan2_);
165 processNonOverlapping(common_, sele2_minus_common_);
166 processOverlapping(common_);
176 void TranslationalOrderParamR::processNonOverlapping(SelectionManager& sman1,
177 SelectionManager& sman2) {
188 for (sd1 = sman1.beginSelected(i); sd1 != NULL;
189 sd1 = sman1.nextSelected(i)) {
190 for (sd2 = sman2.beginSelected(j); sd2 != NULL;
191 sd2 = sman2.nextSelected(j)) {
192 collectHistogram(sd1, sd2);
197 void TranslationalOrderParamR::processOverlapping(SelectionManager& sman) {
208 for (sd1 = sman.beginSelected(i); sd1 != NULL; sd1 = sman.nextSelected(i)) {
209 for (j = i, sd2 = sman.nextSelected(j); sd2 != NULL;
210 sd2 = sman.nextSelected(j)) {
211 collectHistogram(sd1, sd2);
217 void TranslationalOrderParamR::collectHistogram(StuntDouble* sd1,
219 if (sd1 == sd2) {
return; }
221 bool usePeriodicBoundaryConditions_ =
222 info_->getSimParams()->getUsePeriodicBoundaryConditions();
224 Vector3d pos1 = sd1->getPos();
225 Vector3d pos2 = sd2->getPos();
226 Vector3d r12 = pos2 - pos1;
227 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(r12);
231 if (distance < len_) {
232 int whichBin = int(distance / deltaR_);
233 histogram_[whichBin] += 2;
239 for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
240 sd = seleMan1_.nextSelected(isd1)) {
241 myIndex = sd->getGlobalIndex();
243 std::fill(histogram_.begin(), histogram_.end(), 0);
245 for (sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
246 sd2 = seleMan2_.nextSelected(isd2)) {
247 if (sd2->getGlobalIndex() != myIndex) {
248 vec = sd->getPos() - sd2->getPos();
250 if (usePeriodicBoundaryConditions_)
251 currentSnapshot_->wrapVector(vec);
256 int whichBin = int(distance / deltaR_);
257 histogram_[whichBin] += 2;
263 for (
int i = 0; i < nbors - 1; i++) {
264 sdi = myNeighbors[i].second;
267 if (usePeriodicBoundaryConditions_) currentSnapshot_->wrapVector(rik);
271 for (
int j = i + 1; j < nbors; j++) {
272 sdj = myNeighbors[j].second;
275 if (usePeriodicBoundaryConditions_)
276 currentSnapshot_->wrapVector(rkj);
279 cospsi =
dot(rik, rkj);
283 Qk -= (pow(cospsi + 1.0 / 3.0, 2) * 2.25 / nang);
288 RealType shortest = HONKING_LARGE_VALUE;
291 for (sd3 = seleMan3_.beginSelected(isd3); sd3 != NULL;
292 sd3 = seleMan3_.nextSelected(isd3)) {
293 vec = sd->getPos() - sd3->getPos();
295 if (usePeriodicBoundaryConditions_)
296 currentSnapshot_->wrapVector(vec);
300 if (r < shortest) shortest = r;
303 int whichBin = int(shortest / deltaR_);
304 if (whichBin <
int(nBins_)) {
305 sliceT_[whichBin] += Qk;
306 sliceCount_[whichBin] += 1;
314 void TranslationalOrderParamR::writeTr() {
316 std::ofstream tRstream(outputFilename_.c_str());
317 if (tRstream.is_open()) {
318 tRstream <<
"# " << getAnalysisType() <<
"\n";
319 tRstream <<
"# OpenMD " << rev.getFullRevision() <<
"\n";
320 tRstream <<
"# " << rev.getBuildDate() <<
"\n";
321 tRstream <<
"#selection 1: (" << selectionScript1_ <<
")\n";
322 tRstream <<
"#selection 2: (" << selectionScript2_ <<
")\n";
323 tRstream <<
"#selection 3: (" << selectionScript3_ <<
")\n";
324 if (!paramString_.empty())
325 tRstream <<
"# parameters: " << paramString_ <<
"\n";
327 tRstream <<
"#distance"
329 for (
unsigned int i = 0; i < sliceT_.size(); ++i) {
330 RealType Rval = (i + 0.5) * deltaR_;
331 if (sliceCount_[i] != 0) {
332 tRstream << Rval <<
"\t" << sliceT_[i] / (RealType)sliceCount_[i]
338 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
339 "TranslationalOrderParamR: unable to open %s\n",
340 outputFilename_.c_str());
341 painCave.isFatal = 1;
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)
Real distance(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the distance between two DynamicVectors.