ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Stats.cpp
Revision: 1764
Committed: Tue Jul 3 18:32:27 2012 UTC (12 years, 9 months ago) by gezelter
File size: 20328 byte(s)
Log Message:
Refactored Snapshot and Stats to use the Accumulator classes.  Collected
a number of methods into Thermo that belonged there.

File Contents

# Content
1 /*
2 * Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 */
42
43 /**
44 * @file Stats.cpp
45 * @author tlin
46 * @date 11/04/2004
47 * @time 14:26am
48 * @version 1.0
49 */
50
51 #include "brains/Stats.hpp"
52 #include "brains/Thermo.hpp"
53
54 namespace OpenMD {
55
56 Stats::Stats(SimInfo* info) : isInit_(false), info_(info) {
57
58 if (!isInit_) {
59 init();
60 isInit_ = true;
61 }
62 }
63
64 void Stats::init() {
65
66 data_.resize(Stats::ENDINDEX);
67
68 StatsData time;
69 time.units = "fs";
70 time.title = "Time";
71 time.dataType = "RealType";
72 time.accumulator = new Accumulator();
73 data_[TIME] = time;
74 statsMap_["TIME"] = TIME;
75
76 StatsData total_energy;
77 total_energy.units = "kcal/mol";
78 total_energy.title = "Total Energy";
79 total_energy.dataType = "RealType";
80 total_energy.accumulator = new Accumulator();
81 data_[TOTAL_ENERGY] = total_energy;
82 statsMap_["TOTAL_ENERGY"] = TOTAL_ENERGY;
83
84 StatsData potential_energy;
85 potential_energy.units = "kcal/mol";
86 potential_energy.title = "Potential Energy";
87 potential_energy.dataType = "RealType";
88 potential_energy.accumulator = new Accumulator();
89 data_[POTENTIAL_ENERGY] = potential_energy;
90 statsMap_["POTENTIAL_ENERGY"] = POTENTIAL_ENERGY;
91
92 StatsData kinetic_energy;
93 kinetic_energy.units = "kcal/mol";
94 kinetic_energy.title = "Kinetic Energy";
95 kinetic_energy.dataType = "RealType";
96 kinetic_energy.accumulator = new Accumulator();
97 data_[KINETIC_ENERGY] = kinetic_energy;
98 statsMap_["KINETIC_ENERGY"] = KINETIC_ENERGY;
99
100 StatsData temperature;
101 temperature.units = "K";
102 temperature.title = "Temperature";
103 temperature.dataType = "RealType";
104 temperature.accumulator = new Accumulator();
105 data_[TEMPERATURE] = temperature;
106 statsMap_["TEMPERATURE"] = TEMPERATURE;
107
108 StatsData pressure;
109 pressure.units = "atm";
110 pressure.title = "Pressure";
111 pressure.dataType = "RealType";
112 pressure.accumulator = new Accumulator();
113 data_[PRESSURE] = pressure;
114 statsMap_["PRESSURE"] = PRESSURE;
115
116 StatsData volume;
117 volume.units = "A^3";
118 volume.title = "Volume";
119 volume.dataType = "RealType";
120 volume.accumulator = new Accumulator();
121 data_[VOLUME] = volume;
122 statsMap_["VOLUME"] = VOLUME;
123
124 StatsData hullvolume;
125 hullvolume.units = "A^3";
126 hullvolume.title = "Hull Volume";
127 hullvolume.dataType = "RealType";
128 hullvolume.accumulator = new Accumulator();
129 data_[HULLVOLUME] = hullvolume;
130 statsMap_["HULLVOLUME"] = HULLVOLUME;
131
132 StatsData gyrvolume;
133 gyrvolume.units = "A^3";
134 gyrvolume.title = "Gyrational Volume";
135 gyrvolume.dataType = "RealType";
136 gyrvolume.accumulator = new Accumulator();
137 data_[GYRVOLUME] = gyrvolume;
138 statsMap_["GYRVOLUME"] = GYRVOLUME;
139
140 StatsData conserved_quantity;
141 conserved_quantity.units = "kcal/mol";
142 conserved_quantity.title = "Conserved Quantity";
143 conserved_quantity.dataType = "RealType";
144 conserved_quantity.accumulator = new Accumulator();
145 data_[CONSERVED_QUANTITY] = conserved_quantity;
146 statsMap_["CONSERVED_QUANTITY"] = CONSERVED_QUANTITY;
147
148 StatsData translational_kinetic;
149 translational_kinetic.units = "kcal/mol";
150 translational_kinetic.title = "Translational Kinetic";
151 translational_kinetic.dataType = "RealType";
152 translational_kinetic.accumulator = new Accumulator();
153 data_[TRANSLATIONAL_KINETIC] = translational_kinetic;
154 statsMap_["TRANSLATIONAL_KINETIC"] = TRANSLATIONAL_KINETIC;
155
156 StatsData rotational_kinetic;
157 rotational_kinetic.units = "kcal/mol";
158 rotational_kinetic.title = "Rotational Kinetic";
159 rotational_kinetic.dataType = "RealType";
160 rotational_kinetic.accumulator = new Accumulator();
161 data_[ROTATIONAL_KINETIC] = rotational_kinetic;
162 statsMap_["ROTATIONAL_KINETIC"] = ROTATIONAL_KINETIC;
163
164 StatsData long_range_potential;
165 long_range_potential.units = "kcal/mol";
166 long_range_potential.title = "Long Range Potential";
167 long_range_potential.dataType = "RealType";
168 long_range_potential.accumulator = new Accumulator();
169 data_[LONG_RANGE_POTENTIAL] = long_range_potential;
170 statsMap_["LONG_RANGE_POTENTIAL"] = LONG_RANGE_POTENTIAL;
171
172 StatsData vanderwaals_potential;
173 vanderwaals_potential.units = "kcal/mol";
174 vanderwaals_potential.title = "van der waals Potential";
175 vanderwaals_potential.dataType = "RealType";
176 vanderwaals_potential.accumulator = new Accumulator();
177 data_[VANDERWAALS_POTENTIAL] = vanderwaals_potential;
178 statsMap_["VANDERWAALS_POTENTIAL"] = VANDERWAALS_POTENTIAL;
179
180 StatsData electrostatic_potential;
181 electrostatic_potential.units = "kcal/mol";
182 electrostatic_potential.title = "Electrostatic Potential";
183 electrostatic_potential.dataType = "RealType";
184 electrostatic_potential.accumulator = new Accumulator();
185 data_[ELECTROSTATIC_POTENTIAL] = electrostatic_potential;
186 statsMap_["ELECTROSTATIC_POTENTIAL"] = ELECTROSTATIC_POTENTIAL;
187
188 StatsData metallic_potential;
189 metallic_potential.units = "kcal/mol";
190 metallic_potential.title = "Metallic Potential";
191 metallic_potential.dataType = "RealType";
192 metallic_potential.accumulator = new Accumulator();
193 data_[METALLIC_POTENTIAL] = metallic_potential;
194 statsMap_["METALLIC_POTENTIAL"] = METALLIC_POTENTIAL;
195
196 StatsData hydrogenbonding_potential;
197 hydrogenbonding_potential.units = "kcal/mol";
198 hydrogenbonding_potential.title = "Metallic Potential";
199 hydrogenbonding_potential.dataType = "RealType";
200 hydrogenbonding_potential.accumulator = new Accumulator();
201 data_[HYDROGENBONDING_POTENTIAL] = hydrogenbonding_potential;
202 statsMap_["HYDROGENBONDING_POTENTIAL"] = HYDROGENBONDING_POTENTIAL;
203
204 StatsData short_range_potential;
205 short_range_potential.units = "kcal/mol";
206 short_range_potential.title = "Short Range Potential";
207 short_range_potential.dataType = "RealType";
208 short_range_potential.accumulator = new Accumulator();
209 data_[SHORT_RANGE_POTENTIAL] = short_range_potential;
210 statsMap_["SHORT_RANGE_POTENTIAL"] = SHORT_RANGE_POTENTIAL;
211
212 StatsData bond_potential;
213 bond_potential.units = "kcal/mol";
214 bond_potential.title = "Bond Potential";
215 bond_potential.dataType = "RealType";
216 bond_potential.accumulator = new Accumulator();
217 data_[BOND_POTENTIAL] = bond_potential;
218 statsMap_["BOND_POTENTIAL"] = BOND_POTENTIAL;
219
220 StatsData bend_potential;
221 bend_potential.units = "kcal/mol";
222 bend_potential.title = "Bend Potential";
223 bend_potential.dataType = "RealType";
224 bend_potential.accumulator = new Accumulator();
225 data_[BEND_POTENTIAL] = bend_potential;
226 statsMap_["BEND_POTENTIAL"] = BEND_POTENTIAL;
227
228 StatsData dihedral_potential;
229 dihedral_potential.units = "kcal/mol";
230 dihedral_potential.title = "Dihedral Potential";
231 dihedral_potential.dataType = "RealType";
232 dihedral_potential.accumulator = new Accumulator();
233 data_[DIHEDRAL_POTENTIAL] = dihedral_potential;
234 statsMap_["DIHEDRAL_POTENTIAL"] = DIHEDRAL_POTENTIAL;
235
236 StatsData inversion_potential;
237 inversion_potential.units = "kcal/mol";
238 inversion_potential.title = "Inversion Potential";
239 inversion_potential.dataType = "RealType";
240 inversion_potential.accumulator = new Accumulator();
241 data_[INVERSION_POTENTIAL] = inversion_potential;
242 statsMap_["INVERSION_POTENTIAL"] = INVERSION_POTENTIAL;
243
244 StatsData vraw;
245 vraw.units = "kcal/mol";
246 vraw.title = "Raw Potential";
247 vraw.dataType = "RealType";
248 vraw.accumulator = new Accumulator();
249 data_[RAW_POTENTIAL] = vraw;
250 statsMap_["RAW_POTENTIAL"] = RAW_POTENTIAL;
251
252 StatsData vrestraint;
253 vrestraint.units = "kcal/mol";
254 vrestraint.title = "Restraint Potential";
255 vrestraint.dataType = "RealType";
256 vrestraint.accumulator = new Accumulator();
257 data_[RESTRAINT_POTENTIAL] = vrestraint;
258 statsMap_["RESTRAINT_POTENTIAL"] = RESTRAINT_POTENTIAL;
259
260 StatsData pressure_tensor;
261 pressure_tensor.units = "amu*fs^-2*Ang^-1";
262 pressure_tensor.title = "Ptensor";
263 pressure_tensor.dataType = "Mat3x3d";
264 pressure_tensor.accumulator = new MatrixAccumulator();
265 data_[PRESSURE_TENSOR] = pressure_tensor;
266 statsMap_["PRESSURE_TENSOR"] = PRESSURE_TENSOR;
267
268 StatsData system_dipole;
269 system_dipole.units = "C*m";
270 system_dipole.title = "System Dipole";
271 system_dipole.dataType = "Vector3d";
272 system_dipole.accumulator = new VectorAccumulator();
273 data_[SYSTEM_DIPOLE] = system_dipole;
274 statsMap_["SYSTEM_DIPOLE"] = SYSTEM_DIPOLE;
275
276 StatsData tagged_pair_distance;
277 tagged_pair_distance.units = "Ang";
278 tagged_pair_distance.title = "Tagged_Pair_Distance";
279 tagged_pair_distance.dataType = "RealType";
280 tagged_pair_distance.accumulator = new Accumulator();
281 data_[TAGGED_PAIR_DISTANCE] = tagged_pair_distance;
282 statsMap_["TAGGED_PAIR_DISTANCE"] = TAGGED_PAIR_DISTANCE;
283
284 StatsData rnemd_exchange_total;
285 rnemd_exchange_total.units = "Variable";
286 rnemd_exchange_total.title = "RNEMD_exchange_total";
287 rnemd_exchange_total.dataType = "RealType";
288 rnemd_exchange_total.accumulator = new Accumulator();
289 data_[RNEMD_EXCHANGE_TOTAL] = rnemd_exchange_total;
290 statsMap_["RNEMD_EXCHANGE_TOTAL"] = RNEMD_EXCHANGE_TOTAL;
291
292 StatsData shadowh;
293 shadowh.units = "kcal/mol";
294 shadowh.title = "Shadow Hamiltonian";
295 shadowh.dataType = "RealType";
296 shadowh.accumulator = new Accumulator();
297 data_[SHADOWH] = shadowh;
298 statsMap_["SHADOWH"] = SHADOWH;
299
300 StatsData helfandmoment;
301 helfandmoment.units = "Ang*kcal/mol";
302 helfandmoment.title = "Thermal Helfand Moment";
303 helfandmoment.dataType = "Vector3d";
304 helfandmoment.accumulator = new VectorAccumulator();
305 data_[HELFANDMOMENT] = helfandmoment;
306 statsMap_["HELFANDMOMENT"] = HELFANDMOMENT;
307
308 StatsData heatflux;
309 heatflux.units = "amu/fs^3";
310 heatflux.title = "Heat Flux";
311 heatflux.dataType = "Vector3d";
312 heatflux.accumulator = new VectorAccumulator();
313 data_[HEATFLUX] = heatflux;
314 statsMap_["HEATFLUX"] = HEATFLUX;
315
316 StatsData electronic_temperature;
317 electronic_temperature.units = "K";
318 electronic_temperature.title = "Electronic Temperature";
319 electronic_temperature.dataType = "RealType";
320 electronic_temperature.accumulator = new Accumulator();
321 data_[ELECTRONIC_TEMPERATURE] = electronic_temperature;
322 statsMap_["ELECTRONIC_TEMPERATURE"] = ELECTRONIC_TEMPERATURE;
323
324 // Now, set some defaults in the mask:
325
326 Globals* simParams = info_->getSimParams();
327 std::string statFileFormatString = simParams->getStatFileFormat();
328 parseStatFileFormat(statFileFormatString);
329
330 // if we're doing a thermodynamic integration, we'll want the raw
331 // potential as well as the full potential:
332
333 if (simParams->getUseThermodynamicIntegration())
334 statsMask_.set(RAW_POTENTIAL);
335
336 // if we've got restraints turned on, we'll also want a report of the
337 // total harmonic restraints
338 if (simParams->getUseRestraints()){
339 statsMask_.set(RESTRAINT_POTENTIAL);
340 }
341
342 if (simParams->havePrintPressureTensor() &&
343 simParams->getPrintPressureTensor()){
344 statsMask_.set(PRESSURE_TENSOR);
345 }
346
347 // Why do we have both of these?
348 if (simParams->getAccumulateBoxDipole()) {
349 statsMask_.set(SYSTEM_DIPOLE);
350 }
351 if (info_->getCalcBoxDipole()){
352 statsMask_.set(SYSTEM_DIPOLE);
353 }
354
355 if (simParams->havePrintHeatFlux()) {
356 if (simParams->getPrintHeatFlux()){
357 statsMask_.set(HEATFLUX);
358 }
359 }
360
361
362 if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
363 if (simParams->getPrintTaggedPairDistance()) {
364 statsMask_.set(TAGGED_PAIR_DISTANCE);
365 }
366 }
367
368 if (simParams->getRNEMDParameters()->getUseRNEMD()) {
369 statsMask_.set(RNEMD_EXCHANGE_TOTAL);
370 }
371
372
373
374 }
375
376 void Stats::parseStatFileFormat(const std::string& format) {
377 StringTokenizer tokenizer(format, " ,;|\t\n\r");
378
379 while(tokenizer.hasMoreTokens()) {
380 std::string token(tokenizer.nextToken());
381 toUpper(token);
382 StatsMapType::iterator i = statsMap_.find(token);
383 if (i != statsMap_.end()) {
384 statsMask_.set(i->second);
385 } else {
386 sprintf( painCave.errMsg,
387 "Stats::parseStatFileFormat: %s is not a recognized\n"
388 "\tstatFileFormat keyword.\n", token.c_str() );
389 painCave.isFatal = 0;
390 painCave.severity = OPENMD_ERROR;
391 simError();
392 }
393 }
394 }
395
396
397 std::string Stats::getTitle(int index) {
398 assert(index >=0 && index < ENDINDEX);
399 return data_[index].title;
400 }
401
402 std::string Stats::getUnits(int index) {
403 assert(index >=0 && index < ENDINDEX);
404 return data_[index].units;
405 }
406
407 std::string Stats::getDataType(int index) {
408 assert(index >=0 && index < ENDINDEX);
409 return data_[index].dataType;
410 }
411
412 void Stats::collectStats(){
413 Globals* simParams = info_->getSimParams();
414 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
415 Thermo thermo(info_);
416
417 for (int i = 0; i < statsMask_.size(); ++i) {
418 if (statsMask_[i]) {
419 switch (i) {
420 case TIME:
421 data_[i].accumulator->add(snap->getTime());
422 break;
423 case KINETIC_ENERGY:
424 data_[i].accumulator->add(thermo.getKinetic());
425 break;
426 case POTENTIAL_ENERGY:
427 data_[i].accumulator->add(thermo.getPotential());
428 break;
429 case TOTAL_ENERGY:
430 data_[i].accumulator->add(thermo.getTotalEnergy());
431 break;
432 case TEMPERATURE:
433 data_[i].accumulator->add(thermo.getTemperature());
434 break;
435 case PRESSURE:
436 data_[i].accumulator->add(thermo.getPressure());
437 break;
438 case VOLUME:
439 data_[i].accumulator->add(thermo.getVolume());
440 break;
441 case CONSERVED_QUANTITY:
442 data_[i].accumulator->add(snap->getConservedQuantity());
443 break;
444 case PRESSURE_TENSOR:
445 dynamic_cast<MatrixAccumulator *>(data_[i].accumulator)->add(thermo.getPressureTensor());
446 break;
447 case SYSTEM_DIPOLE:
448 dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getSystemDipole());
449 break;
450 case HEATFLUX:
451 dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getHeatFlux());
452 break;
453 case HULLVOLUME:
454 data_[i].accumulator->add(thermo.getHullVolume());
455 break;
456 case GYRVOLUME:
457 data_[i].accumulator->add(thermo.getGyrationalVolume());
458 break;
459 case TRANSLATIONAL_KINETIC:
460 data_[i].accumulator->add(thermo.getTranslationalKinetic());
461 break;
462 case ROTATIONAL_KINETIC:
463 data_[i].accumulator->add(thermo.getRotationalKinetic());
464 break;
465 case LONG_RANGE_POTENTIAL:
466 data_[i].accumulator->add(snap->getLongRangePotential());
467 break;
468 case VANDERWAALS_POTENTIAL:
469 data_[i].accumulator->add(snap->getLongRangePotentials()[VANDERWAALS_FAMILY]);
470 break;
471 case ELECTROSTATIC_POTENTIAL:
472 data_[i].accumulator->add(snap->getLongRangePotentials()[ELECTROSTATIC_FAMILY]);
473 break;
474 case METALLIC_POTENTIAL:
475 data_[i].accumulator->add(snap->getLongRangePotentials()[METALLIC_FAMILY]);
476 break;
477 case HYDROGENBONDING_POTENTIAL:
478 data_[i].accumulator->add(snap->getLongRangePotentials()[HYDROGENBONDING_FAMILY]);
479 break;
480 case SHORT_RANGE_POTENTIAL:
481 data_[i].accumulator->add(snap->getShortRangePotential());
482 break;
483 case BOND_POTENTIAL:
484 data_[i].accumulator->add(snap->getBondPotential());
485 break;
486 case BEND_POTENTIAL:
487 data_[i].accumulator->add(snap->getBendPotential());
488 break;
489 case DIHEDRAL_POTENTIAL:
490 data_[i].accumulator->add(snap->getTorsionPotential());
491 break;
492 case INVERSION_POTENTIAL:
493 data_[i].accumulator->add(snap->getInversionPotential());
494 break;
495 case RAW_POTENTIAL:
496 data_[i].accumulator->add(snap->getRawPotential());
497 break;
498 case RESTRAINT_POTENTIAL:
499 data_[i].accumulator->add(snap->getRestraintPotential());
500 break;
501 case TAGGED_PAIR_DISTANCE:
502 data_[i].accumulator->add(thermo.getTaggedAtomPairDistance());
503 break;
504 /*
505 case RNEMD_EXCHANGE_TOTAL:
506 data_[i].accumulator->add(thermo.get_RNEMD_exchange_total());
507 break;
508 case SHADOWH:
509 data_[i].accumulator->add(thermo.getShadowHamiltionian());
510 break;
511 case HELFANDMOMENT:
512 data_[i].accumulator->add(thermo.getHelfandMoment());
513 break;
514 */
515 case ELECTRONIC_TEMPERATURE:
516 data_[i].accumulator->add(thermo.getElectronicTemperature());
517 break;
518 }
519 }
520 }
521 }
522
523 int Stats::getIntData(int index) {
524 assert(index >=0 && index < ENDINDEX);
525 RealType value;
526 data_[index].accumulator->getLastValue(value);
527 return (int) value;
528 }
529 RealType Stats::getRealData(int index) {
530 assert(index >=0 && index < ENDINDEX);
531 RealType value(0.0);
532 data_[index].accumulator->getLastValue(value);
533 return value;
534 }
535 Vector3d Stats::getVectorData(int index) {
536 assert(index >=0 && index < ENDINDEX);
537 Vector3d value;
538 dynamic_cast<VectorAccumulator*>(data_[index].accumulator)->getLastValue(value);
539 return value;
540 }
541 Mat3x3d Stats::getMatrixData(int index) {
542 assert(index >=0 && index < ENDINDEX);
543 Mat3x3d value;
544 dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)->getLastValue(value);
545 return value;
546 }
547
548 Stats::StatsBitSet Stats::getStatsMask() {
549 return statsMask_;
550 }
551 Stats::StatsMapType Stats::getStatsMap() {
552 return statsMap_;
553 }
554 void Stats::setStatsMask(Stats::StatsBitSet mask) {
555 statsMask_ = mask;
556 }
557
558 }

Properties

Name Value
svn:eol-style native
svn:executable *