src/brains/Stats.cpp

/*
 * Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
  
/**
 * @file Stats.cpp
 * @author tlin
 * @date 11/04/2004
 * @time 14:26am
 * @version 1.0
 */

#include "brains/Stats.hpp"
#include "brains/Thermo.hpp"

namespace OpenMD {

  Stats::Stats(SimInfo* info) : isInit_(false), info_(info) {   

    if (!isInit_) {
      init();
      isInit_ = true;
    }
  }

  void Stats::init() {
   
    data_.resize(Stats::ENDINDEX);

    StatsData time;
    time.units =  "fs";
    time.title =  "Time";
    time.dataType = "RealType";
    time.accumulator = new Accumulator();
    data_[TIME] = time;
    statsMap_["TIME"] = TIME;

    StatsData total_energy;
    total_energy.units =  "kcal/mol";
    total_energy.title =  "Total Energy";
    total_energy.dataType = "RealType";
    total_energy.accumulator = new Accumulator();
    data_[TOTAL_ENERGY] = total_energy;
    statsMap_["TOTAL_ENERGY"] =  TOTAL_ENERGY;
   
    StatsData potential_energy;
    potential_energy.units =  "kcal/mol";
    potential_energy.title =  "Potential Energy";
    potential_energy.dataType = "RealType";
    potential_energy.accumulator = new Accumulator();
    data_[POTENTIAL_ENERGY] = potential_energy;
    statsMap_["POTENTIAL_ENERGY"] =  POTENTIAL_ENERGY;

    StatsData kinetic_energy;
    kinetic_energy.units =  "kcal/mol";
    kinetic_energy.title =  "Kinetic Energy";
    kinetic_energy.dataType = "RealType";
    kinetic_energy.accumulator = new Accumulator();
    data_[KINETIC_ENERGY] = kinetic_energy;
    statsMap_["KINETIC_ENERGY"] =  KINETIC_ENERGY;

    StatsData temperature;
    temperature.units =  "K";
    temperature.title =  "Temperature";
    temperature.dataType = "RealType";
    temperature.accumulator = new Accumulator();
    data_[TEMPERATURE] = temperature;
    statsMap_["TEMPERATURE"] =  TEMPERATURE;

    StatsData pressure;
    pressure.units =  "atm";
    pressure.title =  "Pressure";
    pressure.dataType = "RealType";
    pressure.accumulator = new Accumulator();
    data_[PRESSURE] = pressure;
    statsMap_["PRESSURE"] =  PRESSURE;

    StatsData volume;
    volume.units =  "A^3";
    volume.title =  "Volume";
    volume.dataType = "RealType";
    volume.accumulator = new Accumulator();
    data_[VOLUME] = volume;
    statsMap_["VOLUME"] =  VOLUME;

    StatsData hullvolume;
    hullvolume.units =  "A^3";
    hullvolume.title =  "Hull Volume";
    hullvolume.dataType = "RealType";
    hullvolume.accumulator = new Accumulator();
    data_[HULLVOLUME] = hullvolume;
    statsMap_["HULLVOLUME"] =  HULLVOLUME;

    StatsData gyrvolume;
    gyrvolume.units =  "A^3";
    gyrvolume.title =  "Gyrational Volume";
    gyrvolume.dataType = "RealType";
    gyrvolume.accumulator = new Accumulator();
    data_[GYRVOLUME] = gyrvolume;
    statsMap_["GYRVOLUME"] =  GYRVOLUME;

    StatsData conserved_quantity;
    conserved_quantity.units =  "kcal/mol";             
    conserved_quantity.title =  "Conserved Quantity";             
    conserved_quantity.dataType = "RealType";
    conserved_quantity.accumulator = new Accumulator();
    data_[CONSERVED_QUANTITY] = conserved_quantity;
    statsMap_["CONSERVED_QUANTITY"] =  CONSERVED_QUANTITY;

    StatsData translational_kinetic;
    translational_kinetic.units =  "kcal/mol";
    translational_kinetic.title =  "Translational Kinetic";
    translational_kinetic.dataType = "RealType";
    translational_kinetic.accumulator = new Accumulator();
    data_[TRANSLATIONAL_KINETIC] = translational_kinetic;
    statsMap_["TRANSLATIONAL_KINETIC"] =  TRANSLATIONAL_KINETIC;

    StatsData rotational_kinetic;
    rotational_kinetic.units =  "kcal/mol";
    rotational_kinetic.title =  "Rotational Kinetic";
    rotational_kinetic.dataType = "RealType";
    rotational_kinetic.accumulator = new Accumulator();
    data_[ROTATIONAL_KINETIC] = rotational_kinetic;
    statsMap_["ROTATIONAL_KINETIC"] =  ROTATIONAL_KINETIC;

    StatsData long_range_potential;
    long_range_potential.units =  "kcal/mol";
    long_range_potential.title =  "Long Range Potential";
    long_range_potential.dataType = "RealType";
    long_range_potential.accumulator = new Accumulator();
    data_[LONG_RANGE_POTENTIAL] = long_range_potential;
    statsMap_["LONG_RANGE_POTENTIAL"] =  LONG_RANGE_POTENTIAL;

    StatsData vanderwaals_potential;
    vanderwaals_potential.units =  "kcal/mol";
    vanderwaals_potential.title =  "van der waals Potential";
    vanderwaals_potential.dataType = "RealType";
    vanderwaals_potential.accumulator = new Accumulator();
    data_[VANDERWAALS_POTENTIAL] = vanderwaals_potential;
    statsMap_["VANDERWAALS_POTENTIAL"] =  VANDERWAALS_POTENTIAL;

    StatsData electrostatic_potential;
    electrostatic_potential.units =  "kcal/mol";
    electrostatic_potential.title =  "Electrostatic Potential";    
    electrostatic_potential.dataType = "RealType";
    electrostatic_potential.accumulator = new Accumulator();
    data_[ELECTROSTATIC_POTENTIAL] = electrostatic_potential;
    statsMap_["ELECTROSTATIC_POTENTIAL"] =  ELECTROSTATIC_POTENTIAL;

    StatsData metallic_potential;
    metallic_potential.units =  "kcal/mol";
    metallic_potential.title =  "Metallic Potential";    
    metallic_potential.dataType = "RealType";
    metallic_potential.accumulator = new Accumulator();
    data_[METALLIC_POTENTIAL] = metallic_potential;
    statsMap_["METALLIC_POTENTIAL"] =  METALLIC_POTENTIAL;

    StatsData hydrogenbonding_potential;
    hydrogenbonding_potential.units =  "kcal/mol";
    hydrogenbonding_potential.title =  "Metallic Potential";    
    hydrogenbonding_potential.dataType = "RealType";
    hydrogenbonding_potential.accumulator = new Accumulator();
    data_[HYDROGENBONDING_POTENTIAL] = hydrogenbonding_potential;
    statsMap_["HYDROGENBONDING_POTENTIAL"] =  HYDROGENBONDING_POTENTIAL;

    StatsData short_range_potential;
    short_range_potential.units =  "kcal/mol";
    short_range_potential.title =  "Short Range Potential";
    short_range_potential.dataType = "RealType";
    short_range_potential.accumulator = new Accumulator();
    data_[SHORT_RANGE_POTENTIAL] = short_range_potential;
    statsMap_["SHORT_RANGE_POTENTIAL"] =  SHORT_RANGE_POTENTIAL;

    StatsData bond_potential;
    bond_potential.units =  "kcal/mol";
    bond_potential.title =  "Bond Potential";
    bond_potential.dataType = "RealType";
    bond_potential.accumulator = new Accumulator();
    data_[BOND_POTENTIAL] = bond_potential;
    statsMap_["BOND_POTENTIAL"] =  BOND_POTENTIAL;

    StatsData bend_potential;
    bend_potential.units =  "kcal/mol";
    bend_potential.title =  "Bend Potential";
    bend_potential.dataType = "RealType";
    bend_potential.accumulator = new Accumulator();
    data_[BEND_POTENTIAL] = bend_potential;
    statsMap_["BEND_POTENTIAL"] =  BEND_POTENTIAL;
    
    StatsData dihedral_potential;
    dihedral_potential.units =  "kcal/mol";
    dihedral_potential.title =  "Dihedral Potential";
    dihedral_potential.dataType = "RealType";
    dihedral_potential.accumulator = new Accumulator();
    data_[DIHEDRAL_POTENTIAL] = dihedral_potential;
    statsMap_["DIHEDRAL_POTENTIAL"] =  DIHEDRAL_POTENTIAL;

    StatsData inversion_potential;
    inversion_potential.units =  "kcal/mol";
    inversion_potential.title =  "Inversion Potential";
    inversion_potential.dataType = "RealType";
    inversion_potential.accumulator = new Accumulator();
    data_[INVERSION_POTENTIAL] = inversion_potential;
    statsMap_["INVERSION_POTENTIAL"] =  INVERSION_POTENTIAL;

    StatsData vraw;
    vraw.units =  "kcal/mol";
    vraw.title =  "Raw Potential";
    vraw.dataType = "RealType";
    vraw.accumulator = new Accumulator();
    data_[RAW_POTENTIAL] = vraw;
    statsMap_["RAW_POTENTIAL"] =  RAW_POTENTIAL;

    StatsData vrestraint;
    vrestraint.units =  "kcal/mol";
    vrestraint.title =  "Restraint Potential";
    vrestraint.dataType = "RealType";
    vrestraint.accumulator = new Accumulator();
    data_[RESTRAINT_POTENTIAL] = vrestraint;
    statsMap_["RESTRAINT_POTENTIAL"] =  RESTRAINT_POTENTIAL;

    StatsData pressure_tensor;
    pressure_tensor.units =  "amu*fs^-2*Ang^-1";
    pressure_tensor.title =  "Ptensor";
    pressure_tensor.dataType = "Mat3x3d";
    pressure_tensor.accumulator = new MatrixAccumulator();
    data_[PRESSURE_TENSOR] = pressure_tensor;
    statsMap_["PRESSURE_TENSOR"] =  PRESSURE_TENSOR;

    StatsData system_dipole;
    system_dipole.units =  "C*m";
    system_dipole.title =  "System Dipole";
    system_dipole.dataType = "Vector3d";
    system_dipole.accumulator = new VectorAccumulator();
    data_[SYSTEM_DIPOLE] = system_dipole;
    statsMap_["SYSTEM_DIPOLE"] =  SYSTEM_DIPOLE;

    StatsData tagged_pair_distance;
    tagged_pair_distance.units =  "Ang";
    tagged_pair_distance.title =  "Tagged_Pair_Distance";
    tagged_pair_distance.dataType = "RealType";
    tagged_pair_distance.accumulator = new Accumulator();
    data_[TAGGED_PAIR_DISTANCE] = tagged_pair_distance;
    statsMap_["TAGGED_PAIR_DISTANCE"] =  TAGGED_PAIR_DISTANCE;

    StatsData rnemd_exchange_total;
    rnemd_exchange_total.units =  "Variable";
    rnemd_exchange_total.title =  "RNEMD_exchange_total";
    rnemd_exchange_total.dataType = "RealType";
    rnemd_exchange_total.accumulator = new Accumulator();
    data_[RNEMD_EXCHANGE_TOTAL] = rnemd_exchange_total;
    statsMap_["RNEMD_EXCHANGE_TOTAL"] =  RNEMD_EXCHANGE_TOTAL;

    StatsData shadowh;
    shadowh.units =  "kcal/mol";
    shadowh.title =  "Shadow Hamiltonian";
    shadowh.dataType = "RealType";
    shadowh.accumulator = new Accumulator();
    data_[SHADOWH] = shadowh;
    statsMap_["SHADOWH"] =  SHADOWH;

    StatsData helfandmoment;
    helfandmoment.units =  "Ang*kcal/mol";
    helfandmoment.title =  "Thermal Helfand Moment";
    helfandmoment.dataType = "Vector3d";
    helfandmoment.accumulator = new VectorAccumulator();
    data_[HELFANDMOMENT] = helfandmoment;
    statsMap_["HELFANDMOMENT"] = HELFANDMOMENT;

    StatsData heatflux;
    heatflux.units = "amu/fs^3";
    heatflux.title =  "Heat Flux";  
    heatflux.dataType = "Vector3d";
    heatflux.accumulator = new VectorAccumulator();
    data_[HEATFLUX] = heatflux;
    statsMap_["HEATFLUX"] = HEATFLUX;

    StatsData electronic_temperature;
    electronic_temperature.units = "K";
    electronic_temperature.title =  "Electronic Temperature";  
    electronic_temperature.dataType = "RealType";
    electronic_temperature.accumulator = new Accumulator();
    data_[ELECTRONIC_TEMPERATURE] = electronic_temperature;
    statsMap_["ELECTRONIC_TEMPERATURE"] = ELECTRONIC_TEMPERATURE;

    // Now, set some defaults in the mask:

    Globals* simParams = info_->getSimParams();
    std::string statFileFormatString = simParams->getStatFileFormat();
    parseStatFileFormat(statFileFormatString);

    // if we're doing a thermodynamic integration, we'll want the raw
    // potential as well as the full potential:
    
    if (simParams->getUseThermodynamicIntegration()) 
      statsMask_.set(RAW_POTENTIAL);
    
    // if we've got restraints turned on, we'll also want a report of the
    // total harmonic restraints
    if (simParams->getUseRestraints()){
      statsMask_.set(RESTRAINT_POTENTIAL);
    }
    
    if (simParams->havePrintPressureTensor() && 
        simParams->getPrintPressureTensor()){
      statsMask_.set(PRESSURE_TENSOR);
    }

    // Why do we have both of these?
    if (simParams->getAccumulateBoxDipole()) {
      statsMask_.set(SYSTEM_DIPOLE);
    }
    if (info_->getCalcBoxDipole()){
      statsMask_.set(SYSTEM_DIPOLE);
    }

    if (simParams->havePrintHeatFlux()) {
      if (simParams->getPrintHeatFlux()){
        statsMask_.set(HEATFLUX);
      }
    }    
    
    
    if (simParams->haveTaggedAtomPair() && simParams->havePrintTaggedPairDistance()) {
      if (simParams->getPrintTaggedPairDistance()) {
        statsMask_.set(TAGGED_PAIR_DISTANCE);
      }
    }
    
    if (simParams->getRNEMDParameters()->getUseRNEMD()) {
      statsMask_.set(RNEMD_EXCHANGE_TOTAL);
    }


  }

