OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
AccumulatorView.hpp
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#ifndef OPENMD_UTILS_ACCUMULATORVIEW_HPP
49#define OPENMD_UTILS_ACCUMULATORVIEW_HPP
50
51#include <cmath>
52#include <cstddef>
53#include <iostream>
54#include <string>
55#include <typeinfo>
56#include <vector>
57
58#include "math/SquareMatrix3.hpp"
59#include "math/Vector3.hpp"
60#include "nonbonded/NonBondedInteraction.hpp"
61#include "utils/Accumulator.hpp"
62#include "utils/BaseAccumulator.hpp"
63#include "utils/simError.h"
64
65namespace OpenMD::Utils {
66
67 template<typename T>
68 class AccumulatorView;
69
70 template<>
71 class AccumulatorView<RealAccumulator> :
72 public BaseAccumulator,
73 private RealAccumulator {
74 public:
75 void add(RealType val) override { RealAccumulator::add(val); }
76
77 // Other overrides for invalid entries
78 void add(const std::vector<RealType>&) override {
79 accumulatorFunctionCallMismatch();
80 }
81 void add(const Vector3d&) override { accumulatorFunctionCallMismatch(); }
82 void add(const potVec&) override { accumulatorFunctionCallMismatch(); }
83 void add(const Mat3x3d&) override { accumulatorFunctionCallMismatch(); }
84
85 void writeData(std::ostream& stream, const std::string& errorMessage,
86 DataHandling dataHandling) const override {
87 RealType dat {};
88 std::size_t count = RealAccumulator::getCount();
89
90 switch (dataHandling) {
91 case DataHandling::Average:
92 dat = RealAccumulator::getAverage();
93 break;
94 case DataHandling::Last:
95 dat = RealAccumulator::getLastValue();
96 break;
97 case DataHandling::Max:
98 dat = RealAccumulator::getMax();
99 break;
100 case DataHandling::Min:
101 dat = RealAccumulator::getMin();
102 break;
103 case DataHandling::Total:
104 dat = RealAccumulator::getTotal();
105 break;
106 default:
107 break;
108 }
109
110 if (std::isinf(dat) || std::isnan(dat)) {
111 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
112 errorMessage.c_str());
113 painCave.isFatal = 1;
114 simError();
115 } else {
116 if (count == 0)
117 stream << "\t";
118 else
119 stream << "\t" << dat;
120 }
121 }
122
123 void writeErrorBars(std::ostream& stream, const std::string& errorMessage,
124 ErrorHandling errorHandling) const override {
125 RealType err {};
126 std::size_t count = RealAccumulator::getCount();
127
128 switch (errorHandling) {
129 case ErrorHandling::CI95:
130 err = RealAccumulator::get95percentConfidenceInterval();
131 break;
132 case ErrorHandling::StdDev:
133 err = RealAccumulator::getStdDev();
134 break;
135 case ErrorHandling::Variance:
136 err = RealAccumulator::getVariance();
137 break;
138 default:
139 break;
140 }
141
142 if (count == 0 && errorHandling == ErrorHandling::CI95) {
143 stream << "\t";
144 } else {
145 if (std::isinf(err) || std::isnan(err)) {
146 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
147 errorMessage.c_str());
148 painCave.isFatal = 1;
149 simError();
150 } else {
151 if (count == 0)
152 stream << "\t";
153 else
154 stream << "\t" << err;
155 }
156 }
157 }
158
159 std::type_index getType() const override { return typeid(RealType); }
160
161 std::size_t getCount() const override {
162 return RealAccumulator::getCount();
163 }
164 };
165
166 template<>
167 class AccumulatorView<StdVectorAccumulator> :
168 public BaseAccumulator,
169 private StdVectorAccumulator {
170 public:
171 void add(const std::vector<RealType>& val) override {
172 StdVectorAccumulator::add(val);
173 }
174
175 // Other overrides for invalid entries
176 void add(RealType val) override { accumulatorFunctionCallMismatch(); }
177 void add(const Vector3d&) override { accumulatorFunctionCallMismatch(); }
178 void add(const potVec&) override { accumulatorFunctionCallMismatch(); }
179 void add(const Mat3x3d&) override { accumulatorFunctionCallMismatch(); }
180
181 void writeData(std::ostream& stream, const std::string& errorMessage,
182 DataHandling dataHandling) const override {
183 std::vector<RealType> dat;
184 std::size_t count = StdVectorAccumulator::getCount();
185
186 switch (dataHandling) {
187 case DataHandling::Average:
188 dat = StdVectorAccumulator::getAverage();
189 break;
190 case DataHandling::Last:
191 dat = StdVectorAccumulator::getLastValue();
192 break;
193 case DataHandling::Total:
194 dat = StdVectorAccumulator::getTotal();
195 break;
196 case DataHandling::Max:
197 case DataHandling::Min:
198 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
199 "Max and Min functions are not defined for a "
200 "std::vector Accumulator.");
201 painCave.isFatal = 1;
202 simError();
203 break;
204 default:
205 break;
206 }
207
208 for (int i = 0; i < StdVectorAccumulator::getAverage().size(); i++) {
209 if (std::isinf(dat[i]) || std::isnan(dat[i])) {
210 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
211 errorMessage.c_str());
212 painCave.isFatal = 1;
213 simError();
214 } else {
215 if (count == 0)
216 stream << '\t';
217 else
218 stream << '\t' << dat[i];
219 }
220 }
221 }
222
223 void writeErrorBars(std::ostream& stream, const std::string& errorMessage,
224 ErrorHandling errorHandling) const override {
225 std::vector<RealType> err;
226 std::size_t count = StdVectorAccumulator::getCount();
227
228 switch (errorHandling) {
229 case ErrorHandling::CI95:
230 err = StdVectorAccumulator::get95percentConfidenceInterval();
231 break;
232 case ErrorHandling::StdDev:
233 err = StdVectorAccumulator::getStdDev();
234 break;
235 case ErrorHandling::Variance:
236 err = StdVectorAccumulator::getVariance();
237 break;
238 default:
239 break;
240 }
241
242 for (int i = 0; i < StdVectorAccumulator::getAverage().size(); i++) {
243 if (count == 0 && errorHandling == ErrorHandling::CI95) {
244 stream << "\t";
245 } else {
246 if (std::isinf(err[i]) || std::isnan(err[i])) {
247 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
248 errorMessage.c_str());
249 painCave.isFatal = 1;
250 simError();
251 } else {
252 if (count == 0)
253 stream << "\t";
254 else
255 stream << "\t" << err[i];
256 }
257 }
258 }
259 }
260
261 std::type_index getType() const override {
262 return typeid(std::vector<RealType>);
263 }
264
265 std::size_t getCount() const override {
266 std::size_t count = StdVectorAccumulator::getCount();
267
268 return count;
269 }
270 };
271
272 template<>
273 class AccumulatorView<Vector3dAccumulator> :
274 public BaseAccumulator,
275 private Vector3dAccumulator {
276 public:
277 void add(const Vector3d& val) override { Vector3dAccumulator::add(val); }
278
279 // Other overrides for invalid entries
280 void add(RealType val) override { accumulatorFunctionCallMismatch(); }
281 void add(const std::vector<RealType>&) override {
282 accumulatorFunctionCallMismatch();
283 }
284 void add(const potVec&) override { accumulatorFunctionCallMismatch(); }
285 void add(const Mat3x3d&) override { accumulatorFunctionCallMismatch(); }
286
287 void writeData(std::ostream& stream, const std::string& errorMessage,
288 DataHandling dataHandling) const override {
289 Vector3d dat;
290 std::size_t count = Vector3dAccumulator::getCount();
291
292 switch (dataHandling) {
293 case DataHandling::Average:
294 dat = Vector3dAccumulator::getAverage();
295 break;
296 case DataHandling::Last:
297 dat = Vector3dAccumulator::getLastValue();
298 break;
299 case DataHandling::Total:
300 dat = Vector3dAccumulator::getTotal();
301 break;
302 case DataHandling::Max:
303 case DataHandling::Min:
304 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
305 "Max and Min functions are not defined for a "
306 "std::vector Accumulator.");
307 painCave.isFatal = 1;
308 simError();
309 break;
310 default:
311 break;
312 }
313
314 for (int i = 0; i < Vector3dAccumulator::getAverage().size(); i++) {
315 if (std::isinf(dat[i]) || std::isnan(dat[i])) {
316 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
317 errorMessage.c_str());
318 painCave.isFatal = 1;
319 simError();
320 } else {
321 if (count == 0)
322 stream << "\t";
323 else
324 stream << "\t" << dat[i];
325 }
326 }
327 }
328
329 void writeErrorBars(std::ostream& stream, const std::string& errorMessage,
330 ErrorHandling errorHandling) const override {
331 Vector3d err;
332 std::size_t count = Vector3dAccumulator::getCount();
333
334 switch (errorHandling) {
335 case ErrorHandling::CI95:
336 err = Vector3dAccumulator::get95percentConfidenceInterval();
337 break;
338 case ErrorHandling::StdDev:
339 err = Vector3dAccumulator::getStdDev();
340 break;
341 case ErrorHandling::Variance:
342 err = Vector3dAccumulator::getVariance();
343 break;
344 default:
345 break;
346 }
347
348 for (int i = 0; i < Vector3dAccumulator::getAverage().size(); i++) {
349 if (count == 0 && errorHandling == ErrorHandling::CI95) {
350 stream << "\t";
351 } else {
352 if (std::isinf(err[i]) || std::isnan(err[i])) {
353 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
354 errorMessage.c_str());
355 painCave.isFatal = 1;
356 simError();
357 } else {
358 if (count == 0)
359 stream << "\t";
360 else
361 stream << "\t" << err[i];
362 }
363 }
364 }
365 }
366
367 std::type_index getType() const override { return typeid(Vector3d); }
368
369 std::size_t getCount() const override {
370 return Vector3dAccumulator::getCount();
371 }
372 };
373
374 template<>
375 class AccumulatorView<PotVecAccumulator> :
376 public BaseAccumulator,
377 private PotVecAccumulator {
378 public:
379 void add(const potVec& val) override { PotVecAccumulator::add(val); }
380
381 // Other overrides for invalid entries
382 void add(RealType val) override { accumulatorFunctionCallMismatch(); }
383 void add(const std::vector<RealType>&) override {
384 accumulatorFunctionCallMismatch();
385 }
386 void add(const Vector3d&) override { accumulatorFunctionCallMismatch(); }
387 void add(const Mat3x3d&) override { accumulatorFunctionCallMismatch(); }
388
389 void writeData(std::ostream& stream, const std::string& errorMessage,
390 DataHandling dataHandling) const override {
391 potVec dat;
392 std::size_t count = PotVecAccumulator::getCount();
393
394 switch (dataHandling) {
395 case DataHandling::Average:
396 dat = PotVecAccumulator::getAverage();
397 break;
398 case DataHandling::Last:
399 dat = PotVecAccumulator::getLastValue();
400 break;
401 case DataHandling::Total:
402 dat = PotVecAccumulator::getTotal();
403 break;
404 case DataHandling::Max:
405 case DataHandling::Min:
406 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
407 "Max and Min functions are not defined for a "
408 "std::vector Accumulator.");
409 painCave.isFatal = 1;
410 simError();
411 break;
412 default:
413 break;
414 }
415
416 for (int i = 0; i < PotVecAccumulator::getAverage().size(); i++) {
417 if (std::isinf(dat[i]) || std::isnan(dat[i])) {
418 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
419 errorMessage.c_str());
420 painCave.isFatal = 1;
421 simError();
422 } else {
423 if (count == 0)
424 stream << "\t";
425 else
426 stream << "\t" << dat[i];
427 }
428 }
429 }
430
431 void writeErrorBars(std::ostream& stream, const std::string& errorMessage,
432 ErrorHandling errorHandling) const override {
433 potVec err;
434 std::size_t count = PotVecAccumulator::getCount();
435
436 switch (errorHandling) {
437 case ErrorHandling::CI95:
438 err = PotVecAccumulator::get95percentConfidenceInterval();
439 break;
440 case ErrorHandling::StdDev:
441 err = PotVecAccumulator::getStdDev();
442 break;
443 case ErrorHandling::Variance:
444 err = PotVecAccumulator::getVariance();
445 break;
446 default:
447 break;
448 }
449
450 for (int i = 0; i < PotVecAccumulator::getAverage().size(); i++) {
451 if (count == 0 && errorHandling == ErrorHandling::CI95) {
452 stream << "\t";
453 } else {
454 if (std::isinf(err[i]) || std::isnan(err[i])) {
455 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
456 errorMessage.c_str());
457 painCave.isFatal = 1;
458 simError();
459 } else {
460 if (count == 0)
461 stream << "\t";
462 else
463 stream << "\t" << err[i];
464 }
465 }
466 }
467 }
468
469 std::type_index getType() const override { return typeid(potVec); }
470
471 std::size_t getCount() const override {
472 return PotVecAccumulator::getCount();
473 }
474 };
475
476 template<>
477 class AccumulatorView<Mat3x3dAccumulator> :
478 public BaseAccumulator,
479 private Mat3x3dAccumulator {
480 public:
481 void add(const Mat3x3d& val) override { Mat3x3dAccumulator::add(val); }
482
483 // Other overrides for invalid entries
484 void add(RealType val) override { accumulatorFunctionCallMismatch(); }
485 void add(const std::vector<RealType>&) override {
486 accumulatorFunctionCallMismatch();
487 }
488 void add(const Vector3d&) override { accumulatorFunctionCallMismatch(); }
489 void add(const potVec&) override { accumulatorFunctionCallMismatch(); }
490
491 void writeData(std::ostream& stream, const std::string& errorMessage,
492 DataHandling dataHandling) const override {
493 Mat3x3d dat;
494 std::size_t count = Mat3x3dAccumulator::getCount();
495
496 switch (dataHandling) {
497 case DataHandling::Average:
498 dat = Mat3x3dAccumulator::getAverage();
499 break;
500 case DataHandling::Last:
501 dat = Mat3x3dAccumulator::getLastValue();
502 break;
503 case DataHandling::Total:
504 dat = Mat3x3dAccumulator::getTotal();
505 break;
506 case DataHandling::Max:
507 case DataHandling::Min:
508 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
509 "Max and Min functions are not defined for a "
510 "std::vector Accumulator.");
511 painCave.isFatal = 1;
512 simError();
513 break;
514 default:
515 break;
516 }
517
518 for (unsigned int i = 0; i < 3; i++) {
519 for (unsigned int j = 0; j < 3; j++) {
520 if (std::isinf(dat(i, j)) || std::isnan(dat(i, j))) {
521 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
522 errorMessage.c_str());
523 painCave.isFatal = 1;
524 simError();
525 } else {
526 if (count == 0)
527 stream << "\t";
528 else
529 stream << "\t" << dat(i, j);
530 }
531 }
532 }
533 }
534
535 void writeErrorBars(std::ostream& stream, const std::string& errorMessage,
536 ErrorHandling errorHandling) const override {
537 Mat3x3d err;
538 std::size_t count = Mat3x3dAccumulator::getCount();
539
540 switch (errorHandling) {
541 case ErrorHandling::CI95:
542 err = Mat3x3dAccumulator::get95percentConfidenceInterval();
543 break;
544 case ErrorHandling::StdDev:
545 err = Mat3x3dAccumulator::getStdDev();
546 break;
547 case ErrorHandling::Variance:
548 err = Mat3x3dAccumulator::getVariance();
549 break;
550 default:
551 break;
552 }
553
554 for (unsigned int i = 0; i < 3; i++) {
555 for (unsigned int j = 0; j < 3; j++) {
556 if (count == 0 && errorHandling == ErrorHandling::CI95) {
557 stream << "\t";
558 } else {
559 if (std::isinf(err(i, j)) || std::isnan(err(i, j))) {
560 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH, "%s",
561 errorMessage.c_str());
562 painCave.isFatal = 1;
563 simError();
564 } else {
565 if (count == 0)
566 stream << "\t";
567 else
568 stream << "\t" << err(i, j);
569 }
570 }
571 }
572 }
573 }
574
575 std::type_index getType() const override { return typeid(Mat3x3d); }
576
577 std::size_t getCount() const override {
578 return Mat3x3dAccumulator::getCount();
579 }
580 };
581} // namespace OpenMD::Utils
582
583#endif