OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Stats.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 "brains/Stats.hpp"
49
50#include <iomanip>
51#include <sstream>
52
53#include "brains/Thermo.hpp"
54#include "utils/MemoryUtils.hpp"
55
56namespace OpenMD {
57
58 Stats::Stats(SimInfo* info) : info_(info), isInit_(false) {
59 if (!isInit_) {
60 init();
61 isInit_ = true;
62 }
63 }
64
65 Stats::~Stats() {
66 for (auto& data : data_) {
67 if (!data.accumulatorArray2d.empty())
68 Utils::deletePointers(data.accumulatorArray2d);
69 else
70 delete data.accumulator;
71 }
72 }
73
74 void Stats::init() {
75 data_.resize(Stats::ENDINDEX);
76
77 StatsData time;
78 time.units = "fs";
79 time.title = "Time";
80 time.dataType = "RealType";
81 time.accumulator = new Accumulator();
82 data_[TIME] = time;
83 statsMap_["TIME"] = TIME;
84
85 StatsData total_energy;
86 total_energy.units = "kcal/mol";
87 total_energy.title = "Total Energy";
88 total_energy.dataType = "RealType";
89 total_energy.accumulator = new Accumulator();
90 data_[TOTAL_ENERGY] = total_energy;
91 statsMap_["TOTAL_ENERGY"] = TOTAL_ENERGY;
92
93 StatsData potential_energy;
94 potential_energy.units = "kcal/mol";
95 potential_energy.title = "Potential Energy";
96 potential_energy.dataType = "RealType";
97 potential_energy.accumulator = new Accumulator();
98 data_[POTENTIAL_ENERGY] = potential_energy;
99 statsMap_["POTENTIAL_ENERGY"] = POTENTIAL_ENERGY;
100
101 StatsData kinetic_energy;
102 kinetic_energy.units = "kcal/mol";
103 kinetic_energy.title = "Kinetic Energy";
104 kinetic_energy.dataType = "RealType";
105 kinetic_energy.accumulator = new Accumulator();
106 data_[KINETIC_ENERGY] = kinetic_energy;
107 statsMap_["KINETIC_ENERGY"] = KINETIC_ENERGY;
108
109 StatsData temperature;
110 temperature.units = "K";
111 temperature.title = "Temperature";
112 temperature.dataType = "RealType";
113 temperature.accumulator = new Accumulator();
114 data_[TEMPERATURE] = temperature;
115 statsMap_["TEMPERATURE"] = TEMPERATURE;
116
117 StatsData pressure;
118 pressure.units = "atm";
119 pressure.title = "Pressure";
120 pressure.dataType = "RealType";
121 pressure.accumulator = new Accumulator();
122 data_[PRESSURE] = pressure;
123 statsMap_["PRESSURE"] = PRESSURE;
124
125 StatsData volume;
126 volume.units = "A^3";
127 volume.title = "Volume";
128 volume.dataType = "RealType";
129 volume.accumulator = new Accumulator();
130 data_[VOLUME] = volume;
131 statsMap_["VOLUME"] = VOLUME;
132
133 StatsData hullvolume;
134 hullvolume.units = "A^3";
135 hullvolume.title = "Hull Volume";
136 hullvolume.dataType = "RealType";
137 hullvolume.accumulator = new Accumulator();
138 data_[HULLVOLUME] = hullvolume;
139 statsMap_["HULLVOLUME"] = HULLVOLUME;
140
141 StatsData gyrvolume;
142 gyrvolume.units = "A^3";
143 gyrvolume.title = "Gyrational Volume";
144 gyrvolume.dataType = "RealType";
145 gyrvolume.accumulator = new Accumulator();
146 data_[GYRVOLUME] = gyrvolume;
147 statsMap_["GYRVOLUME"] = GYRVOLUME;
148
149 StatsData conserved_quantity;
150 conserved_quantity.units = "kcal/mol";
151 conserved_quantity.title = "Conserved Quantity";
152 conserved_quantity.dataType = "RealType";
153 conserved_quantity.accumulator = new Accumulator();
154 data_[CONSERVED_QUANTITY] = conserved_quantity;
155 statsMap_["CONSERVED_QUANTITY"] = CONSERVED_QUANTITY;
156
157 StatsData translational_kinetic;
158 translational_kinetic.units = "kcal/mol";
159 translational_kinetic.title = "Translational Kinetic";
160 translational_kinetic.dataType = "RealType";
161 translational_kinetic.accumulator = new Accumulator();
162 data_[TRANSLATIONAL_KINETIC] = translational_kinetic;
163 statsMap_["TRANSLATIONAL_KINETIC"] = TRANSLATIONAL_KINETIC;
164
165 StatsData rotational_kinetic;
166 rotational_kinetic.units = "kcal/mol";
167 rotational_kinetic.title = "Rotational Kinetic";
168 rotational_kinetic.dataType = "RealType";
169 rotational_kinetic.accumulator = new Accumulator();
170 data_[ROTATIONAL_KINETIC] = rotational_kinetic;
171 statsMap_["ROTATIONAL_KINETIC"] = ROTATIONAL_KINETIC;
172
173 StatsData electronic_kinetic;
174 electronic_kinetic.units = "kcal/mol";
175 electronic_kinetic.title = "Electronic Kinetic";
176 electronic_kinetic.dataType = "RealType";
177 electronic_kinetic.accumulator = new Accumulator();
178 data_[ELECTRONIC_KINETIC] = electronic_kinetic;
179 statsMap_["ELECTRONIC_KINETIC"] = ELECTRONIC_KINETIC;
180
181 StatsData long_range_potential;
182 long_range_potential.units = "kcal/mol";
183 long_range_potential.title = "Long Range Potential";
184 long_range_potential.dataType = "RealType";
185 long_range_potential.accumulator = new Accumulator();
186 data_[LONG_RANGE_POTENTIAL] = long_range_potential;
187 statsMap_["LONG_RANGE_POTENTIAL"] = LONG_RANGE_POTENTIAL;
188
189 StatsData vanderwaals_potential;
190 vanderwaals_potential.units = "kcal/mol";
191 vanderwaals_potential.title = "van der waals Potential";
192 vanderwaals_potential.dataType = "RealType";
193 vanderwaals_potential.accumulator = new Accumulator();
194 data_[VANDERWAALS_POTENTIAL] = vanderwaals_potential;
195 statsMap_["VANDERWAALS_POTENTIAL"] = VANDERWAALS_POTENTIAL;
196
197 StatsData electrostatic_potential;
198 electrostatic_potential.units = "kcal/mol";
199 electrostatic_potential.title = "Electrostatic Potential";
200 electrostatic_potential.dataType = "RealType";
201 electrostatic_potential.accumulator = new Accumulator();
202 data_[ELECTROSTATIC_POTENTIAL] = electrostatic_potential;
203 statsMap_["ELECTROSTATIC_POTENTIAL"] = ELECTROSTATIC_POTENTIAL;
204
205 StatsData metallic_potential;
206 metallic_potential.units = "kcal/mol";
207 metallic_potential.title = "Metallic Potential";
208 metallic_potential.dataType = "RealType";
209 metallic_potential.accumulator = new Accumulator();
210 data_[METALLIC_POTENTIAL] = metallic_potential;
211 statsMap_["METALLIC_POTENTIAL"] = METALLIC_POTENTIAL;
212
213 StatsData metallic_embedding;
214 metallic_embedding.units = "kcal/mol";
215 metallic_embedding.title = "Metallic Embedding";
216 metallic_embedding.dataType = "RealType";
217 metallic_embedding.accumulator = new Accumulator();
218 data_[METALLIC_EMBEDDING] = metallic_embedding;
219 statsMap_["METALLIC_EMBEDDING"] = METALLIC_EMBEDDING;
220
221 StatsData metallic_pair;
222 metallic_pair.units = "kcal/mol";
223 metallic_pair.title = "Metallic Pair";
224 metallic_pair.dataType = "RealType";
225 metallic_pair.accumulator = new Accumulator();
226 data_[METALLIC_PAIR] = metallic_pair;
227 statsMap_["METALLIC_PAIR"] = METALLIC_PAIR;
228
229 StatsData hydrogenbonding_potential;
230 hydrogenbonding_potential.units = "kcal/mol";
231 hydrogenbonding_potential.title = "Hydrogen Bonding Pot.";
232 hydrogenbonding_potential.dataType = "RealType";
233 hydrogenbonding_potential.accumulator = new Accumulator();
234 data_[HYDROGENBONDING_POTENTIAL] = hydrogenbonding_potential;
235 statsMap_["HYDROGENBONDING_POTENTIAL"] = HYDROGENBONDING_POTENTIAL;
236
237 StatsData reciprocal_potential;
238 reciprocal_potential.units = "kcal/mol";
239 reciprocal_potential.title = "Reciprocal Space Pot.";
240 reciprocal_potential.dataType = "RealType";
241 reciprocal_potential.accumulator = new Accumulator();
242 data_[RECIPROCAL_POTENTIAL] = reciprocal_potential;
243 statsMap_["RECIPROCAL_POTENTIAL"] = RECIPROCAL_POTENTIAL;
244
245 StatsData surface_potential;
246 surface_potential.units = "kcal/mol";
247 surface_potential.title = "Surface Potential";
248 surface_potential.dataType = "RealType";
249 surface_potential.accumulator = new Accumulator();
250 data_[SURFACE_POTENTIAL] = surface_potential;
251 statsMap_["SURFACE_POTENTIAL"] = SURFACE_POTENTIAL;
252
253 StatsData short_range_potential;
254 short_range_potential.units = "kcal/mol";
255 short_range_potential.title = "Short Range Potential";
256 short_range_potential.dataType = "RealType";
257 short_range_potential.accumulator = new Accumulator();
258 data_[SHORT_RANGE_POTENTIAL] = short_range_potential;
259 statsMap_["SHORT_RANGE_POTENTIAL"] = SHORT_RANGE_POTENTIAL;
260
261 StatsData bond_potential;
262 bond_potential.units = "kcal/mol";
263 bond_potential.title = "Bond Potential";
264 bond_potential.dataType = "RealType";
265 bond_potential.accumulator = new Accumulator();
266 data_[BOND_POTENTIAL] = bond_potential;
267 statsMap_["BOND_POTENTIAL"] = BOND_POTENTIAL;
268
269 StatsData bend_potential;
270 bend_potential.units = "kcal/mol";
271 bend_potential.title = "Bend Potential";
272 bend_potential.dataType = "RealType";
273 bend_potential.accumulator = new Accumulator();
274 data_[BEND_POTENTIAL] = bend_potential;
275 statsMap_["BEND_POTENTIAL"] = BEND_POTENTIAL;
276
277 StatsData dihedral_potential;
278 dihedral_potential.units = "kcal/mol";
279 dihedral_potential.title = "Dihedral Potential";
280 dihedral_potential.dataType = "RealType";
281 dihedral_potential.accumulator = new Accumulator();
282 data_[DIHEDRAL_POTENTIAL] = dihedral_potential;
283 statsMap_["DIHEDRAL_POTENTIAL"] = DIHEDRAL_POTENTIAL;
284
285 StatsData inversion_potential;
286 inversion_potential.units = "kcal/mol";
287 inversion_potential.title = "Inversion Potential";
288 inversion_potential.dataType = "RealType";
289 inversion_potential.accumulator = new Accumulator();
290 data_[INVERSION_POTENTIAL] = inversion_potential;
291 statsMap_["INVERSION_POTENTIAL"] = INVERSION_POTENTIAL;
292
293 StatsData vraw;
294 vraw.units = "kcal/mol";
295 vraw.title = "Raw Potential";
296 vraw.dataType = "RealType";
297 vraw.accumulator = new Accumulator();
298 data_[RAW_POTENTIAL] = vraw;
299 statsMap_["RAW_POTENTIAL"] = RAW_POTENTIAL;
300
301 StatsData vrestraint;
302 vrestraint.units = "kcal/mol";
303 vrestraint.title = "Restraint Potential";
304 vrestraint.dataType = "RealType";
305 vrestraint.accumulator = new Accumulator();
306 data_[RESTRAINT_POTENTIAL] = vrestraint;
307 statsMap_["RESTRAINT_POTENTIAL"] = RESTRAINT_POTENTIAL;
308
309 StatsData vexcluded;
310 vexcluded.units = "kcal/mol";
311 vexcluded.title = "Excluded Potential";
312 vexcluded.dataType = "RealType";
313 vexcluded.accumulator = new Accumulator();
314 data_[EXCLUDED_POTENTIAL] = vexcluded;
315 statsMap_["EXCLUDED_POTENTIAL"] = EXCLUDED_POTENTIAL;
316
317 StatsData pressure_tensor;
318 pressure_tensor.units = "amu/fs^2/A";
319 pressure_tensor.title = "Pressure Tensor";
320 pressure_tensor.dataType = "Mat3x3d";
321 pressure_tensor.accumulator = new MatrixAccumulator();
322 data_[PRESSURE_TENSOR] = pressure_tensor;
323 statsMap_["PRESSURE_TENSOR"] = PRESSURE_TENSOR;
324
325 // virial tensor added
326 StatsData virial_tensor;
327 virial_tensor.units = "kcal/mol";
328 virial_tensor.title = "Virial Tensor";
329 virial_tensor.dataType = "Mat3x3d";
330 virial_tensor.accumulator = new MatrixAccumulator();
331 data_[VIRIAL_TENSOR] = virial_tensor;
332 statsMap_["VIRIAL_TENSOR"] = VIRIAL_TENSOR;
333
334 StatsData system_dipole;
335 system_dipole.units = "C m";
336 system_dipole.title = "System Dipole";
337 system_dipole.dataType = "Vector3d";
338 system_dipole.accumulator = new VectorAccumulator();
339 data_[SYSTEM_DIPOLE] = system_dipole;
340 statsMap_["SYSTEM_DIPOLE"] = SYSTEM_DIPOLE;
341
342 StatsData system_quadrupole;
343 system_quadrupole.units = "C m^2";
344 system_quadrupole.title = "System Quadrupole";
345 system_quadrupole.dataType = "Mat3x3d";
346 system_quadrupole.accumulator = new MatrixAccumulator();
347 data_[SYSTEM_QUADRUPOLE] = system_quadrupole;
348 statsMap_["SYSTEM_QUADRUPOLE"] = SYSTEM_QUADRUPOLE;
349
350 StatsData tagged_pair_distance;
351 tagged_pair_distance.units = "A";
352 tagged_pair_distance.title = "Tagged Pair Distance";
353 tagged_pair_distance.dataType = "RealType";
354 tagged_pair_distance.accumulator = new Accumulator();
355 data_[TAGGED_PAIR_DISTANCE] = tagged_pair_distance;
356 statsMap_["TAGGED_PAIR_DISTANCE"] = TAGGED_PAIR_DISTANCE;
357
358 StatsData shadowh;
359 shadowh.units = "kcal/mol";
360 shadowh.title = "Shadow Hamiltonian";
361 shadowh.dataType = "RealType";
362 shadowh.accumulator = new Accumulator();
363 data_[SHADOWH] = shadowh;
364 statsMap_["SHADOWH"] = SHADOWH;
365
366 StatsData helfandmoment;
367 helfandmoment.units = "A*kcal/mol";
368 helfandmoment.title = "Thermal Helfand Moment";
369 helfandmoment.dataType = "Vector3d";
370 helfandmoment.accumulator = new VectorAccumulator();
371 data_[HELFANDMOMENT] = helfandmoment;
372 statsMap_["HELFANDMOMENT"] = HELFANDMOMENT;
373
374 StatsData heatflux;
375 heatflux.units = "amu/fs^3";
376 heatflux.title = "Heat Flux";
377 heatflux.dataType = "Vector3d";
378 heatflux.accumulator = new VectorAccumulator();
379 data_[HEATFLUX] = heatflux;
380 statsMap_["HEATFLUX"] = HEATFLUX;
381
382 StatsData electronic_temperature;
383 electronic_temperature.units = "K";
384 electronic_temperature.title = "Electronic Temperature";
385 electronic_temperature.dataType = "RealType";
386 electronic_temperature.accumulator = new Accumulator();
387 data_[ELECTRONIC_TEMPERATURE] = electronic_temperature;
388 statsMap_["ELECTRONIC_TEMPERATURE"] = ELECTRONIC_TEMPERATURE;
389
390 StatsData com;
391 com.units = "A";
392 com.title = "Center of Mass";
393 com.dataType = "Vector3d";
394 com.accumulator = new VectorAccumulator();
395 data_[COM] = com;
396 statsMap_["COM"] = COM;
397
398 StatsData comVel;
399 comVel.units = "A/fs";
400 comVel.title = "COM Velocity";
401 comVel.dataType = "Vector3d";
402 comVel.accumulator = new VectorAccumulator();
403 data_[COM_VELOCITY] = comVel;
404 statsMap_["COM_VELOCITY"] = COM_VELOCITY;
405
406 StatsData angMom;
407 angMom.units = "amu A^2/fs";
408 angMom.title = "Angular Momentum";
409 angMom.dataType = "Vector3d";
410 angMom.accumulator = new VectorAccumulator();
411 data_[ANGULAR_MOMENTUM] = angMom;
412 statsMap_["ANGULAR_MOMENTUM"] = ANGULAR_MOMENTUM;
413
414 StatsData potSelection;
415 potSelection.units = "kcal/mol";
416 potSelection.title = "Selection Potentials";
417 potSelection.dataType = "potVec";
418 potSelection.accumulator = new PotVecAccumulator();
419 data_[POTENTIAL_SELECTION] = potSelection;
420 statsMap_["POTENTIAL_SELECTION"] = POTENTIAL_SELECTION;
421
422 StatsData netCharge;
423 netCharge.units = "e";
424 netCharge.title = "Net Charge";
425 netCharge.dataType = "RealType";
426 netCharge.accumulator = new Accumulator();
427 data_[NET_CHARGE] = netCharge;
428 statsMap_["NET_CHARGE"] = NET_CHARGE;
429
430 StatsData chargeMomentum;
431 chargeMomentum.units = "kcal fs / e / mol";
432 chargeMomentum.title = "Charge Momentum";
433 chargeMomentum.dataType = "RealType";
434 chargeMomentum.accumulator = new Accumulator();
435 data_[CHARGE_MOMENTUM] = chargeMomentum;
436 statsMap_["CHARGE_MOMENTUM"] = CHARGE_MOMENTUM;
437
438 StatsData currentDensity;
439 currentDensity.units = "amps / m^2";
440 currentDensity.title = "Current Density: total, then by AtomType";
441 currentDensity.dataType = "Array2d";
442 unsigned int nCols = 3 * (info_->getSimulatedAtomTypes().size()) + 3;
443 currentDensity.accumulatorArray2d.resize(nCols);
444 for (unsigned int j = 0; j < nCols; j++) {
445 currentDensity.accumulatorArray2d[j] = new Accumulator();
446 }
447 data_[CURRENT_DENSITY] = currentDensity;
448 statsMap_["CURRENT_DENSITY"] = CURRENT_DENSITY;
449
450 // Now, set some defaults in the mask:
451
452 Globals* simParams = info_->getSimParams();
453 std::string statFileFormatString = simParams->getStatFileFormat();
454 parseStatFileFormat(statFileFormatString);
455
456 // if we're doing a thermodynamic integration, we'll want the raw
457 // potential as well as the full potential:
458
459 if (simParams->getUseThermodynamicIntegration()) {
460 statsMask_.set(RAW_POTENTIAL);
461 statsMask_.set(RESTRAINT_POTENTIAL);
462 }
463
464 // if we've got restraints turned on, we'll also want a report of the
465 // total harmonic restraints
466 if (simParams->getUseRestraints()) {
467 statsMask_.set(RAW_POTENTIAL);
468 statsMask_.set(RESTRAINT_POTENTIAL);
469 }
470
471 if (simParams->havePrintPressureTensor() &&
472 simParams->getPrintPressureTensor()) {
473 statsMask_.set(PRESSURE_TENSOR);
474 }
475
476 if (simParams->havePrintVirialTensor() &&
477 simParams->getPrintVirialTensor()) {
478 statsMask_.set(VIRIAL_TENSOR);
479 }
480
481 // Why do we have both of these?
482 if (simParams->getAccumulateBoxDipole()) { statsMask_.set(SYSTEM_DIPOLE); }
483 if (info_->getCalcBoxDipole()) { statsMask_.set(SYSTEM_DIPOLE); }
484
485 // Why do we have both of these?
486 if (simParams->getAccumulateBoxQuadrupole()) {
487 statsMask_.set(SYSTEM_QUADRUPOLE);
488 }
489 if (info_->getCalcBoxQuadrupole()) { statsMask_.set(SYSTEM_QUADRUPOLE); }
490
491 if (simParams->havePrintHeatFlux()) {
492 if (simParams->getPrintHeatFlux()) { statsMask_.set(HEATFLUX); }
493 }
494
495 if (simParams->haveTaggedAtomPair() &&
496 simParams->havePrintTaggedPairDistance()) {
497 if (simParams->getPrintTaggedPairDistance()) {
498 statsMask_.set(TAGGED_PAIR_DISTANCE);
499 }
500 }
501
502 if (simParams->havePotentialSelection()) {
503 statsMask_.set(POTENTIAL_SELECTION);
504 }
505 }
506
507 int Stats::getPrecision() {
508 Globals* simParams = info_->getSimParams();
509 int statFilePrecision = simParams->getStatFilePrecision();
510 return statFilePrecision;
511 }
512
513 void Stats::parseStatFileFormat(const std::string& format) {
514 StringTokenizer tokenizer(format, " ,;|\t\n\r");
515
516 while (tokenizer.hasMoreTokens()) {
517 std::string token(tokenizer.nextToken());
518 toUpper(token);
519 StatsMapType::iterator i = statsMap_.find(token);
520 if (i != statsMap_.end()) {
521 statsMask_.set(i->second);
522 } else {
523 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
524 "Stats::parseStatFileFormat: %s is not a recognized\n"
525 "\tstatFileFormat keyword.\n",
526 token.c_str());
527 painCave.isFatal = 0;
528 painCave.severity = OPENMD_ERROR;
529 simError();
530 }
531 }
532 }
533
534 std::string Stats::getTitle(int index) {
535 assert(index >= 0 && index < ENDINDEX);
536 return data_[index].title;
537 }
538
539 std::string Stats::getUnits(int index) {
540 assert(index >= 0 && index < ENDINDEX);
541 return data_[index].units;
542 }
543
544 std::string Stats::getDataType(int index) {
545 assert(index >= 0 && index < ENDINDEX);
546 return data_[index].dataType;
547 }
548
549 void Stats::collectStats() {
550 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
551 Thermo thermo(info_);
552
553 for (unsigned int i = 0; i < statsMask_.size(); ++i) {
554 if (statsMask_[i]) {
555 switch (i) {
556 case TIME:
557 dynamic_cast<Accumulator*>(data_[i].accumulator)
558 ->add(snap->getTime());
559 break;
560 case KINETIC_ENERGY:
561 dynamic_cast<Accumulator*>(data_[i].accumulator)
562 ->add(thermo.getKinetic());
563 break;
564 case POTENTIAL_ENERGY:
565 dynamic_cast<Accumulator*>(data_[i].accumulator)
566 ->add(thermo.getPotential());
567 break;
568 case TOTAL_ENERGY:
569 dynamic_cast<Accumulator*>(data_[i].accumulator)
570 ->add(thermo.getTotalEnergy());
571 break;
572 case TEMPERATURE:
573 dynamic_cast<Accumulator*>(data_[i].accumulator)
574 ->add(thermo.getTemperature());
575 break;
576 case PRESSURE:
577 dynamic_cast<Accumulator*>(data_[i].accumulator)
578 ->add(thermo.getPressure());
579 break;
580 case VOLUME:
581 dynamic_cast<Accumulator*>(data_[i].accumulator)
582 ->add(thermo.getVolume());
583 break;
584 case CONSERVED_QUANTITY:
585 dynamic_cast<Accumulator*>(data_[i].accumulator)
586 ->add(snap->getConservedQuantity());
587 break;
588 case PRESSURE_TENSOR:
589 dynamic_cast<MatrixAccumulator*>(data_[i].accumulator)
590 ->add(thermo.getPressureTensor());
591 break;
592 // virial Tensor
593 case VIRIAL_TENSOR:
594 dynamic_cast<MatrixAccumulator*>(data_[i].accumulator)
595 ->add(snap->getVirialTensor());
596 break;
597
598 case SYSTEM_DIPOLE:
599 dynamic_cast<VectorAccumulator*>(data_[i].accumulator)
600 ->add(thermo.getSystemDipole());
601 break;
602 case SYSTEM_QUADRUPOLE:
603 dynamic_cast<MatrixAccumulator*>(data_[i].accumulator)
604 ->add(thermo.getSystemQuadrupole());
605 break;
606 case HEATFLUX:
607 dynamic_cast<VectorAccumulator*>(data_[i].accumulator)
608 ->add(thermo.getHeatFlux());
609 break;
610 case HULLVOLUME:
611 dynamic_cast<Accumulator*>(data_[i].accumulator)
612 ->add(thermo.getHullVolume());
613 break;
614 case GYRVOLUME:
615 dynamic_cast<Accumulator*>(data_[i].accumulator)
616 ->add(thermo.getGyrationalVolume());
617 break;
618 case TRANSLATIONAL_KINETIC:
619 dynamic_cast<Accumulator*>(data_[i].accumulator)
620 ->add(thermo.getTranslationalKinetic());
621 break;
622 case ROTATIONAL_KINETIC:
623 dynamic_cast<Accumulator*>(data_[i].accumulator)
624 ->add(thermo.getRotationalKinetic());
625 break;
626 case ELECTRONIC_KINETIC:
627 dynamic_cast<Accumulator*>(data_[i].accumulator)
628 ->add(thermo.getElectronicKinetic());
629 break;
630 case LONG_RANGE_POTENTIAL:
631 dynamic_cast<Accumulator*>(data_[i].accumulator)
632 ->add(snap->getLongRangePotential());
633 break;
634 case VANDERWAALS_POTENTIAL:
635 dynamic_cast<Accumulator*>(data_[i].accumulator)
636 ->add(snap->getLongRangePotentials()[VANDERWAALS_FAMILY]);
637 break;
638 case ELECTROSTATIC_POTENTIAL:
639 dynamic_cast<Accumulator*>(data_[i].accumulator)
640 ->add(snap->getLongRangePotentials()[ELECTROSTATIC_FAMILY] +
641 snap->getSelfPotentials()[ELECTROSTATIC_FAMILY]);
642 break;
643 case METALLIC_POTENTIAL:
644 dynamic_cast<Accumulator*>(data_[i].accumulator)
645 ->add(snap->getSelfPotentials()[METALLIC_EMBEDDING_FAMILY] +
646 snap->getLongRangePotentials()[METALLIC_PAIR_FAMILY]);
647 break;
648 case METALLIC_EMBEDDING:
649 dynamic_cast<Accumulator*>(data_[i].accumulator)
650 ->add(snap->getSelfPotentials()[METALLIC_EMBEDDING_FAMILY]);
651 break;
652 case METALLIC_PAIR:
653 dynamic_cast<Accumulator*>(data_[i].accumulator)
654 ->add(snap->getLongRangePotentials()[METALLIC_PAIR_FAMILY]);
655 break;
656 case HYDROGENBONDING_POTENTIAL:
657 dynamic_cast<Accumulator*>(data_[i].accumulator)
658 ->add(snap->getLongRangePotentials()[HYDROGENBONDING_FAMILY]);
659 break;
660 case RECIPROCAL_POTENTIAL:
661 dynamic_cast<Accumulator*>(data_[i].accumulator)
662 ->add(snap->getReciprocalPotential());
663 break;
664 case SURFACE_POTENTIAL:
665 dynamic_cast<Accumulator*>(data_[i].accumulator)
666 ->add(snap->getSurfacePotential());
667 break;
668 case SHORT_RANGE_POTENTIAL:
669 dynamic_cast<Accumulator*>(data_[i].accumulator)
670 ->add(snap->getShortRangePotential());
671 break;
672 case BOND_POTENTIAL:
673 dynamic_cast<Accumulator*>(data_[i].accumulator)
674 ->add(snap->getBondPotential());
675 break;
676 case BEND_POTENTIAL:
677 dynamic_cast<Accumulator*>(data_[i].accumulator)
678 ->add(snap->getBendPotential());
679 break;
680 case DIHEDRAL_POTENTIAL:
681 dynamic_cast<Accumulator*>(data_[i].accumulator)
682 ->add(snap->getTorsionPotential());
683 break;
684 case INVERSION_POTENTIAL:
685 dynamic_cast<Accumulator*>(data_[i].accumulator)
686 ->add(snap->getInversionPotential());
687 break;
688 case RAW_POTENTIAL:
689 dynamic_cast<Accumulator*>(data_[i].accumulator)
690 ->add(snap->getRawPotential());
691 break;
692 case RESTRAINT_POTENTIAL:
693 dynamic_cast<Accumulator*>(data_[i].accumulator)
694 ->add(snap->getRestraintPotential());
695 break;
696 case EXCLUDED_POTENTIAL:
697 dynamic_cast<Accumulator*>(data_[i].accumulator)
698 ->add(snap->getExcludedPotential());
699 break;
700 case TAGGED_PAIR_DISTANCE:
701 dynamic_cast<Accumulator*>(data_[i].accumulator)
702 ->add(thermo.getTaggedAtomPairDistance());
703 break;
704 case ELECTRONIC_TEMPERATURE:
705 dynamic_cast<Accumulator*>(data_[i].accumulator)
706 ->add(thermo.getElectronicTemperature());
707 break;
708 case COM:
709 dynamic_cast<VectorAccumulator*>(data_[i].accumulator)
710 ->add(thermo.getCom());
711 break;
712 case COM_VELOCITY:
713 dynamic_cast<VectorAccumulator*>(data_[i].accumulator)
714 ->add(thermo.getComVel());
715 break;
716 case ANGULAR_MOMENTUM:
717 dynamic_cast<VectorAccumulator*>(data_[i].accumulator)
718 ->add(thermo.getAngularMomentum());
719 break;
720 case POTENTIAL_SELECTION:
721 dynamic_cast<PotVecAccumulator*>(data_[i].accumulator)
722 ->add(thermo.getSelectionPotentials());
723 break;
724 case NET_CHARGE:
725 dynamic_cast<Accumulator*>(data_[i].accumulator)
726 ->add(thermo.getNetCharge());
727 break;
728 case CHARGE_MOMENTUM:
729 dynamic_cast<Accumulator*>(data_[i].accumulator)
730 ->add(thermo.getChargeMomentum());
731 break;
732 case CURRENT_DENSITY:
733 unsigned int k = 0;
734 std::vector<Vector3d> Jc = thermo.getCurrentDensity();
735 for (unsigned int j = 0; j < Jc.size(); ++j) {
736 dynamic_cast<Accumulator*>(data_[i].accumulatorArray2d[k++])
737 ->add(Jc[j][0]);
738 dynamic_cast<Accumulator*>(data_[i].accumulatorArray2d[k++])
739 ->add(Jc[j][1]);
740 dynamic_cast<Accumulator*>(data_[i].accumulatorArray2d[k++])
741 ->add(Jc[j][2]);
742 }
743 break;
744
745 /*
746 case SHADOWH:
747 dynamic_cast<Accumulator
748 *>(data_[i].accumulator)->add(thermo.getShadowHamiltionian());
749 break; case HELFANDMOMENT: dynamic_cast<Accumulator
750 *>(data_[i].accumulator)->add(thermo.getHelfandMoment()); break;
751 */
752 }
753 }
754 }
755 }
756
757 int Stats::getIntData(int index) {
758 assert(index >= 0 && index < ENDINDEX);
759 RealType value;
760 dynamic_cast<Accumulator*>(data_[index].accumulator)->getLastValue(value);
761 return (int)value;
762 }
763 RealType Stats::getRealData(int index) {
764 assert(index >= 0 && index < ENDINDEX);
765 RealType value(0.0);
766 dynamic_cast<Accumulator*>(data_[index].accumulator)->getLastValue(value);
767 return value;
768 }
769 Vector3d Stats::getVectorData(int index) {
770 assert(index >= 0 && index < ENDINDEX);
771 Vector3d value;
772 dynamic_cast<VectorAccumulator*>(data_[index].accumulator)
773 ->getLastValue(value);
774 return value;
775 }
776 potVec Stats::getPotVecData(int index) {
777 assert(index >= 0 && index < ENDINDEX);
778 potVec value;
779 dynamic_cast<PotVecAccumulator*>(data_[index].accumulator)
780 ->getLastValue(value);
781 return value;
782 }
783 Mat3x3d Stats::getMatrixData(int index) {
784 assert(index >= 0 && index < ENDINDEX);
785 Mat3x3d value;
786 dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)
787 ->getLastValue(value);
788 return value;
789 }
790 std::vector<RealType> Stats::getArrayData(int index) {
791 assert(index >= 0 && index < ENDINDEX);
792 std::vector<RealType> value;
793 RealType v;
794 for (unsigned int i = 0; i < data_[index].accumulatorArray2d.size(); ++i) {
795 dynamic_cast<Accumulator*>(data_[index].accumulatorArray2d[i])
796 ->getLastValue(v);
797 value.push_back(v);
798 }
799 return value;
800 }
801
802 int Stats::getIntAverage(int index) {
803 assert(index >= 0 && index < ENDINDEX);
804 RealType value;
805 dynamic_cast<Accumulator*>(data_[index].accumulator)->getAverage(value);
806 return (int)value;
807 }
808 RealType Stats::getRealAverage(int index) {
809 assert(index >= 0 && index < ENDINDEX);
810 RealType value(0.0);
811 dynamic_cast<Accumulator*>(data_[index].accumulator)->getAverage(value);
812 return value;
813 }
814 Vector3d Stats::getVectorAverage(int index) {
815 assert(index >= 0 && index < ENDINDEX);
816 Vector3d value;
817 dynamic_cast<VectorAccumulator*>(data_[index].accumulator)
818 ->getAverage(value);
819 return value;
820 }
821 potVec Stats::getPotVecAverage(int index) {
822 assert(index >= 0 && index < ENDINDEX);
823 potVec value;
824 dynamic_cast<PotVecAccumulator*>(data_[index].accumulator)
825 ->getAverage(value);
826 return value;
827 }
828 Mat3x3d Stats::getMatrixAverage(int index) {
829 assert(index >= 0 && index < ENDINDEX);
830 Mat3x3d value;
831 dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)
832 ->getAverage(value);
833 return value;
834 }
835 std::vector<RealType> Stats::getArrayAverage(int index) {
836 assert(index >= 0 && index < ENDINDEX);
837 std::vector<RealType> value;
838 RealType v;
839 for (unsigned int i = 0; i < data_[index].accumulatorArray2d.size(); ++i) {
840 dynamic_cast<Accumulator*>(data_[index].accumulatorArray2d[i])
841 ->getAverage(v);
842 value.push_back(v);
843 }
844 return value;
845 }
846
847 int Stats::getIntError(int index) {
848 assert(index >= 0 && index < ENDINDEX);
849 RealType value;
850 dynamic_cast<Accumulator*>(data_[index].accumulator)
851 ->get95percentConfidenceInterval(value);
852 return (int)value;
853 }
854 RealType Stats::getRealError(int index) {
855 assert(index >= 0 && index < ENDINDEX);
856 RealType value(0.0);
857 dynamic_cast<Accumulator*>(data_[index].accumulator)
858 ->get95percentConfidenceInterval(value);
859 return value;
860 }
861 Vector3d Stats::getVectorError(int index) {
862 assert(index >= 0 && index < ENDINDEX);
863 Vector3d value;
864 dynamic_cast<VectorAccumulator*>(data_[index].accumulator)
865 ->get95percentConfidenceInterval(value);
866 return value;
867 }
868 potVec Stats::getPotVecError(int index) {
869 assert(index >= 0 && index < ENDINDEX);
870 potVec value;
871 dynamic_cast<PotVecAccumulator*>(data_[index].accumulator)
872 ->get95percentConfidenceInterval(value);
873 return value;
874 }
875 Mat3x3d Stats::getMatrixError(int index) {
876 assert(index >= 0 && index < ENDINDEX);
877 Mat3x3d value;
878 dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)
879 ->get95percentConfidenceInterval(value);
880 return value;
881 }
882 std::vector<RealType> Stats::getArrayError(int index) {
883 assert(index >= 0 && index < ENDINDEX);
884 std::vector<RealType> value;
885 RealType v;
886 for (unsigned int i = 0; i < data_[index].accumulatorArray2d.size(); ++i) {
887 dynamic_cast<Accumulator*>(data_[index].accumulatorArray2d[i])
888 ->get95percentConfidenceInterval(v);
889 value.push_back(v);
890 }
891 return value;
892 }
893
894 Stats::StatsBitSet Stats::getStatsMask() { return statsMask_; }
895 Stats::StatsMapType Stats::getStatsMap() { return statsMap_; }
896 void Stats::setStatsMask(Stats::StatsBitSet mask) { statsMask_ = mask; }
897
898 std::string Stats::getStatsReport() {
899 std::stringstream report;
900#if defined(_MSC_VER)
901 std::string pm = " +/- ";
902 std::string luc = "[";
903 std::string lex = "|";
904 std::string llc = "[";
905 std::string ruc = "]";
906 std::string rex = "|";
907 std::string rlc = "]";
908#else
909 std::string pm = " \u00B1 ";
910 std::string luc = "\u23A1";
911 std::string lex = "\u23A2";
912 std::string llc = "\u23A3";
913 std::string ruc = "\u23A4";
914 std::string rex = "\u23A5";
915 std::string rlc = "\u23A6";
916#endif
917
918 int nSamp = dynamic_cast<Accumulator*>(data_[TIME].accumulator)->count();
919
920 std::string head(79, '#');
921 report << head << std::endl;
922 report << "# Status Report:" << std::string(62, ' ') << "#" << std::endl;
923 report << "# " << right << setw(24) << "Total Time:";
924 report << setw(12) << getRealData(TIME);
925 report << " " << setw(17) << left << getUnits(TIME)
926 << " #" << std::endl;
927 report << "# " << right << setw(24) << "Number of Samples:";
928 report << setw(12) << nSamp;
929 report << " #" << std::endl;
930
931 for (unsigned int i = 0; i < statsMask_.size(); ++i) {
932 if (statsMask_[i] && i != TIME) {
933 if (getDataType(i) == "RealType") {
934 report << "# " << right << setw(23) << getTitle(i) << ":";
935 report << right << setw(12) << getRealAverage(i);
936 report << pm << left << setw(12) << getRealError(i);
937 report << " " << left << setw(17) << getUnits(i);
938 report << " #" << std::endl;
939
940 } else if (getDataType(i) == "Vector3d") {
941 Vector3d s = getVectorAverage(i);
942 Vector3d e = getVectorError(i);
943
944 report << "# ";
945 report << luc << right << setw(12) << s(0) << ruc << " ";
946 report << luc << right << setw(12) << e(0);
947 report << ruc << " #" << std::endl;
948
949 report << "# " << right << setw(23) << getTitle(i) << ":";
950
951 report << lex << right << setw(12) << s(1) << rex << pm;
952 report << lex << right << setw(12) << e(1) << rex << " ";
953 report << left << setw(17) << getUnits(i) << " #";
954 report << std::endl;
955
956 report << "# ";
957 report << llc << right << setw(12) << s(2) << rlc << " ";
958 report << llc << right << setw(12) << e(2);
959 report << rlc << " #" << std::endl;
960
961 } else if (getDataType(i) == "potVec") {
962 potVec s = getPotVecAverage(i);
963 potVec e = getPotVecError(i);
964
965 report << "# " << right << setw(23) << getTitle(i);
966 report << ": #";
967 report << std::endl;
968
969 for (unsigned int j = 1; j < N_INTERACTION_FAMILIES; j++) {
970 switch (j) {
971 case VANDERWAALS_FAMILY:
972 report << "# " << right << setw(24) << "van der Waals:";
973 break;
974 case ELECTROSTATIC_FAMILY:
975 report << "# " << right << setw(24) << "Electrostatic:";
976 break;
977 case METALLIC_EMBEDDING:
978 report << "# " << right << setw(24) << "Metallic Embedding:";
979 break;
980 case METALLIC_PAIR:
981 report << "# " << right << setw(24) << "Metallic Pair:";
982 break;
983 case HYDROGENBONDING_FAMILY:
984 report << "# " << right << setw(24) << "Hydrogen Bonding:";
985 break;
986 case BONDED_FAMILY:
987 report << "# " << right << setw(24) << "Bonded (1-2,1-3,1-4):";
988 break;
989 default:
990 report << "# " << right << setw(24) << "Unknown:";
991 break;
992 }
993 report << right << setw(12) << s[j];
994 report << pm << left << setw(12) << e[j];
995 report << " " << left << setw(17) << getUnits(i) << " #";
996 report << std::endl;
997 }
998
999 } else if (getDataType(i) == "Mat3x3d") {
1000 Mat3x3d s = getMatrixAverage(i);
1001 Mat3x3d e = getMatrixError(i);
1002
1003 report << "# ";
1004 report << luc << right << setw(12) << s(0, 0) << " ";
1005 report << right << setw(12) << s(0, 1) << " ";
1006 report << right << setw(12) << s(0, 2) << ruc << " #";
1007 report << std::endl;
1008
1009 report << "# " << right << setw(23) << getTitle(i) << ":";
1010
1011 report << lex << right << setw(12) << s(1, 0) << " ";
1012 report << right << setw(12) << s(1, 1) << " ";
1013 report << right << setw(12) << s(1, 2) << rex << " ";
1014 report << left << setw(11) << getUnits(i) << "#" << std::endl;
1015
1016 report << "# ";
1017 report << llc << right << setw(12) << s(2, 0) << " ";
1018 report << right << setw(12) << s(2, 1) << " ";
1019 report << right << setw(12) << s(2, 2) << rlc << " #";
1020 report << std::endl;
1021
1022 report << "# ";
1023 report << luc << right << setw(12) << e(0, 0) << " ";
1024 report << right << setw(12) << e(0, 1) << " ";
1025 report << right << setw(12) << e(0, 2) << ruc << " #";
1026 report << std::endl;
1027
1028 report << "# " << pm;
1029
1030 report << lex << right << setw(12) << e(1, 0) << " ";
1031 report << right << setw(12) << e(1, 1) << " ";
1032 report << right << setw(12) << e(1, 2) << rex << " #";
1033 report << std::endl;
1034
1035 report << "# ";
1036 report << llc << right << setw(12) << e(2, 0) << " ";
1037 report << right << setw(12) << e(2, 1) << " ";
1038 report << right << setw(12) << e(2, 2) << rlc << " #";
1039 report << std::endl;
1040 }
1041 }
1042 }
1043 report << head << std::endl;
1044
1045 return report.str();
1046 }
1047
1048} // 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
OpenMD::HYDROGENBONDING_FAMILY
@ HYDROGENBONDING_FAMILY
Short-range directional interactions.
Definition NonBondedInteraction.hpp:74
OpenMD::VANDERWAALS_FAMILY
@ VANDERWAALS_FAMILY
Long-range dispersion and short-range pauli repulsion.
Definition NonBondedInteraction.hpp:67
OpenMD::BONDED_FAMILY
@ BONDED_FAMILY
directly bonded 1-2, 1-3, or 1-4 interactions
Definition NonBondedInteraction.hpp:75
OpenMD::ELECTROSTATIC_FAMILY
@ ELECTROSTATIC_FAMILY
Coulombic and point-multipole interactions.
Definition NonBondedInteraction.hpp:69