OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
TimeCorrFunc.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the following paper when you publish your work:
33 *
34 * [1] Drisko et al., J. Open Source Softw. 9, 7004 (2024).
35 *
36 * Good starting points for code and simulation methodology are:
37 *
38 * [2] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
39 * [3] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
40 * [4] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
41 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
42 * [6] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
43 * [7] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
44 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
45 * [9] Drisko & Gezelter, J. Chem. Theory Comput. 20, 4986-4997 (2024).
46 */
47
48#include "applications/dynamicProps/TimeCorrFunc.hpp"
49
50#include <memory>
51
52#include "math/DynamicVector.hpp"
53#include "primitives/Molecule.hpp"
54#include "utils/Revision.hpp"
55#include "utils/simError.h"
56
57using namespace std;
58namespace OpenMD {
59
60 template<typename T>
61 TimeCorrFunc<T>::TimeCorrFunc(SimInfo* info, const string& filename,
62 const string& sele1, const string& sele2) :
63 info_(info),
64 dumpFilename_(filename), seleMan1_(info_), seleMan2_(info_),
65 selectionScript1_(sele1), selectionScript2_(sele2), evaluator1_(info_),
66 evaluator2_(info_), autoCorrFunc_(false), doSystemProperties_(false),
67 doMolecularProperties_(false), doObjectProperties_(false),
68 doBondProperties_(false), allowTimeFuzz_(false) {
69 reader_ = new DumpReader(info_, dumpFilename_);
70
71 uniqueSelections_ = (sele1.compare(sele2) != 0) ? true : false;
72
73 Globals* simParams = info_->getSimParams();
74 if (simParams->haveSampleTime()) {
75 deltaTime_ = simParams->getSampleTime();
76 } else {
77 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
78 "TimeCorrFunc Error: can not figure out deltaTime\n");
79 painCave.isFatal = 1;
80 simError();
81 }
82
83 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "Scanning for frames.");
84 painCave.isFatal = 0;
85 painCave.severity = OPENMD_INFO;
86 simError();
87
88 nFrames_ = reader_->getNFrames();
89
90 if (nFrames_ < 2) {
91 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
92 "Not enough frames for a meaningful time correlation"
93 " (need >= 2, have %d).\n",
94 nFrames_);
95 painCave.isFatal = 0;
96 painCave.severity = OPENMD_WARNING;
97 simError();
98 }
99
100 nTimeBins_ = nFrames_;
101
102 T zeroType(0.0);
103 histogram_.resize(nTimeBins_, zeroType);
104 count_.resize(nTimeBins_, 0);
105
106 times_.resize(nFrames_);
107 sele1ToIndex_.resize(nFrames_);
108 GIDtoSele1_.resize(nFrames_);
109 selection1StartFrame_.resize(nFrames_);
110
111 if (uniqueSelections_) {
112 sele2ToIndex_.resize(nFrames_);
113 GIDtoSele2_.resize(nFrames_);
114 selection2StartFrame_.resize(nFrames_);
115 }
116
117 progressBar_ = std::make_unique<ProgressBar>();
118 }
119
120 template<typename T>
121 void TimeCorrFunc<T>::preCorrelate() {
122 evaluator1_.loadScriptString(selectionScript1_);
123 // if selection is static, we only need to evaluate it once
124 if (!evaluator1_.isDynamic()) {
125 seleMan1_.setSelectionSet(evaluator1_.evaluate());
126 validateSelection(seleMan1_);
127 }
128
129 if (uniqueSelections_) {
130 evaluator2_.loadScriptString(selectionScript2_);
131 if (!evaluator2_.isDynamic()) {
132 seleMan2_.setSelectionSet(evaluator2_.evaluate());
133 validateSelection(seleMan2_);
134 }
135 }
136
137 progressBar_->clear();
138
139 for (int istep = 0; istep < nFrames_; istep++) {
140 reader_->readFrame(istep);
141 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
142 times_[istep] = currentSnapshot_->getTime();
143
144 progressBar_->setStatus(istep + 1, nFrames_);
145 progressBar_->update();
146
147 computeFrame(istep);
148 }
149
150 dtMean_ = 0.0;
151 RealType dt2Avg(0.0);
152 dtSigma_ = 0.0;
153
154 if (nFrames_ > 1) {
155 for (int istep = 1; istep < nFrames_; istep++) {
156 RealType dt = times_[istep] - times_[istep - 1];
157 dtMean_ += dt;
158 dt2Avg += dt * dt;
159 }
160 dtMean_ /= RealType(nFrames_ - 1);
161 dt2Avg /= RealType(nFrames_ - 1);
162 dtSigma_ = std::sqrt(dt2Avg - dtMean_ * dtMean_);
163 }
164
165 if (dtSigma_ > 1.0e-6) {
166 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
167 "TimeCorrFunc::preCorrelate: sampleTime (%f) does not match\n"
168 "\tthe mean spacing between configurations (%f), with\n"
169 "\tsigma (%f). Proceeding with the mean value.\n",
170 deltaTime_, dtMean_, dtSigma_);
171 allowTimeFuzz_ = true;
172 painCave.isFatal = 0;
173 painCave.severity = OPENMD_INFO;
174 simError();
175 }
176
177 // Compute number of windows that fit in the trajectory
178 if (useWindowing_) {
179 // Need: nStart_ + navg_*nTimeBins_ + (navg_-1)*nSep_ <= nFrames_
180 // Solve: navg_ <= (nFrames_ - nStart_ - nSep_) / (nTimeBins_ + nSep_)
181 if (nStart_ >= nFrames_) {
182 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
183 "TimeCorrFunc: nStart (%d) >= nFrames (%d).\n",
184 nStart_, nFrames_);
185 painCave.isFatal = 1;
186 simError();
187 }
188 navg_ = (nFrames_ - nStart_ - nSep_) / nStride_;
189 if (navg_ < 1) {
190 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
191 "TimeCorrFunc: trajectory too short for even one window.\n"
192 "\tnFrames=%d, nStart=%d, nTimeBins=%d, nSep=%d\n",
193 nFrames_, nStart_, nTimeBins_, nSep_);
194 painCave.isFatal = 1;
195 simError();
196 }
197 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
198 "TimeCorrFunc: windowing enabled.\n"
199 "\tnStart=%d, nTimeBins=%d, nSep=%d, nStride=%d, navg=%d\n",
200 nStart_, nTimeBins_, nSep_, nStride_, navg_);
201 painCave.isFatal = 0;
202 painCave.severity = OPENMD_INFO;
203 simError();
204 } else {
205 // Original behavior: all frames are origins, nTimeBins_ = nFrames_
206 navg_ = nFrames_;
207 }
208 }
209
210 template<typename T>
211 void TimeCorrFunc<T>::computeFrame(int istep) {
212 Molecule* mol;
213 StuntDouble* sd;
214 Bond* bond;
215 int imol1, imol2, isd1, isd2, ibond1, ibond2;
216 unsigned int index;
217
218 if (evaluator1_.isDynamic()) {
219 seleMan1_.setSelectionSet(evaluator1_.evaluate());
220 validateSelection(seleMan1_);
221 }
222
223 if (uniqueSelections_ && evaluator2_.isDynamic()) {
224 seleMan2_.setSelectionSet(evaluator2_.evaluate());
225 validateSelection(seleMan2_);
226 }
227
228 if (doSystemProperties_) {
229 computeProperty1(istep);
230 if (!autoCorrFunc_) computeProperty2(istep);
231 }
232
233 if (doMolecularProperties_) {
234 // Map of molecular global IDs to selections
235 if (selectionModeRestart_)
236 GIDtoSele1_[istep].resize(info_->getNGlobalMolecules(), -1);
237
238 for (mol = seleMan1_.beginSelectedMolecule(imol1); mol != NULL;
239 mol = seleMan1_.nextSelectedMolecule(imol1)) {
240 index = computeProperty1(istep, mol);
241
242 if (index == sele1ToIndex_[istep].size()) {
243 sele1ToIndex_[istep].push_back(mol->getGlobalIndex());
244 } else {
245 sele1ToIndex_[istep].resize(index + 1);
246 sele1ToIndex_[istep][index] = mol->getGlobalIndex();
247 }
248 if (selectionModeRestart_) {
249 GIDtoSele1_[istep][mol->getGlobalIndex()] = index;
250
251 if (istep == 0) {
252 selection1StartFrame_[istep].push_back(istep);
253 } else {
254 int prevIndex = GIDtoSele1_[istep-1][mol->getGlobalIndex()];
255 if (prevIndex == -1)
256 selection1StartFrame_[istep].push_back(istep);
257 else
258 selection1StartFrame_[istep].push_back(selection1StartFrame_[istep-1][prevIndex]);
259 }
260 }
261
262 if (!uniqueSelections_) { index = computeProperty2(istep, mol); }
263 }
264
265 if (uniqueSelections_) {
266 // Map of molecular global IDs to selections
267 if (selectionModeRestart_)
268 GIDtoSele2_[istep].resize(info_->getNGlobalMolecules(), -1);
269
270 for (mol = seleMan2_.beginSelectedMolecule(imol2); mol != NULL;
271 mol = seleMan2_.nextSelectedMolecule(imol2)) {
272 if (autoCorrFunc_) {
273 index = computeProperty1(istep, mol);
274 } else {
275 index = computeProperty2(istep, mol);
276 }
277 if (index == sele2ToIndex_[istep].size()) {
278 sele2ToIndex_[istep].push_back(mol->getGlobalIndex());
279 } else {
280 sele2ToIndex_[istep].resize(index + 1);
281 sele2ToIndex_[istep][index] = mol->getGlobalIndex();
282 }
283 if (selectionModeRestart_) {
284
285 GIDtoSele2_[istep][mol->getGlobalIndex()] = index;
286 if (istep == 0) {
287 selection2StartFrame_[istep].push_back(istep);
288 } else {
289 int prevIndex = GIDtoSele2_[istep-1][mol->getGlobalIndex()];
290 if (prevIndex == -1)
291 selection2StartFrame_[istep].push_back(istep);
292 else
293 selection2StartFrame_[istep].push_back(selection2StartFrame_[istep-1][prevIndex]);
294 }
295 }
296 }
297 }
298 }
299
300 if (doObjectProperties_) {
301 // Map of StuntDouble global IDs to selections
302 if (selectionModeRestart_)
303 GIDtoSele1_[istep].resize(info_->getNGlobalIntegrableObjects(), -1);
304
305 for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
306 sd = seleMan1_.nextSelected(isd1)) {
307 index = computeProperty1(istep, sd);
308
309 if (index == sele1ToIndex_[istep].size()) {
310 sele1ToIndex_[istep].push_back(sd->getGlobalIndex());
311 } else {
312 sele1ToIndex_[istep].resize(index + 1);
313 sele1ToIndex_[istep][index] = sd->getGlobalIndex();
314 }
315
316 if (selectionModeRestart_) {
317
318 GIDtoSele1_[istep][sd->getGlobalIndex()] = index;
319 if (istep == 0) {
320 selection1StartFrame_[istep].push_back(istep);
321 } else {
322 int prevIndex = GIDtoSele1_[istep-1][sd->getGlobalIndex()];
323 if (prevIndex == -1)
324 selection1StartFrame_[istep].push_back(istep);
325 else
326 selection1StartFrame_[istep].push_back(selection1StartFrame_[istep-1][prevIndex]);
327 }
328 }
329
330 if (!uniqueSelections_) { index = computeProperty2(istep, sd); }
331 }
332
333 if (uniqueSelections_) {
334 // Map of StuntDouble global IDs to selections
335 if (selectionModeRestart_)
336 GIDtoSele2_[istep].resize(info_->getNGlobalIntegrableObjects(), -1);
337
338 for (sd = seleMan2_.beginSelected(isd2); sd != NULL;
339 sd = seleMan2_.nextSelected(isd2)) {
340 if (autoCorrFunc_) {
341 index = computeProperty1(istep, sd);
342 } else {
343 index = computeProperty2(istep, sd);
344 }
345 if (index == sele2ToIndex_[istep].size()) {
346 sele2ToIndex_[istep].push_back(sd->getGlobalIndex());
347 } else {
348 sele2ToIndex_[istep].resize(index + 1);
349 sele2ToIndex_[istep][index] = sd->getGlobalIndex();
350 }
351 if (selectionModeRestart_) {
352
353 GIDtoSele2_[istep][sd->getGlobalIndex()] = index;
354 if (istep == 0) {
355 selection2StartFrame_[istep].push_back(istep);
356 } else {
357 int prevIndex = GIDtoSele2_[istep-1][sd->getGlobalIndex()];
358 if (prevIndex == -1)
359 selection2StartFrame_[istep].push_back(istep);
360 else
361 selection2StartFrame_[istep].push_back(selection2StartFrame_[istep-1][prevIndex]);
362 }
363 }
364 }
365 }
366 }
367
368 if (doBondProperties_) {
369 // Map of bond global IDs to selections
370 if (selectionModeRestart_)
371 GIDtoSele1_[istep].resize(info_->getNGlobalBonds(), -1);
372
373 for (bond = seleMan1_.beginSelectedBond(ibond1); bond != NULL;
374 bond = seleMan1_.nextSelectedBond(ibond1)) {
375 index = computeProperty1(istep, bond);
376
377 if (index == sele1ToIndex_[istep].size()) {
378 sele1ToIndex_[istep].push_back(bond->getGlobalIndex());
379 } else {
380 sele1ToIndex_[istep].resize(index + 1);
381 sele1ToIndex_[istep][index] = bond->getGlobalIndex();
382 }
383
384 if (selectionModeRestart_) {
385
386 GIDtoSele1_[istep][bond->getGlobalIndex()] = index;
387 if (istep == 0) {
388 selection1StartFrame_[istep].push_back(istep);
389 } else {
390 int prevIndex = GIDtoSele1_[istep-1][bond->getGlobalIndex()];
391 if (prevIndex == -1)
392 selection1StartFrame_[istep].push_back(istep);
393 else
394 selection1StartFrame_[istep].push_back(selection1StartFrame_[istep-1][prevIndex]);
395 }
396 }
397
398 if (!uniqueSelections_) { index = computeProperty2(istep, bond); }
399 }
400
401 if (uniqueSelections_) {
402 // Map of bond global IDs to selections
403 if (selectionModeRestart_)
404 GIDtoSele2_[istep].resize(info_->getNGlobalBonds(), -1);
405
406 for (bond = seleMan2_.beginSelectedBond(ibond2); bond != NULL;
407 bond = seleMan2_.nextSelectedBond(ibond2)) {
408 if (autoCorrFunc_) {
409 index = computeProperty1(istep, bond);
410 } else {
411 index = computeProperty2(istep, bond);
412 }
413 if (index == sele2ToIndex_[istep].size()) {
414 sele2ToIndex_[istep].push_back(bond->getGlobalIndex());
415 } else {
416 sele2ToIndex_[istep].resize(index + 1);
417 sele2ToIndex_[istep][index] = bond->getGlobalIndex();
418 }
419 if (selectionModeRestart_) {
420
421 GIDtoSele2_[istep][bond->getGlobalIndex()] = index;
422 if (istep == 0) {
423 selection2StartFrame_[istep].push_back(istep);
424 } else {
425 int prevIndex = GIDtoSele2_[istep-1][bond->getGlobalIndex()];
426 if (prevIndex == -1)
427 selection2StartFrame_[istep].push_back(istep);
428 else
429 selection2StartFrame_[istep].push_back(selection2StartFrame_[istep-1][prevIndex]);
430 }
431 }
432 }
433 }
434 }
435 }
436
437 // Setters — call before doCorrelate()
438 template<typename T>
439 void TimeCorrFunc<T>::setWindowingParameters(RealType tcorr_fs, int nStart, RealType tsep_fs) {
440 // deltaTime_ and nFrames_ are already set by the constructor
441 // so we can compute frame counts here.
442 nStart_ = nStart;
443
444 // tsep_fs > 0: gap between end of one window and start of the next
445 // (non-overlapping with decorrelation gap)
446 // tsep_fs = 0: tight non-overlapping windows (nStride = nTimeBins)
447 // tsep_fs < 0: overlapping windows.
448 // |tsep_fs| is the overlap in femtoseconds, i.e. the
449 // amount by which successive windows share frames.
450 // nStride = nTimeBins + nSep, where nSep < 0.
451 // nStride is clamped to a minimum of 1 so that window
452 // origins always advance by at least one frame.
453
454 nSep_ = static_cast<int>(tsep_fs / deltaTime_);
455 int nBins = (tcorr_fs > 0)
456 ? static_cast<int>(tcorr_fs / deltaTime_)
457 : nFrames_;
458 nStride_ = std::max(1, nBins + nSep_);
459 useWindowing_ = (nStride_ != nBins || nStart_ > 0);
460 nTimeBins_ = useWindowing_ ? nBins : nFrames_;
461
462 T zeroType(0.0);
463 histogram_.assign(nTimeBins_, zeroType);
464 count_.assign(nTimeBins_, 0);
465 }
466
467 template<typename T>
468 void TimeCorrFunc<T>::doCorrelate() {
469 painCave.isFatal = 0;
470 painCave.severity = OPENMD_INFO;
471 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
472 "Starting pre-correlate scan.");
473 simError();
474 preCorrelate();
475
476 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
477 "Calculating correlation function.");
478 simError();
479 correlation();
480
481 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
482 "Doing post-correlation calculations.");
483 simError();
484 postCorrelate();
485
486 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "Writing output.");
487 simError();
488 writeCorrelate();
489 }
490
491 template<typename T>
492 void TimeCorrFunc<T>::correlation() {
493 T zeroType(0.0);
494 for (unsigned int i = 0; i < nTimeBins_; ++i) {
495 histogram_[i] = zeroType;
496 count_[i] = 0;
497 }
498
499 progressBar_->clear();
500
501 if (useWindowing_) {
502 // -------------------------------------------------------------------
503 // Windowed mode: non-overlapping (or strided) origins.
504 //
505 // For each window ii, the origin frame is:
506 // frame0 = nStart_ + ii * nStride_
507 // The window covers frames [frame0, frame0 + nTimeBins_).
508 // -------------------------------------------------------------------
509 RealType samples = static_cast<RealType>(navg_ * nTimeBins_);
510 int visited = 0;
511
512 for (int ii = 0; ii < navg_; ++ii) {
513 int frame0 = nStart_ + ii * nStride_;
514
515 for (int tt = 0; tt < nTimeBins_; ++tt) {
516 int frame2 = frame0 + tt;
517 if (frame2 >= nFrames_) break;
518
519 visited++;
520 progressBar_->setStatus(visited, samples);
521 progressBar_->update();
522
523 // timeBin is just tt — lag from this window's origin
524 correlateFrames(frame0, frame2, tt);
525 }
526 }
527
528 } else {
529 // -------------------------------------------------------------------
530 // Original overlapping-origin mode (nStride_ == 1).
531 // All frame pairs (i, j) with j >= i are correlated.
532 // -------------------------------------------------------------------
533 RealType samples = 0.5 * (nFrames_ + 1) * nFrames_;
534 int visited = 0;
535
536 for (int i = 0; i < nFrames_; ++i) {
537 RealType time1 = times_[i];
538
539 for (int j = i; j < nFrames_; ++j) {
540 visited++;
541 progressBar_->setStatus(visited, samples);
542 progressBar_->update();
543
544 RealType time2 = times_[j];
545
546 if (std::fabs((time2 - time1) - (j - i) * dtMean_) >
547 6 * dtSigma_ * (j - i)) {
548 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
549 "TimeCorrFunc::correlation: mean sampleTime (%f)\n"
550 "\tin %s does not match actual time-spacing between\n"
551 "\tconfigurations %d (t = %f) and %d (t = %f).\n",
552 dtMean_, dumpFilename_.c_str(), i, time1, j, time2);
553 if (allowTimeFuzz_) {
554 painCave.isFatal = 0;
555 painCave.severity = OPENMD_INFO;
556 } else {
557 painCave.isFatal = 1;
558 painCave.severity = OPENMD_ERROR;
559 }
560 simError();
561 }
562
563 int timeBin = (nFrames_== 1) ? 0 : int((time2 - time1)/dtMean_ + 0.5);
564 correlateFrames(i, j, timeBin);
565 }
566 }
567 }
568 }
569
570
571 /*
572template<typename T>
573void TimeCorrFunc<T>::validateSelection(SelectionManager& seleMan) {
574}
575*/
576
577 template<typename T>
578 void TimeCorrFunc<T>::correlateFrames(int frame1, int frame2, int timeBin) {
579 std::vector<int> s1;
580 std::vector<int> s2;
581
582 std::vector<int>::iterator i1;
583 std::vector<int>::iterator i2;
584
585 T corrVal(0.0);
586 T zeroType(0.0);
587
588 if (doSystemProperties_) {
589 corrVal = calcCorrVal(frame1, frame2);
590 histogram_[timeBin] += corrVal;
591 count_[timeBin]++;
592
593 } else {
594 s1 = sele1ToIndex_[frame1];
595
596 if (uniqueSelections_)
597 s2 = sele2ToIndex_[frame2];
598 else
599 s2 = sele1ToIndex_[frame2];
600
601 for (i1 = s1.begin(), i2 = s2.begin(); i1 != s1.end() && i2 != s2.end();
602 ++i1, ++i2) {
603 // If the selections are dynamic, they might not have the
604 // same objects in both frames, so we need to roll either of
605 // the selections until we have the same object to
606 // correlate.
607
608 while (i1 != s1.end() && *i1 < *i2) {
609 ++i1;
610 }
611
612 while (i2 != s2.end() && *i2 < *i1) {
613 ++i2;
614 }
615
616 if (i1 == s1.end() || i2 == s2.end()) break;
617
618 if (selectionModeRestart_) {
619
620 int ssf1 = selection1StartFrame_[frame1][i1 - s1.begin()];
621 int ssf2 = uniqueSelections_ ?
622 selection2StartFrame_[frame2][i2 - s2.begin()] :
623 selection1StartFrame_[frame2][i2 - s2.begin()];
624
625
626 if (ssf1 != ssf2) break;
627 corrVal = calcCorrVal(frame1, frame2,
628 i1 - s1.begin(), i2 - s2.begin());
629 } else {
630
631 corrVal = calcCorrVal(frame1, frame2,
632 i1 - s1.begin(), i2 - s2.begin());
633 }
634 histogram_[timeBin] += corrVal;
635 count_[timeBin]++;
636 }
637 }
638 }
639
640 template<typename T>
641 void TimeCorrFunc<T>::postCorrelate() {
642 T zeroType(0.0);
643 for (unsigned int i = 0; i < nTimeBins_; ++i) {
644 if (count_[i] > 0) {
645 histogram_[i] /= count_[i];
646 } else {
647 histogram_[i] = zeroType;
648 }
649 }
650 }
651
652 template<typename T>
653 void TimeCorrFunc<T>::validateSelection(SelectionManager&) {}
654
655 template<typename T>
656 void TimeCorrFunc<T>::writeCorrelate() {
657 ofstream ofs(outputFilename_.c_str());
658
659 if (ofs.is_open()) {
660 Revision r;
661
662 ofs << "# " << getCorrFuncType() << "\n";
663 ofs << "# OpenMD " << r.getFullRevision() << "\n";
664 ofs << "# " << r.getBuildDate() << "\n";
665 ofs << "# selection script1: \"" << selectionScript1_;
666 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
667 if (!paramString_.empty())
668 ofs << "# parameters: " << paramString_ << "\n";
669 if (!labelString_.empty())
670 ofs << "#time\t" << labelString_ << "\n";
671 else
672 ofs << "#time\tcorrVal\n";
673
674 for (unsigned int i = 0; i < nTimeBins_; ++i) {
675 if (count_[i] > 0) {
676 RealType lagTime = useWindowing_ ?
677 static_cast<RealType>(i) * dtMean_ :
678 times_[i] - times_[0];
679 ofs << lagTime << "\t" << histogram_[i] << "\n";
680 }
681 }
682 } else {
683 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
684 "TimeCorrFunc::writeCorrelate Error: fail to open %s\n",
685 outputFilename_.c_str());
686 painCave.isFatal = 1;
687 simError();
688 }
689
690 ofs.close();
691 }
692
693 // Template specialization of writeCorrelate for Vector3d
694 template<>
695 void TimeCorrFunc<Vector3d>::writeCorrelate() {
696 ofstream ofs(outputFilename_.c_str());
697
698 if (ofs.is_open()) {
699 Revision r;
700
701 ofs << "# " << getCorrFuncType() << "\n";
702 ofs << "# OpenMD " << r.getFullRevision() << "\n";
703 ofs << "# " << r.getBuildDate() << "\n";
704 ofs << "# selection script1: \"" << selectionScript1_;
705 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
706 if (!paramString_.empty())
707 ofs << "# parameters: " << paramString_ << "\n";
708 if (!labelString_.empty())
709 ofs << "#time\t" << labelString_ << "\n";
710 else
711 ofs << "#time\tcorrVal\n";
712
713 for (unsigned int i = 0; i < nTimeBins_; ++i) {
714 if (count_[i] > 0) {
715
716 ofs << times_[i] - times_[0] << "\t";
717 for (int j = 0; j < 3; j++) {
718 ofs << histogram_[i](j) << '\t';
719 }
720 ofs << '\n';
721
722 }
723 }
724
725 } else {
726 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
727 "TimeCorrFunc::writeCorrelate Error: fail to open %s\n",
728 outputFilename_.c_str());
729 painCave.isFatal = 1;
730 simError();
731 }
732
733 ofs.close();
734 }
735
736 // Template specialization of writeCorrelate for Mat3x3d
737 template<>
738 void TimeCorrFunc<Mat3x3d>::writeCorrelate() {
739 ofstream ofs(outputFilename_.c_str());
740
741 if (ofs.is_open()) {
742 Revision r;
743
744 ofs << "# " << getCorrFuncType() << "\n";
745 ofs << "# OpenMD " << r.getFullRevision() << "\n";
746 ofs << "# " << r.getBuildDate() << "\n";
747 ofs << "# selection script1: \"" << selectionScript1_;
748 ofs << "\"\tselection script2: \"" << selectionScript2_ << "\"\n";
749 if (!paramString_.empty())
750 ofs << "# parameters: " << paramString_ << "\n";
751 if (!labelString_.empty())
752 ofs << "#time\t" << labelString_ << "\n";
753 else
754 ofs << "#time\tcorrVal\n";
755
756 for (unsigned int i = 0; i < nTimeBins_; ++i) {
757
758 if (count_[i] > 0) {
759
760 ofs << times_[i] - times_[0] << "\t";
761 for (int j = 0; j < 3; j++) {
762 for (int k = 0; k < 3; k++) {
763 ofs << histogram_[i](j, k) << '\t';
764 }
765 }
766 ofs << '\n';
767 }
768 }
769
770 } else {
771 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
772 "TimeCorrFunc::writeCorrelate Error: fail to open %s\n",
773 outputFilename_.c_str());
774 painCave.isFatal = 1;
775 simError();
776 }
777
778 ofs.close();
779 }
780
781 // it is necessary to keep the constructor definitions here or the code wont
782 // be generated and linking issues will occur. Blame templating
783 template<typename T>
784 CrossCorrFunc<T>::CrossCorrFunc(SimInfo* info, const std::string& filename,
785 const std::string& sele1,
786 const std::string& sele2) :
787 TimeCorrFunc<T>(info, filename, sele1, sele2) {
788 this->autoCorrFunc_ = false;
789 }
790
791 template<typename T>
792 AutoCorrFunc<T>::AutoCorrFunc(SimInfo* info, const std::string& filename,
793 const std::string& sele1,
794 const std::string& sele2) :
795 TimeCorrFunc<T>(info, filename, sele1, sele2) {
796 this->autoCorrFunc_ = true;
797 }
798
799 template<typename T>
800 SystemACF<T>::SystemACF(SimInfo* info, const std::string& filename,
801 const std::string& sele1, const std::string& sele2) :
802 AutoCorrFunc<T>(info, filename, sele1, sele2) {
803 this->autoCorrFunc_ = true;
804 this->doSystemProperties_ = true;
805 this->doMolecularProperties_ = false;
806 this->doObjectProperties_ = false;
807 this->doBondProperties_ = false;
808 }
809
810 template<typename T>
811 SystemCCF<T>::SystemCCF(SimInfo* info, const std::string& filename,
812 const std::string& sele1, const std::string& sele2) :
813 CrossCorrFunc<T>(info, filename, sele1, sele2) {
814 this->autoCorrFunc_ = false;
815 this->doSystemProperties_ = true;
816 this->doMolecularProperties_ = false;
817 this->doObjectProperties_ = false;
818 this->doBondProperties_ = false;
819 }
820
821 template<typename T>
822 ObjectACF<T>::ObjectACF(SimInfo* info, const std::string& filename,
823 const std::string& sele1, const std::string& sele2) :
824 AutoCorrFunc<T>(info, filename, sele1, sele2) {
825 this->autoCorrFunc_ = true;
826 this->doSystemProperties_ = false;
827 this->doMolecularProperties_ = false;
828 this->doObjectProperties_ = true;
829 this->doBondProperties_ = false;
830 }
831
832 template<typename T>
833 ObjectCCF<T>::ObjectCCF(SimInfo* info, const std::string& filename,
834 const std::string& sele1, const std::string& sele2) :
835 CrossCorrFunc<T>(info, filename, sele1, sele2) {
836 this->autoCorrFunc_ = false;
837 this->doSystemProperties_ = false;
838 this->doMolecularProperties_ = false;
839 this->doObjectProperties_ = true;
840 this->doBondProperties_ = false;
841 }
842
843 template<typename T>
844 MoleculeACF<T>::MoleculeACF(SimInfo* info, const std::string& filename,
845 const std::string& sele1,
846 const std::string& sele2) :
847 AutoCorrFunc<T>(info, filename, sele1, sele2) {
848 this->autoCorrFunc_ = true;
849 this->doSystemProperties_ = false;
850 this->doMolecularProperties_ = true;
851 this->doObjectProperties_ = false;
852 this->doBondProperties_ = false;
853 }
854
855 template<typename T>
856 MoleculeCCF<T>::MoleculeCCF(SimInfo* info, const std::string& filename,
857 const std::string& sele1,
858 const std::string& sele2) :
859 CrossCorrFunc<T>(info, filename, sele1, sele2) {
860 this->autoCorrFunc_ = false;
861 this->doSystemProperties_ = false;
862 this->doMolecularProperties_ = true;
863 this->doObjectProperties_ = false;
864 this->doBondProperties_ = false;
865 }
866
867 template class AutoCorrFunc<RealType>;
868 template class TimeCorrFunc<RealType>;
869 template class CrossCorrFunc<RealType>;
870
871 template class AutoCorrFunc<Vector3d>;
872 template class TimeCorrFunc<Vector3d>;
873 template class CrossCorrFunc<Vector3d>;
874
875 template class AutoCorrFunc<Mat3x3d>;
876 template class TimeCorrFunc<Mat3x3d>;
877 template class CrossCorrFunc<Mat3x3d>;
878
879 template class TimeCorrFunc<DynamicVector<RealType>>;
880
881 template class AutoCorrFunc<Vector<RealType, 4>>;
882 template class TimeCorrFunc<Vector<RealType, 4>>;
883 template class CrossCorrFunc<Vector<RealType, 4>>;
884
885 template class SystemACF<RealType>;
886 template class SystemACF<Vector3d>;
887 template class SystemACF<Mat3x3d>;
888
889 template class SystemCCF<RealType>;
890 template class SystemCCF<Vector3d>;
891 template class SystemCCF<Mat3x3d>;
892
893 template class ObjectACF<RealType>;
894 template class ObjectACF<Vector3d>;
895 template class ObjectACF<Mat3x3d>;
896
897 template class ObjectCCF<RealType>;
898 template class ObjectCCF<Vector3d>;
899 template class ObjectCCF<Mat3x3d>;
900 template class ObjectCCF<int>;
901
902 template class MoleculeACF<RealType>;
903 template class MoleculeACF<Vector3d>;
904 template class MoleculeACF<Mat3x3d>;
905 template class MoleculeACF<Vector<RealType, 4>>;
906
907 template class MoleculeCCF<RealType>;
908 template class MoleculeCCF<Vector3d>;
909 template class MoleculeCCF<Mat3x3d>;
910} // namespace OpenMD
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:96
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:63