  void Stats::parseStatFileFormat(const std::string& format) {
    StringTokenizer tokenizer(format, " ,;|\t\n\r");

    while(tokenizer.hasMoreTokens()) {
      std::string token(tokenizer.nextToken());
      toUpper(token);
      StatsMapType::iterator i = statsMap_.find(token);
      if (i != statsMap_.end()) {
        statsMask_.set(i->second);
      } else {
        sprintf( painCave.errMsg,
                 "Stats::parseStatFileFormat: %s is not a recognized\n"
                 "\tstatFileFormat keyword.\n", token.c_str() );
        painCave.isFatal = 0;
        painCave.severity = OPENMD_ERROR;
        simError();            
      }
    }   
  }


  std::string Stats::getTitle(int index) {
    assert(index >=0 && index < ENDINDEX);
    return data_[index].title;
  }

  std::string Stats::getUnits(int index) {
    assert(index >=0 && index < ENDINDEX);
    return data_[index].units;
  }

  std::string Stats::getDataType(int index) {
    assert(index >=0 && index < ENDINDEX);
    return data_[index].dataType;
  }

  void Stats::collectStats(){
    Globals* simParams = info_->getSimParams();
    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
    Thermo thermo(info_);
   
    for (int i = 0; i < statsMask_.size(); ++i) {
      if (statsMask_[i]) {
        switch (i) {
        case TIME:
          data_[i].accumulator->add(snap->getTime());
          break;
        case KINETIC_ENERGY:
          data_[i].accumulator->add(thermo.getKinetic());
          break;
        case POTENTIAL_ENERGY:
          data_[i].accumulator->add(thermo.getPotential());
          break;
        case TOTAL_ENERGY:
          data_[i].accumulator->add(thermo.getTotalEnergy());
          break;
        case TEMPERATURE:
          data_[i].accumulator->add(thermo.getTemperature());
          break;
        case PRESSURE:
          data_[i].accumulator->add(thermo.getPressure());
          break;
        case VOLUME:
          data_[i].accumulator->add(thermo.getVolume());
          break;
        case CONSERVED_QUANTITY:
          data_[i].accumulator->add(snap->getConservedQuantity());
          break;
        case PRESSURE_TENSOR:
          dynamic_cast<MatrixAccumulator *>(data_[i].accumulator)->add(thermo.getPressureTensor());
          break;
        case SYSTEM_DIPOLE:
          dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getSystemDipole());
          break;
        case HEATFLUX:
          dynamic_cast<VectorAccumulator *>(data_[i].accumulator)->add(thermo.getHeatFlux());
          break;
        case HULLVOLUME:
          data_[i].accumulator->add(thermo.getHullVolume());
          break;
        case GYRVOLUME:
          data_[i].accumulator->add(thermo.getGyrationalVolume());
          break;
        case TRANSLATIONAL_KINETIC:
          data_[i].accumulator->add(thermo.getTranslationalKinetic());
          break;
        case ROTATIONAL_KINETIC:
          data_[i].accumulator->add(thermo.getRotationalKinetic());
          break;
        case LONG_RANGE_POTENTIAL:
          data_[i].accumulator->add(snap->getLongRangePotential());
          break;
        case VANDERWAALS_POTENTIAL:
          data_[i].accumulator->add(snap->getLongRangePotentials()[VANDERWAALS_FAMILY]);
          break;
        case ELECTROSTATIC_POTENTIAL:
          data_[i].accumulator->add(snap->getLongRangePotentials()[ELECTROSTATIC_FAMILY]);
          break;
        case METALLIC_POTENTIAL:
          data_[i].accumulator->add(snap->getLongRangePotentials()[METALLIC_FAMILY]);
          break;
        case HYDROGENBONDING_POTENTIAL:
          data_[i].accumulator->add(snap->getLongRangePotentials()[HYDROGENBONDING_FAMILY]);
          break;
        case SHORT_RANGE_POTENTIAL:
          data_[i].accumulator->add(snap->getShortRangePotential());
          break;
        case BOND_POTENTIAL:
          data_[i].accumulator->add(snap->getBondPotential());
          break;
        case BEND_POTENTIAL:
          data_[i].accumulator->add(snap->getBendPotential());
          break;
        case DIHEDRAL_POTENTIAL:
          data_[i].accumulator->add(snap->getTorsionPotential());
          break;
        case INVERSION_POTENTIAL:
          data_[i].accumulator->add(snap->getInversionPotential());
          break;
        case RAW_POTENTIAL:
          data_[i].accumulator->add(snap->getRawPotential());
          break;
        case RESTRAINT_POTENTIAL:
          data_[i].accumulator->add(snap->getRestraintPotential());
          break;
        case TAGGED_PAIR_DISTANCE:
          data_[i].accumulator->add(thermo.getTaggedAtomPairDistance());
          break;
          /*
        case RNEMD_EXCHANGE_TOTAL:
          data_[i].accumulator->add(thermo.get_RNEMD_exchange_total());
          break;
        case SHADOWH:
          data_[i].accumulator->add(thermo.getShadowHamiltionian());
          break;
        case HELFANDMOMENT:
          data_[i].accumulator->add(thermo.getHelfandMoment());
          break;
          */
        case ELECTRONIC_TEMPERATURE:
          data_[i].accumulator->add(thermo.getElectronicTemperature());
          break; 
        }
      }
    }
  }

  int Stats::getIntData(int index) { 
    assert(index >=0 && index < ENDINDEX);
    RealType value;
    data_[index].accumulator->getLastValue(value);
    return (int) value;
  }
  RealType Stats::getRealData(int index) {
    assert(index >=0 && index < ENDINDEX);
    RealType value(0.0);
    data_[index].accumulator->getLastValue(value);
    return value;
  }
  Vector3d Stats::getVectorData(int index) {
    assert(index >=0 && index < ENDINDEX);
    Vector3d value;
    dynamic_cast<VectorAccumulator*>(data_[index].accumulator)->getLastValue(value);
    return value;
  }
  Mat3x3d Stats::getMatrixData(int index) {
    assert(index >=0 && index < ENDINDEX);
    Mat3x3d value;
    dynamic_cast<MatrixAccumulator*>(data_[index].accumulator)->getLastValue(value);
    return value;
  }
 
  Stats::StatsBitSet Stats::getStatsMask() {
    return statsMask_;
  }
  Stats::StatsMapType Stats::getStatsMap() {
    return statsMap_;
  }
  void Stats::setStatsMask(Stats::StatsBitSet mask) {
    statsMask_ = mask;
  }

}
Revision:	1764
Committed:	Tue Jul 3 18:32:27 2012 UTC (12 years, 10 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.
#	User	Rev	Content
1	gezelter	1710	/*
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	gezelter	1764	#include "brains/Thermo.hpp"
53	gezelter	1710
54			namespace OpenMD {
55
56	gezelter	1764	Stats::Stats(SimInfo* info) : isInit_(false), info_(info) {
57	gezelter	1710
58			if (!isInit_) {
59			init();
60			isInit_ = true;
61			}
62			}
63
64			void Stats::init() {
65	gezelter	1764
66			data_.resize(Stats::ENDINDEX);
67	gezelter	1710
68	gezelter	1764	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	gezelter	1723
76	gezelter	1764	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	gezelter	1710
92	gezelter	1764	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 = "amufs^-2Ang^-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	gezelter	1710	}
375
376	gezelter	1764	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	gezelter	1710	}