ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Stats.cpp
Revision: 1808
Committed: Mon Oct 22 20:42:10 2012 UTC (12 years, 6 months ago) by gezelter
File size: 19719 byte(s)
Log Message:
A bug fix in the electric field for the new electrostatic code.  Also comment fixes for Doxygen 

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 * @version 1.0
48 */
49
50 #include "brains/Stats.hpp"
51 #include "brains/Thermo.hpp"
52
53 namespace OpenMD {
54
55 Stats::Stats(SimInfo* info) : isInit_(false), info_(info) {
56
57 if (!isInit_) {
58 init();
59 isInit_ = true;
60 }
61 }
62
63 void Stats::init() {
64
65 data_.resize(Stats::ENDINDEX);
66
67 StatsData time;
68 time.units = "fs";
69 time.title = "Time";
70 time.dataType = "RealType";
71 time.accumulator = new Accumulator();
72 data_[TIME] = time;
73 statsMap_["TIME"] = TIME;
74
75 StatsData total_energy;
76 total_energy.units = "kcal/mol";
77 total_energy.title = "Total Energy";
78 total_energy.dataType = "RealType";
79 total_energy.accumulator = new Accumulator();
80 data_[TOTAL_ENERGY] = total_energy;
81 statsMap_["TOTAL_ENERGY"] = TOTAL_ENERGY;
82
83 StatsData potential_energy;
84 potential_energy.units = "kcal/mol";
85 potential_energy.title = "Potential Energy";
86 potential_energy.dataType = "RealType";
87 potential_energy.accumulator = new Accumulator();
88 data_[POTENTIAL_ENERGY] = potential_energy;
89 statsMap_["POTENTIAL_ENERGY"] = POTENTIAL_ENERGY;
90
91 StatsData kinetic_energy;
92 kinetic_energy.units = "kcal/mol";
93 kinetic_energy.title = "Kinetic Energy";
94 kinetic_energy.dataType = "RealType";
95 kinetic_energy.accumulator = new Accumulator();
96 data_[KINETIC_ENERGY] = kinetic_energy;
97 statsMap_["KINETIC_ENERGY"] = KINETIC_ENERGY;
98
99 StatsData temperature;
100 temperature.units = "K";
101 temperature.title = "Temperature";
102 temperature.dataType = "RealType";
103 temperature.accumulator = new Accumulator();
104 data_[TEMPERATURE] = temperature;
105 statsMap_["TEMPERATURE"] = TEMPERATURE;
106
107 StatsData pressure;
108 pressure.units = "atm";
109 pressure.title = "Pressure";
110 pressure.dataType = "RealType";
111 pressure.accumulator = new Accumulator();
112 data_[PRESSURE] = pressure;
113 statsMap_["PRESSURE"] = PRESSURE;
114
115 StatsData volume;
116 volume.units = "A^3";
117 volume.title = "Volume";
118 volume.dataType = "RealType";
119 volume.accumulator = new Accumulator();
120 data_[VOLUME] = volume;
121 statsMap_["VOLUME"] = VOLUME;
122
123 StatsData hullvolume;
124 hullvolume.units = "A^3";
125 hullvolume.title = "Hull Volume";
126 hullvolume.dataType = "RealType";
127 hullvolume.accumulator = new Accumulator();
128 data_[HULLVOLUME] = hullvolume;
129 statsMap_["HULLVOLUME"] = HULLVOLUME;
130
131 StatsData gyrvolume;
132 gyrvolume.units = "A^3";
133 gyrvolume.title = "Gyrational Volume";
134 gyrvolume.dataType = "RealType";
135 gyrvolume.accumulator = new Accumulator();
136 data_[GYRVOLUME] = gyrvolume;
137 statsMap_["GYRVOLUME"] = GYRVOLUME;
138
139 StatsData conserved_quantity;
140 conserved_quantity.units = "kcal/mol";
141 conserved_quantity.title = "Conserved Quantity";
142 conserved_quantity.dataType = "RealType";
143 conserved_quantity.accumulator = new Accumulator();
144 data_[CONSERVED_QUANTITY] = conserved_quantity;
145 statsMap_["CONSERVED_QUANTITY"] = CONSERVED_QUANTITY;
146
147 StatsData translational_kinetic;
148 translational_kinetic.units = "kcal/mol";
149 translational_kinetic.title = "Translational Kinetic";
150 translational_kinetic.dataType = "RealType";
151 translational_kinetic.accumulator = new Accumulator();
152 data_[TRANSLATIONAL_KINETIC] = translational_kinetic;
153 statsMap_["TRANSLATIONAL_KINETIC"] = TRANSLATIONAL_KINETIC;
154
155 StatsData rotational_kinetic;
156 rotational_kinetic.units = "kcal/mol";
157 rotational_kinetic.title = "Rotational Kinetic";
158 rotational_kinetic.dataType = "RealType";
159 rotational_kinetic.accumulator = new Accumulator();
160 data_[ROTATIONAL_KINETIC] = rotational_kinetic;
161 statsMap_["ROTATIONAL_KINETIC"] = ROTATIONAL_KINETIC;
162
163 StatsData long_range_potential;
164 long_range_potential.units = "kcal/mol";
165 long_range_potential.title = "Long Range Potential";
166 long_range_potential.dataType = "RealType";
167 long_range_potential.accumulator = new Accumulator();
168 data_[LONG_RANGE_POTENTIAL] = long_range_potential;
169 statsMap_["LONG_RANGE_POTENTIAL"] = LONG_RANGE_POTENTIAL;
170
171 StatsData vanderwaals_potential;
172 vanderwaals_potential.units = "kcal/mol";
173 vanderwaals_potential.title = "van der waals Potential";
174 vanderwaals_potential.dataType = "RealType";
175 vanderwaals_potential.accumulator = new Accumulator();
176 data_[VANDERWAALS_POTENTIAL] = vanderwaals_potential;
177 statsMap_["VANDERWAALS_POTENTIAL"] = VANDERWAALS_POTENTIAL;
178
179 StatsData electrostatic_potential;
180 electrostatic_potential.units = "kcal/mol";
181 electrostatic_potential.title = "Electrostatic Potential";
182 electrostatic_potential.dataType = "RealType";
183 electrostatic_potential.accumulator = new Accumulator();
184 data_[ELECTROSTATIC_POTENTIAL] = electrostatic_potential;
185 statsMap_["ELECTROSTATIC_POTENTIAL"] = ELECTROSTATIC_POTENTIAL;
186
187 StatsData metallic_potential;
188 metallic_potential.units = "kcal/mol";
189 metallic_potential.title = "Metallic Potential";
190 metallic_potential.dataType = "RealType";
191 metallic_potential.accumulator = new Accumulator();
192 data_[METALLIC_POTENTIAL] = metallic_potential;
193 statsMap_["METALLIC_POTENTIAL"] = METALLIC_POTENTIAL;
194
195 StatsData hydrogenbonding_potential;
196 hydrogenbonding_potential.units = "kcal/mol";
197 hydrogenbonding_potential.title = "Metallic Potential";
198 hydrogenbonding_potential.dataType = "RealType";
199 hydrogenbonding_potential.accumulator = new Accumulator();
200 data_[HYDROGENBONDING_POTENTIAL] = hydrogenbonding_potential;
201 statsMap_["HYDROGENBONDING_POTENTIAL"] = HYDROGENBONDING_POTENTIAL;
202
203 StatsData short_range_potential;
204 short_range_potential.units = "kcal/mol";
205 short_range_potential.title = "Short Range Potential";
206 short_range_potential.dataType = "RealType";
207 short_range_potential.accumulator = new Accumulator();
208 data_[SHORT_RANGE_POTENTIAL] = short_range_potential;
209 statsMap_["SHORT_RANGE_POTENTIAL"] = SHORT_RANGE_POTENTIAL;
210
211 StatsData bond_potential;
212 bond_potential.units = "kcal/mol";
213 bond_potential.title = "Bond Potential";
214 bond_potential.dataType = "RealType";
215 bond_potential.accumulator = new Accumulator();
216 data_[BOND_POTENTIAL] = bond_potential;
217 statsMap_["BOND_POTENTIAL"] = BOND_POTENTIAL;
218
219 StatsData bend_potential;
220 bend_potential.units = "kcal/mol";
221 bend_potential.title = "Bend Potential";
222 bend_potential.dataType = "RealType";
223 bend_potential.accumulator = new Accumulator();
224 data_[BEND_POTENTIAL] = bend_potential;
225 statsMap_["BEND_POTENTIAL"] = BEND_POTENTIAL;
226
227 StatsData dihedral_potential;
228 dihedral_potential.units = "kcal/mol";
229 dihedral_potential.title = "Dihedral Potential";
230 dihedral_potential.dataType = "RealType";
231 dihedral_potential.accumulator = new Accumulator();
232 data_[DIHEDRAL_POTENTIAL] = dihedral_potential;
233 statsMap_["DIHEDRAL_POTENTIAL"] = DIHEDRAL_POTENTIAL;
234
235 StatsData inversion_potential;
236 inversion_potential.units = "kcal/mol";
237 inversion_potential.title = "Inversion Potential";
238 inversion_potential.dataType = "RealType";
239 inversion_potential.accumulator = new Accumulator();
240 data_[INVERSION_POTENTIAL] = inversion_potential;
241 statsMap_["INVERSION_POTENTIAL"] = INVERSION_POTENTIAL;
242
243 StatsData vraw;
244 vraw.units = "kcal/mol";
245 vraw.title = "Raw Potential";
246 vraw.dataType = "RealType";
247 vraw.accumulator = new Accumulator();
248 data_[RAW_POTENTIAL] = vraw;
249 statsMap_["RAW_POTENTIAL"] = RAW_POTENTIAL;
250
251 StatsData vrestraint;
252 vrestraint.units = "kcal/mol";
253 vrestraint.title = "Restraint Potential";
254 vrestraint.dataType = "RealType";
255 vrestraint.accumulator = new Accumulator();
256 data_[RESTRAINT_POTENTIAL] = vrestraint;
257 statsMap_["RESTRAINT_POTENTIAL"] = RESTRAINT_POTENTIAL;
258
259 StatsData pressure_tensor;
260 pressure_tensor.units = "amu*fs^-2*Ang^-1";
261 pressure_tensor.title = "Ptensor";
262 pressure_tensor.dataType = "Mat3x3d";
263 pressure_tensor.accumulator = new MatrixAccumulator();
264 data_[PRESSURE_TENSOR] = pressure_tensor;
265 statsMap_["PRESSURE_TENSOR"] = PRESSURE_TENSOR;
266
267 StatsData system_dipole;
268 system_dipole.units = "C*m";
269 system_dipole.title = "System Dipole";
270 system_dipole.dataType = "Vector3d";
271 system_dipole.accumulator = new VectorAccumulator();
272 data_[SYSTEM_DIPOLE] = system_dipole;
273 statsMap_["SYSTEM_DIPOLE"] = SYSTEM_DIPOLE;
274
275 StatsData tagged_pair_distance;
276 tagged_pair_distance.units = "Ang";
277 tagged_pair_distance.title = "Tagged_Pair_Distance";
278 tagged_pair_distance.dataType = "RealType";
279 tagged_pair_distance.accumulator = new Accumulator();
280 data_[TAGGED_PAIR_DISTANCE] = tagged_pair_distance;
281 statsMap_["TAGGED_PAIR_DISTANCE"] = TAGGED_PAIR_DISTANCE;
282
283 StatsData shadowh;
284 shadowh.units = "kcal/mol";
285 shadowh.title = "Shadow Hamiltonian";
286 shadowh.dataType = "RealType";
287 shadowh.accumulator = new Accumulator();
288 data_[SHADOWH] = shadowh;
289 statsMap_["SHADOWH"] = SHADOWH;
290
291 StatsData helfandmoment;
292 helfandmoment.units = "Ang*kcal/mol";
293 helfandmoment.title = "Thermal Helfand Moment";
294 helfandmoment.dataType = "Vector3d";
295 helfandmoment.accumulator = new VectorAccumulator();
296 data_[HELFANDMOMENT] = helfandmoment;
297 statsMap_["HELFANDMOMENT"] = HELFANDMOMENT;
298
299 StatsData heatflux;
300 heatflux.units = "amu/fs^3";
301 heatflux.title = "Heat Flux";
302 heatflux.dataType = "Vector3d";
303 heatflux.accumulator = new VectorAccumulator();
304 data_[HEATFLUX] = heatflux;
305 statsMap_["HEATFLUX"] = HEATFLUX;
306
307 StatsData electronic_temperature;
308 electronic_temperature.units = "K";
309 electronic_temperature.title = "Electronic Temperature";
310 electronic_temperature.dataType = "RealType";
311 electronic_temperature.accumulator = new Accumulator();
312 data_[ELECTRONIC_TEMPERATURE] = electronic_temperature;
313 statsMap_["ELECTRONIC_TEMPERATURE"] = ELECTRONIC_TEMPERATURE;
314
315 // Now, set some defaults in the mask:
316
317 Globals* simParams = info_->getSimParams();
318 std::string statFileFormatString = simParams->getStatFileFormat();
319 parseStatFileFormat(statFileFormatString);
320
321 // if we're doing a thermodynamic integration, we'll want the raw
322 // potential as well as the full potential:
323
324 if (simParams->getUseThermodynamicIntegration())
325 statsMask_.set(RAW_POTENTIAL);
326
327 // if we've got restraints turned on, we'll also want a report of the
328 // total harmonic restraints
329 if (simParams->getUseRestraints()){
330 statsMask_.set(RESTRAINT_POTENTIAL);
331 }
332
333 if (simParams->havePrintPressureTensor() &&
334 simParams->getPrintPressureTensor()){
335 statsMask_.set(PRESSURE_TENSOR);
336 }
337
338 // Why do we have both of these?
339 if (simParams->getAccumulateBoxDipole()) {
340 statsMask_.set(SYSTEM_DIPOLE);
341 }
342 if (info_->getCalcBoxDipole()){
343 statsMask_.set(SYSTEM_DIPOLE);
344 }
345
346 if (simParams->havePrintHeatFlux()) {
347 if (simParams->getPrintHeatFlux()){
348 statsMask_.set(HEATFLUX);
349 }
350 }
351
352
353 if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
354 if (simParams->getPrintTaggedPairDistance()) {
355 statsMask_.set(TAGGED_PAIR_DISTANCE);
356 }
357 }
358
359 }
360
361 void Stats::parseStatFileFormat(const std::string& format) {
362 StringTokenizer tokenizer(format, " ,;|\t\n\r");
363
364 while(tokenizer.hasMoreTokens()) {
365 std::string token(tokenizer.nextToken());
366 toUpper(token);
367 StatsMapType::iterator i = statsMap_.find(token);
368 if (i != statsMap_.end()) {
369 statsMask_.set(i->second);
370 } else {
371 sprintf( painCave.errMsg,
372 "Stats::parseStatFileFormat: %s is not a recognized\n"
373 "\tstatFileFormat keyword.\n", token.c_str() );
374 painCave.isFatal = 0;
375 painCave.severity = OPENMD_ERROR;
376 simError();
377 }
378 }
379 }
380
381
382 std::string Stats::getTitle(int index) {
383 assert(index >=0 && index < ENDINDEX);
384 return data_[index].title;
385 }
386
387 std::string Stats::getUnits(int index) {
388 assert(index >=0 && index < ENDINDEX);
389 return data_[index].units;
390 }
391
392 std::string Stats::getDataType(int index) {
393 assert(index >=0 && index < ENDINDEX);
394 return data_[index].dataType;
395 }
396
397 void Stats::collectStats(){
398 Globals* simParams = info_->getSimParams();
399 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
400 Thermo thermo(info_);
401
402 for (unsigned int i = 0; i < statsMask_.size(); ++i) {
403 if (statsMask_[i]) {
404 switch (i) {
405 case TIME:
406 data_[i].accumulator->add(snap->getTime());
407 break;
408 case KINETIC_ENERGY:
409 data_[i].accumulator->add(thermo.getKinetic());
410 break;
411 case POTENTIAL_ENERGY:
412 data_[i].accumulator->add(thermo.getPotential());
413 break;
414 case TOTAL_ENERGY:
415 data_[i].accumulator->add(thermo.getTotalEnergy());
416 break;
417 case TEMPERATURE:
418 data_[i].accumulator->add(thermo.getTemperature());
419 break;
420 case PRESSURE:
421 data_[i].accumulator->add(thermo.getPressure());
422 break;
423 case VOLUME:
424 data_[i].accumulator->add(thermo.getVolume());
425 break;
426 case CONSERVED_QUANTITY:
427 data_[i].accumulator->add(snap->getConservedQuantity());
428 break;
429 case PRESSURE_TENSOR:
430 dynamic_cast<MatrixAccumulator *>(data_[i].accumulator)->add(thermo.getPressureTensor());
431 break;
432 case SYSTEM_DIPOLE:
433 dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getSystemDipole());
434 break;
435 case HEATFLUX:
436 dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getHeatFlux());
437 break;
438 case HULLVOLUME:
439 data_[i].accumulator->add(thermo.getHullVolume());
440 break;
441 case GYRVOLUME:
442 data_[i].accumulator->add(thermo.getGyrationalVolume());
443 break;
444 case TRANSLATIONAL_KINETIC:
445 data_[i].accumulator->add(thermo.getTranslationalKinetic());
446 break;
447 case ROTATIONAL_KINETIC:
448 data_[i].accumulator->add(thermo.getRotationalKinetic());
449 break;
450 case LONG_RANGE_POTENTIAL:
451 data_[i].accumulator->add(snap->getLongRangePotential());
452 break;
453 case VANDERWAALS_POTENTIAL:
454 data_[i].accumulator->add(snap->getLongRangePotentials()[VANDERWAALS_FAMILY]);
455 break;
456 case ELECTROSTATIC_POTENTIAL:
457 data_[i].accumulator->add(snap->getLongRangePotentials()[ELECTROSTATIC_FAMILY]);
458 break;
459 case METALLIC_POTENTIAL:
460 data_[i].accumulator->add(snap->getLongRangePotentials()[METALLIC_FAMILY]);
461 break;
462 case HYDROGENBONDING_POTENTIAL:
463 data_[i].accumulator->add(snap->getLongRangePotentials()[HYDROGENBONDING_FAMILY]);
464 break;
465 case SHORT_RANGE_POTENTIAL:
466 data_[i].accumulator->add(snap->getShortRangePotential());
467 break;
468 case BOND_POTENTIAL:
469 data_[i].accumulator->add(snap->getBondPotential());
470 break;
471 case BEND_POTENTIAL:
472 data_[i].accumulator->add(snap->getBendPotential());
473 break;
474 case DIHEDRAL_POTENTIAL:
475 data_[i].accumulator->add(snap->getTorsionPotential());
476 break;
477 case INVERSION_POTENTIAL:
478 data_[i].accumulator->add(snap->getInversionPotential());
479 break;
480 case RAW_POTENTIAL:
481 data_[i].accumulator->add(snap->getRawPotential());
482 break;
483 case RESTRAINT_POTENTIAL:
484 data_[i].accumulator->add(snap->getRestraintPotential());
485 break;
486 case TAGGED_PAIR_DISTANCE:
487 data_[i].accumulator->add(thermo.getTaggedAtomPairDistance());
488 break;
489 /*
490 case SHADOWH:
491 data_[i].accumulator->add(thermo.getShadowHamiltionian());
492 break;
493 case HELFANDMOMENT:
494 data_[i].accumulator->add(thermo.getHelfandMoment());
495 break;
496 */
497 case ELECTRONIC_TEMPERATURE:
498 data_[i].accumulator->add(thermo.getElectronicTemperature());
499 break;
500 }
501 }
502 }
503 }
504
505 int Stats::getIntData(int index) {
506 assert(index >=0 && index < ENDINDEX);
507 RealType value;
508 data_[index].accumulator->getLastValue(value);
509 return (int) value;
510 }
511 RealType Stats::getRealData(int index) {
512 assert(index >=0 && index < ENDINDEX);
513 RealType value(0.0);
514 data_[index].accumulator->getLastValue(value);
515 return value;
516 }
517 Vector3d Stats::getVectorData(int index) {
518 assert(index >=0 && index < ENDINDEX);
519 Vector3d value;
520 dynamic_cast<VectorAccumulator*>(data_[index].accumulator)->getLastValue(value);
521 return value;
522 }
523 Mat3x3d Stats::getMatrixData(int index) {
524 assert(index >=0 && index < ENDINDEX);
525 Mat3x3d value;
526 dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)->getLastValue(value);
527 return value;
528 }
529
530 Stats::StatsBitSet Stats::getStatsMask() {
531 return statsMask_;
532 }
533 Stats::StatsMapType Stats::getStatsMap() {
534 return statsMap_;
535 }
536 void Stats::setStatsMask(Stats::StatsBitSet mask) {
537 statsMask_ = mask;
538 }
539
540 }

Properties

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