ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Stats.cpp
Revision: 1833
Committed: Tue Jan 15 16:28:22 2013 UTC (12 years, 3 months ago) by gezelter
File size: 20614 byte(s)
Log Message:
Bug fix in Sticky potential, column label in stats

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 = "Hydrogen Bonding 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 shadowh;
285 shadowh.units = "kcal/mol";
286 shadowh.title = "Shadow Hamiltonian";
287 shadowh.dataType = "RealType";
288 shadowh.accumulator = new Accumulator();
289 data_[SHADOWH] = shadowh;
290 statsMap_["SHADOWH"] = SHADOWH;
291
292 StatsData helfandmoment;
293 helfandmoment.units = "Ang*kcal/mol";
294 helfandmoment.title = "Thermal Helfand Moment";
295 helfandmoment.dataType = "Vector3d";
296 helfandmoment.accumulator = new VectorAccumulator();
297 data_[HELFANDMOMENT] = helfandmoment;
298 statsMap_["HELFANDMOMENT"] = HELFANDMOMENT;
299
300 StatsData heatflux;
301 heatflux.units = "amu/fs^3";
302 heatflux.title = "Heat Flux";
303 heatflux.dataType = "Vector3d";
304 heatflux.accumulator = new VectorAccumulator();
305 data_[HEATFLUX] = heatflux;
306 statsMap_["HEATFLUX"] = HEATFLUX;
307
308 StatsData electronic_temperature;
309 electronic_temperature.units = "K";
310 electronic_temperature.title = "Electronic Temperature";
311 electronic_temperature.dataType = "RealType";
312 electronic_temperature.accumulator = new Accumulator();
313 data_[ELECTRONIC_TEMPERATURE] = electronic_temperature;
314 statsMap_["ELECTRONIC_TEMPERATURE"] = ELECTRONIC_TEMPERATURE;
315
316 // Now, set some defaults in the mask:
317
318 Globals* simParams = info_->getSimParams();
319 std::string statFileFormatString = simParams->getStatFileFormat();
320 parseStatFileFormat(statFileFormatString);
321
322 // if we're doing a thermodynamic integration, we'll want the raw
323 // potential as well as the full potential:
324
325 if (simParams->getUseThermodynamicIntegration())
326 statsMask_.set(RAW_POTENTIAL);
327
328 // if we've got restraints turned on, we'll also want a report of the
329 // total harmonic restraints
330 if (simParams->getUseRestraints()){
331 statsMask_.set(RESTRAINT_POTENTIAL);
332 }
333
334 if (simParams->havePrintPressureTensor() &&
335 simParams->getPrintPressureTensor()){
336 statsMask_.set(PRESSURE_TENSOR);
337 }
338
339 // Why do we have both of these?
340 if (simParams->getAccumulateBoxDipole()) {
341 statsMask_.set(SYSTEM_DIPOLE);
342 }
343 if (info_->getCalcBoxDipole()){
344 statsMask_.set(SYSTEM_DIPOLE);
345 }
346
347 if (simParams->havePrintHeatFlux()) {
348 if (simParams->getPrintHeatFlux()){
349 statsMask_.set(HEATFLUX);
350 }
351 }
352
353
354 if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
355 if (simParams->getPrintTaggedPairDistance()) {
356 statsMask_.set(TAGGED_PAIR_DISTANCE);
357 }
358 }
359
360 }
361
362 void Stats::parseStatFileFormat(const std::string& format) {
363 StringTokenizer tokenizer(format, " ,;|\t\n\r");
364
365 while(tokenizer.hasMoreTokens()) {
366 std::string token(tokenizer.nextToken());
367 toUpper(token);
368 StatsMapType::iterator i = statsMap_.find(token);
369 if (i != statsMap_.end()) {
370 statsMask_.set(i->second);
371 } else {
372 sprintf( painCave.errMsg,
373 "Stats::parseStatFileFormat: %s is not a recognized\n"
374 "\tstatFileFormat keyword.\n", token.c_str() );
375 painCave.isFatal = 0;
376 painCave.severity = OPENMD_ERROR;
377 simError();
378 }
379 }
380 }
381
382
383 std::string Stats::getTitle(int index) {
384 assert(index >=0 && index < ENDINDEX);
385 return data_[index].title;
386 }
387
388 std::string Stats::getUnits(int index) {
389 assert(index >=0 && index < ENDINDEX);
390 return data_[index].units;
391 }
392
393 std::string Stats::getDataType(int index) {
394 assert(index >=0 && index < ENDINDEX);
395 return data_[index].dataType;
396 }
397
398 void Stats::collectStats(){
399 Globals* simParams = info_->getSimParams();
400 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
401 Thermo thermo(info_);
402
403 for (unsigned int i = 0; i < statsMask_.size(); ++i) {
404 if (statsMask_[i]) {
405 switch (i) {
406 case TIME:
407 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getTime());
408 break;
409 case KINETIC_ENERGY:
410 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getKinetic());
411 break;
412 case POTENTIAL_ENERGY:
413 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getPotential());
414 break;
415 case TOTAL_ENERGY:
416 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getTotalEnergy());
417 break;
418 case TEMPERATURE:
419 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getTemperature());
420 break;
421 case PRESSURE:
422 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getPressure());
423 break;
424 case VOLUME:
425 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getVolume());
426 break;
427 case CONSERVED_QUANTITY:
428 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getConservedQuantity());
429 break;
430 case PRESSURE_TENSOR:
431 dynamic_cast<MatrixAccumulator *>(data_[i].accumulator)->add(thermo.getPressureTensor());
432 break;
433 case SYSTEM_DIPOLE:
434 dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getSystemDipole());
435 break;
436 case HEATFLUX:
437 dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getHeatFlux());
438 break;
439 case HULLVOLUME:
440 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getHullVolume());
441 break;
442 case GYRVOLUME:
443 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getGyrationalVolume());
444 break;
445 case TRANSLATIONAL_KINETIC:
446 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getTranslationalKinetic());
447 break;
448 case ROTATIONAL_KINETIC:
449 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getRotationalKinetic());
450 break;
451 case LONG_RANGE_POTENTIAL:
452 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getLongRangePotential());
453 break;
454 case VANDERWAALS_POTENTIAL:
455 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getLongRangePotentials()[VANDERWAALS_FAMILY]);
456 break;
457 case ELECTROSTATIC_POTENTIAL:
458 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getLongRangePotentials()[ELECTROSTATIC_FAMILY]);
459 break;
460 case METALLIC_POTENTIAL:
461 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getLongRangePotentials()[METALLIC_FAMILY]);
462 break;
463 case HYDROGENBONDING_POTENTIAL:
464 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getLongRangePotentials()[HYDROGENBONDING_FAMILY]);
465 break;
466 case SHORT_RANGE_POTENTIAL:
467 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getShortRangePotential());
468 break;
469 case BOND_POTENTIAL:
470 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getBondPotential());
471 break;
472 case BEND_POTENTIAL:
473 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getBendPotential());
474 break;
475 case DIHEDRAL_POTENTIAL:
476 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getTorsionPotential());
477 break;
478 case INVERSION_POTENTIAL:
479 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getInversionPotential());
480 break;
481 case RAW_POTENTIAL:
482 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getRawPotential());
483 break;
484 case RESTRAINT_POTENTIAL:
485 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(snap->getRestraintPotential());
486 break;
487 case TAGGED_PAIR_DISTANCE:
488 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getTaggedAtomPairDistance());
489 break;
490 /*
491 case SHADOWH:
492 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getShadowHamiltionian());
493 break;
494 case HELFANDMOMENT:
495 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getHelfandMoment());
496 break;
497 */
498 case ELECTRONIC_TEMPERATURE:
499 dynamic_cast<Accumulator *>(data_[i].accumulator)->add(thermo.getElectronicTemperature());
500 break;
501 }
502 }
503 }
504 }
505
506 int Stats::getIntData(int index) {
507 assert(index >=0 && index < ENDINDEX);
508 RealType value;
509 dynamic_cast<Accumulator *>(data_[index].accumulator)->getLastValue(value);
510 return (int) value;
511 }
512 RealType Stats::getRealData(int index) {
513 assert(index >=0 && index < ENDINDEX);
514 RealType value(0.0);
515 dynamic_cast<Accumulator *>(data_[index].accumulator)->getLastValue(value);
516 return value;
517 }
518 Vector3d Stats::getVectorData(int index) {
519 assert(index >=0 && index < ENDINDEX);
520 Vector3d value;
521 dynamic_cast<VectorAccumulator*>(data_[index].accumulator)->getLastValue(value);
522 return value;
523 }
524 Mat3x3d Stats::getMatrixData(int index) {
525 assert(index >=0 && index < ENDINDEX);
526 Mat3x3d value;
527 dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)->getLastValue(value);
528 return value;
529 }
530
531 Stats::StatsBitSet Stats::getStatsMask() {
532 return statsMask_;
533 }
534 Stats::StatsMapType Stats::getStatsMap() {
535 return statsMap_;
536 }
537 void Stats::setStatsMask(Stats::StatsBitSet mask) {
538 statsMask_ = mask;
539 }
540
541 }

Properties

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