src/brains/Snapshot.cpp

/*
 * Copyright (c) 2005 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 Snapshot.cpp
 * @author tlin
 * @date 11/11/2004
 * @time 10:56am
 * @version 1.0
 */

#include "brains/Snapshot.hpp"
#include "utils/NumericConstant.hpp"
#include "utils/simError.h"
#include "utils/Utility.hpp"
#include <cstdio>

namespace OpenMD {

  Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups) : 
    atomData(nAtoms), rigidbodyData(nRigidbodies),
    cgData(nCutoffGroups, DataStorage::dslPosition), 
    orthoTolerance_(1e-6) {
    
    frameData.id = -1;                   
    frameData.currentTime = 0;     
    frameData.hmat = Mat3x3d(0.0);             
    frameData.invHmat = Mat3x3d(0.0);          
    frameData.orthoRhombic = false;        
    frameData.bondPotential = 0.0;      
    frameData.bendPotential = 0.0;      
    frameData.torsionPotential = 0.0;   
    frameData.inversionPotential = 0.0; 
    frameData.lrPotentials = potVec(0.0);
    frameData.excludedPotentials = potVec(0.0); 
    frameData.restraintPotential = 0.0; 
    frameData.rawPotential = 0.0;       
    frameData.volume = 0.0;          
    frameData.thermostat = make_pair(0.0, 0.0);
    frameData.electronicThermostat = make_pair(0.0, 0.0);
    frameData.barostat = Mat3x3d(0.0);              
    frameData.stressTensor = Mat3x3d(0.0);              
    frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);

    clearDerivedProperties();
  }
  
  Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups, 
                     int storageLayout) : 
    atomData(nAtoms, storageLayout), 
    rigidbodyData(nRigidbodies, storageLayout),
    cgData(nCutoffGroups, DataStorage::dslPosition),
    orthoTolerance_(1e-6) {
    
    frameData.id = -1;                   
    frameData.currentTime = 0;     
    frameData.hmat = Mat3x3d(0.0);             
    frameData.invHmat = Mat3x3d(0.0);          
    frameData.orthoRhombic = false;        
    frameData.bondPotential = 0.0;      
    frameData.bendPotential = 0.0;      
    frameData.torsionPotential = 0.0;   
    frameData.inversionPotential = 0.0; 
    frameData.lrPotentials = potVec(0.0);
    frameData.excludedPotentials = potVec(0.0); 
    frameData.restraintPotential = 0.0; 
    frameData.rawPotential = 0.0;       
    frameData.volume = 0.0;          
    frameData.thermostat = make_pair(0.0, 0.0);
    frameData.electronicThermostat = make_pair(0.0, 0.0);
    frameData.barostat = Mat3x3d(0.0);              
    frameData.stressTensor = Mat3x3d(0.0);              
    frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);

    clearDerivedProperties();
  }

  void Snapshot::clearDerivedProperties() {
    frameData.totalEnergy = 0.0;     
    frameData.translationalKinetic = 0.0;   
    frameData.rotationalKinetic = 0.0;   
    frameData.kineticEnergy = 0.0;   
    frameData.potentialEnergy = 0.0; 
    frameData.shortRangePotential = 0.0;
    frameData.longRangePotential = 0.0; 
    frameData.pressure = 0.0;        
    frameData.temperature = 0.0;
    frameData.pressureTensor = Mat3x3d(0.0);   
    frameData.systemDipole = Vector3d(0.0);            
    frameData.convectiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
    frameData.electronicTemperature = 0.0;
    frameData.COM = V3Zero;             
    frameData.COMvel = V3Zero;          
    frameData.COMw = V3Zero;            

    hasTotalEnergy = false;         
    hasTranslationalKineticEnergy = false;       
    hasRotationalKineticEnergy = false;       
    hasKineticEnergy = false;       
    hasShortRangePotential = false;
    hasLongRangePotential = false;
    hasPotentialEnergy = false;     
    hasVolume = false;         
    hasPressure = false;       
    hasTemperature = false;    
    hasElectronicTemperature = false;
    hasCOM = false;
    hasCOMvel = false;
    hasCOMw = false;
    hasPressureTensor = false;    
    hasSystemDipole = false;      
    hasConvectiveHeatFlux = false;  
    hasInertiaTensor = false;
    hasGyrationalVolume = false;   
    hasHullVolume = false;
    hasConservedQuantity = false; 
  }

  /** Returns the id of this Snapshot */
  int Snapshot::getID() {
    return frameData.id;
  }
  
  /** Sets the id of this Snapshot */
  void Snapshot::setID(int id) {
    frameData.id = id;
  }
  
  int Snapshot::getSize() {
    return atomData.getSize() + rigidbodyData.getSize();
  }
  
  /** Returns the number of atoms */
  int Snapshot::getNumberOfAtoms() {
    return atomData.getSize();
  }
  
  /** Returns the number of rigid bodies */
  int Snapshot::getNumberOfRigidBodies() {
    return rigidbodyData.getSize();
  }
  
  /** Returns the number of rigid bodies */
  int Snapshot::getNumberOfCutoffGroups() {
    return cgData.getSize();
  }
  
  /** Returns the H-Matrix */
  Mat3x3d Snapshot::getHmat() {
    return frameData.hmat;
  }

  /** Sets the H-Matrix */  
  void Snapshot::setHmat(const Mat3x3d& m) {
    hasVolume = false;
    frameData.hmat = m;
    frameData.invHmat = frameData.hmat.inverse();
    
    //determine whether the box is orthoTolerance or not
    bool oldOrthoRhombic = frameData.orthoRhombic;
    
    RealType smallDiag = fabs(frameData.hmat(0, 0));
    if(smallDiag > fabs(frameData.hmat(1, 1))) smallDiag = fabs(frameData.hmat(1, 1));
    if(smallDiag > fabs(frameData.hmat(2, 2))) smallDiag = fabs(frameData.hmat(2, 2));    
    RealType tol = smallDiag * orthoTolerance_;

    frameData.orthoRhombic = true;

    for (int i = 0; i < 3; i++ ) {
      for (int j = 0 ; j < 3; j++) {
        if (i != j) {
          if (frameData.orthoRhombic) {
            if ( fabs(frameData.hmat(i, j)) >= tol)
              frameData.orthoRhombic = false;
          }        
        }
      }
    }
    
    if( oldOrthoRhombic != frameData.orthoRhombic){
      
      if( frameData.orthoRhombic ) {
        sprintf( painCave.errMsg,
                 "OpenMD is switching from the default Non-Orthorhombic\n"
                 "\tto the faster Orthorhombic periodic boundary computations.\n"
                 "\tThis is usually a good thing, but if you want the\n"
                 "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
                 "\tvariable ( currently set to %G ) smaller.\n",
                 orthoTolerance_);
        painCave.severity = OPENMD_INFO;
        simError();
      }
      else {
        sprintf( painCave.errMsg,
                 "OpenMD is switching from the faster Orthorhombic to the more\n"
                 "\tflexible Non-Orthorhombic periodic boundary computations.\n"
                 "\tThis is usually because the box has deformed under\n"
                 "\tNPTf integration. If you want to live on the edge with\n"
                 "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
                 "\tvariable ( currently set to %G ) larger.\n",
                 orthoTolerance_);
        painCave.severity = OPENMD_WARNING;
        simError();
      }
    }    
  }
  
  /** Returns the inverse H-Matrix */
  Mat3x3d Snapshot::getInvHmat() {
    return frameData.invHmat;
  }

  RealType Snapshot::getVolume() {
    if (!hasVolume) {
      frameData.volume = frameData.hmat.determinant();
      hasVolume = true;
    }
    return frameData.volume;
  }

  void Snapshot::setVolume(RealType vol) {
    hasVolume = true;
    frameData.volume = vol;
  }

  /** Wrap a vector according to periodic boundary conditions */
  void Snapshot::wrapVector(Vector3d& pos) {
    
    Vector3d scaled = scaleVector(pos);
    
    for (int i = 0; i < 3; i++) 
      scaled[i] -= roundMe(scaled[i]);
    
    if( !frameData.orthoRhombic )
      pos = frameData.hmat * scaled;    
    else {
      
      // calc the wrapped real coordinates from the wrapped scaled coordinates
      for (int i=0; i<3; i++) {
        pos[i] = scaled[i] * frameData.hmat(i, i);
      }   
    }
  }

  /** Scaling a vector to multiples of the periodic box */
  inline Vector3d Snapshot::scaleVector(Vector3d& pos) {   
    
    Vector3d scaled;

    if( !frameData.orthoRhombic )
      scaled = frameData.invHmat * pos;
    else {
      // calc the scaled coordinates.
      for (int i=0; i<3; i++) 
        scaled[i] = pos[i] * frameData.invHmat(i, i);
    }

    return scaled;
  }

  void Snapshot::setCOM(const Vector3d& com) {
    frameData.COM = com;
    hasCOM = true;
  }
  
  void Snapshot::setCOMvel(const Vector3d& comVel) {
    frameData.COMvel = comVel;
    hasCOMvel = true;
  }
  
  void Snapshot::setCOMw(const Vector3d& comw) {
    frameData.COMw = comw;
    hasCOMw = true;
  }
  
  Vector3d Snapshot::getCOM() {
    return frameData.COM;
  }
  
  Vector3d Snapshot::getCOMvel() {
    return frameData.COMvel;
  }
  
  Vector3d Snapshot::getCOMw() {
    return frameData.COMw;
  }
  
  RealType Snapshot::getTime() {
    return frameData.currentTime;
  }
  
  void Snapshot::increaseTime(RealType dt) {
    setTime(getTime() + dt);
  }
  
  void Snapshot::setTime(RealType time) {
    frameData.currentTime = time;
  }
 
  void Snapshot::setBondPotential(RealType bp) {
    frameData.bondPotential = bp;
  }
  
  void Snapshot::setBendPotential(RealType bp) {
    frameData.bendPotential = bp;
  }
  
  void Snapshot::setTorsionPotential(RealType tp) {
    frameData.torsionPotential = tp;
  }
  
  void Snapshot::setInversionPotential(RealType ip) {
    frameData.inversionPotential = ip;
  }


  RealType Snapshot::getBondPotential() {
    return frameData.bondPotential;
  }
  RealType Snapshot::getBendPotential() {
    return frameData.bendPotential;
  }
  RealType Snapshot::getTorsionPotential() {
    return frameData.torsionPotential;
  }
  RealType Snapshot::getInversionPotential() {
    return frameData.inversionPotential;
  }

  RealType Snapshot::getShortRangePotential() {
    if (!hasShortRangePotential) {
      frameData.shortRangePotential = frameData.bondPotential;
      frameData.shortRangePotential += frameData.bendPotential;
      frameData.shortRangePotential += frameData.torsionPotential;
      frameData.shortRangePotential += frameData.inversionPotential;
      hasShortRangePotential = true;
    }
    return frameData.shortRangePotential;
  }

  void Snapshot::setLongRangePotential(potVec lrPot) {
    frameData.lrPotentials = lrPot;
  }
    
  RealType Snapshot::getLongRangePotential() {
    if (!hasLongRangePotential) {
      for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
        frameData.longRangePotential += frameData.lrPotentials[i];
      }
      hasLongRangePotential = true;
    }   
    return frameData.longRangePotential;
  }

  potVec Snapshot::getLongRangePotentials() {
    return frameData.lrPotentials;
  }

  RealType Snapshot::getPotentialEnergy() {
    if (!hasPotentialEnergy) {
      frameData.potentialEnergy = this->getLongRangePotential();
      frameData.potentialEnergy += this->getShortRangePotential();
      hasPotentialEnergy = true;
    }
    return frameData.potentialEnergy;
  }
    
  void Snapshot::setExcludedPotentials(potVec exPot) {
    frameData.excludedPotentials = exPot;
  }

  potVec Snapshot::getExcludedPotentials() {
    return frameData.excludedPotentials;
  }
      
  void Snapshot::setRestraintPotential(RealType rp) {
    frameData.restraintPotential = rp;
  }
  
  RealType Snapshot::getRestraintPotential() {
    return frameData.restraintPotential;
  }
  
  void Snapshot::setRawPotential(RealType rp) {
    frameData.rawPotential = rp;
  }
  
  RealType Snapshot::getRawPotential() {
    return frameData.rawPotential;
  }

  RealType Snapshot::getTranslationalKineticEnergy() {
    return frameData.translationalKinetic;
  }

  RealType Snapshot::getRotationalKineticEnergy() {
    return frameData.rotationalKinetic;
  }

  RealType Snapshot::getKineticEnergy() {
    return frameData.kineticEnergy;
  }

  void Snapshot::setTranslationalKineticEnergy(RealType tke) {
    hasTranslationalKineticEnergy = true;
    frameData.translationalKinetic = tke;
  }

  void Snapshot::setRotationalKineticEnergy(RealType rke) {
    hasRotationalKineticEnergy = true;
    frameData.rotationalKinetic = rke;
  }

  void Snapshot::setKineticEnergy(RealType ke) {
    hasKineticEnergy = true;
    frameData.kineticEnergy = ke;
  }

  RealType Snapshot::getTotalEnergy() {
    return frameData.totalEnergy;
  }

  void Snapshot::setTotalEnergy(RealType te) {
    hasTotalEnergy = true;
    frameData.totalEnergy = te;
  }

  RealType Snapshot::getConservedQuantity() {
    return frameData.conservedQuantity;
  }

  void Snapshot::setConservedQuantity(RealType cq) {
    hasConservedQuantity = true;
    frameData.conservedQuantity = cq;
  }

  RealType Snapshot::getTemperature() {
    return frameData.temperature;
  }

  void Snapshot::setTemperature(RealType temp) {
    hasTemperature = true;
    frameData.temperature = temp;
  }

  RealType Snapshot::getElectronicTemperature() {
    return frameData.electronicTemperature;
  }

  void Snapshot::setElectronicTemperature(RealType eTemp) {
    hasElectronicTemperature = true;
    frameData.electronicTemperature = eTemp;
  }

  RealType Snapshot::getPressure() {
    return frameData.pressure;
  }

  void Snapshot::setPressure(RealType pressure) {
    hasPressure = true;
    frameData.pressure = pressure;
  }

  Mat3x3d Snapshot::getPressureTensor() {
    return frameData.pressureTensor;
  }


  void Snapshot::setPressureTensor(const Mat3x3d& pressureTensor) {
    hasPressureTensor = true;
    frameData.pressureTensor = pressureTensor;
  }

  void Snapshot::setStressTensor(const Mat3x3d& stressTensor) {
    frameData.stressTensor = stressTensor;
  }

  Mat3x3d  Snapshot::getStressTensor() {
    return frameData.stressTensor;
  }

  void Snapshot::setConductiveHeatFlux(const Vector3d& chf) {
    frameData.conductiveHeatFlux = chf;
  }

  Vector3d Snapshot::getConductiveHeatFlux() {
    return frameData.conductiveHeatFlux;
  }
  
  Vector3d Snapshot::getConvectiveHeatFlux() {
    return frameData.convectiveHeatFlux;
  }

  void Snapshot::setConvectiveHeatFlux(const Vector3d& chf) {    
    hasConvectiveHeatFlux = true;
    frameData.convectiveHeatFlux = chf;
  }

  Vector3d Snapshot::getHeatFlux() {
    // BE CAREFUL WITH UNITS
    return getConductiveHeatFlux() + getConvectiveHeatFlux();
  }

  Vector3d Snapshot::getSystemDipole() {
    return frameData.systemDipole;
  }

  void Snapshot::setSystemDipole(const Vector3d& bd) {    
    hasSystemDipole = true;
    frameData.systemDipole = bd;
  }

  void Snapshot::setThermostat(const pair<RealType, RealType>& thermostat) {
    frameData.thermostat = thermostat;
  }

  pair<RealType, RealType> Snapshot::getThermostat() {
    return frameData.thermostat;
  }

  void Snapshot::setElectronicThermostat(const pair<RealType, RealType>& eTherm) {
    frameData.electronicThermostat = eTherm;
  }

  pair<RealType, RealType> Snapshot::getElectronicThermostat() {
    return frameData.electronicThermostat;
  }
 
  void Snapshot::setBarostat(const Mat3x3d& barostat) {
    frameData.barostat = barostat;
  }

  Mat3x3d Snapshot::getBarostat() {
    return frameData.barostat;
  }

  void Snapshot::setInertiaTensor(const Mat3x3d& inertiaTensor) {
    frameData.inertiaTensor = inertiaTensor;
    hasInertiaTensor = true;
  }

  Mat3x3d Snapshot::getInertiaTensor() {
    return frameData.inertiaTensor;
  }

  void Snapshot::setGyrationalVolume(const RealType gyrationalVolume) {
    frameData.gyrationalVolume = gyrationalVolume;
    hasGyrationalVolume = true;
  }

  RealType Snapshot::getGyrationalVolume() {
    return frameData.gyrationalVolume;
  }

  void Snapshot::setHullVolume(const RealType hullVolume) {
    frameData.hullVolume = hullVolume;
    hasHullVolume = true;
  }

  RealType Snapshot::getHullVolume() {
    return frameData.hullVolume;
  }

  void Snapshot::setOrthoTolerance(RealType ot) {
    orthoTolerance_ = ot;
  }
}
Revision:	1764
Committed:	Tue Jul 3 18:32:27 2012 UTC (12 years, 10 months ago) by gezelter
File size:	17681 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	507	/*
2	gezelter	246	* Copyright (c) 2005 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	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* 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	gezelter	1390	*
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	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43	gezelter	507	/**
44			* @file Snapshot.cpp
45			* @author tlin
46			* @date 11/11/2004
47			* @time 10:56am
48			* @version 1.0
49			*/
50	gezelter	246
51			#include "brains/Snapshot.hpp"
52			#include "utils/NumericConstant.hpp"
53			#include "utils/simError.h"
54			#include "utils/Utility.hpp"
55	jmarr	1401	#include <cstdio>
56
57	gezelter	1390	namespace OpenMD {
58	gezelter	246
59	gezelter	1764	Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups) :
60			atomData(nAtoms), rigidbodyData(nRigidbodies),
61			cgData(nCutoffGroups, DataStorage::dslPosition),
62			orthoTolerance_(1e-6) {
63
64			frameData.id = -1;
65			frameData.currentTime = 0;
66			frameData.hmat = Mat3x3d(0.0);
67			frameData.invHmat = Mat3x3d(0.0);
68			frameData.orthoRhombic = false;
69			frameData.bondPotential = 0.0;
70			frameData.bendPotential = 0.0;
71			frameData.torsionPotential = 0.0;
72			frameData.inversionPotential = 0.0;
73			frameData.lrPotentials = potVec(0.0);
74			frameData.excludedPotentials = potVec(0.0);
75			frameData.restraintPotential = 0.0;
76			frameData.rawPotential = 0.0;
77			frameData.volume = 0.0;
78			frameData.thermostat = make_pair(0.0, 0.0);
79			frameData.electronicThermostat = make_pair(0.0, 0.0);
80			frameData.barostat = Mat3x3d(0.0);
81			frameData.stressTensor = Mat3x3d(0.0);
82			frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
83
84			clearDerivedProperties();
85			}
86
87			Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups,
88			int storageLayout) :
89			atomData(nAtoms, storageLayout),
90			rigidbodyData(nRigidbodies, storageLayout),
91			cgData(nCutoffGroups, DataStorage::dslPosition),
92			orthoTolerance_(1e-6) {
93
94			frameData.id = -1;
95			frameData.currentTime = 0;
96			frameData.hmat = Mat3x3d(0.0);
97			frameData.invHmat = Mat3x3d(0.0);
98			frameData.orthoRhombic = false;
99			frameData.bondPotential = 0.0;
100			frameData.bendPotential = 0.0;
101			frameData.torsionPotential = 0.0;
102			frameData.inversionPotential = 0.0;
103			frameData.lrPotentials = potVec(0.0);
104			frameData.excludedPotentials = potVec(0.0);
105			frameData.restraintPotential = 0.0;
106			frameData.rawPotential = 0.0;
107			frameData.volume = 0.0;
108			frameData.thermostat = make_pair(0.0, 0.0);
109			frameData.electronicThermostat = make_pair(0.0, 0.0);
110			frameData.barostat = Mat3x3d(0.0);
111			frameData.stressTensor = Mat3x3d(0.0);
112			frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
113
114			clearDerivedProperties();
115			}
116
117			void Snapshot::clearDerivedProperties() {
118			frameData.totalEnergy = 0.0;
119			frameData.translationalKinetic = 0.0;
120			frameData.rotationalKinetic = 0.0;
121			frameData.kineticEnergy = 0.0;
122			frameData.potentialEnergy = 0.0;
123			frameData.shortRangePotential = 0.0;
124			frameData.longRangePotential = 0.0;
125			frameData.pressure = 0.0;
126			frameData.temperature = 0.0;
127			frameData.pressureTensor = Mat3x3d(0.0);
128			frameData.systemDipole = Vector3d(0.0);
129			frameData.convectiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
130			frameData.electronicTemperature = 0.0;
131			frameData.COM = V3Zero;
132			frameData.COMvel = V3Zero;
133			frameData.COMw = V3Zero;
134
135			hasTotalEnergy = false;
136			hasTranslationalKineticEnergy = false;
137			hasRotationalKineticEnergy = false;
138			hasKineticEnergy = false;
139			hasShortRangePotential = false;
140			hasLongRangePotential = false;
141			hasPotentialEnergy = false;
142			hasVolume = false;
143			hasPressure = false;
144			hasTemperature = false;
145			hasElectronicTemperature = false;
146			hasCOM = false;
147			hasCOMvel = false;
148			hasCOMw = false;
149			hasPressureTensor = false;
150			hasSystemDipole = false;
151			hasConvectiveHeatFlux = false;
152			hasInertiaTensor = false;
153			hasGyrationalVolume = false;
154			hasHullVolume = false;
155			hasConservedQuantity = false;
156			}
157
158			/** Returns the id of this Snapshot */
159			int Snapshot::getID() {
160			return frameData.id;
161			}
162
163			/** Sets the id of this Snapshot */
164			void Snapshot::setID(int id) {
165			frameData.id = id;
166			}
167
168			int Snapshot::getSize() {
169			return atomData.getSize() + rigidbodyData.getSize();
170			}
171
172			/** Returns the number of atoms */
173			int Snapshot::getNumberOfAtoms() {
174			return atomData.getSize();
175			}
176
177			/** Returns the number of rigid bodies */
178			int Snapshot::getNumberOfRigidBodies() {
179			return rigidbodyData.getSize();
180			}
181
182			/** Returns the number of rigid bodies */
183			int Snapshot::getNumberOfCutoffGroups() {
184			return cgData.getSize();
185			}
186
187			/** Returns the H-Matrix */
188			Mat3x3d Snapshot::getHmat() {
189			return frameData.hmat;
190			}
191
192			/** Sets the H-Matrix */
193			void Snapshot::setHmat(const Mat3x3d& m) {
194			hasVolume = false;
195	gezelter	1715	frameData.hmat = m;
196			frameData.invHmat = frameData.hmat.inverse();
197	gezelter	246
198			//determine whether the box is orthoTolerance or not
199	gezelter	1715	bool oldOrthoRhombic = frameData.orthoRhombic;
200	gezelter	246
201	gezelter	1715	RealType smallDiag = fabs(frameData.hmat(0, 0));
202			if(smallDiag > fabs(frameData.hmat(1, 1))) smallDiag = fabs(frameData.hmat(1, 1));
203			if(smallDiag > fabs(frameData.hmat(2, 2))) smallDiag = fabs(frameData.hmat(2, 2));
204	gezelter	1021	RealType tol = smallDiag * orthoTolerance_;
205	gezelter	246
206	gezelter	1715	frameData.orthoRhombic = true;
207	gezelter	246
208			for (int i = 0; i < 3; i++ ) {
209	gezelter	507	for (int j = 0 ; j < 3; j++) {
210			if (i != j) {
211	gezelter	1715	if (frameData.orthoRhombic) {
212			if ( fabs(frameData.hmat(i, j)) >= tol)
213			frameData.orthoRhombic = false;
214	gezelter	507	}
215			}
216			}
217	gezelter	246	}
218	gezelter	1764
219	gezelter	1715	if( oldOrthoRhombic != frameData.orthoRhombic){
220	gezelter	1764
221	gezelter	1715	if( frameData.orthoRhombic ) {
222	gezelter	507	sprintf( painCave.errMsg,
223	gezelter	1390	"OpenMD is switching from the default Non-Orthorhombic\n"
224	gezelter	507	"\tto the faster Orthorhombic periodic boundary computations.\n"
225	gezelter	890	"\tThis is usually a good thing, but if you want the\n"
226	gezelter	507	"\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
227			"\tvariable ( currently set to %G ) smaller.\n",
228	gezelter	1021	orthoTolerance_);
229	gezelter	1390	painCave.severity = OPENMD_INFO;
230	gezelter	507	simError();
231			}
232			else {
233			sprintf( painCave.errMsg,
234	gezelter	1390	"OpenMD is switching from the faster Orthorhombic to the more\n"
235	gezelter	507	"\tflexible Non-Orthorhombic periodic boundary computations.\n"
236			"\tThis is usually because the box has deformed under\n"
237	gezelter	890	"\tNPTf integration. If you want to live on the edge with\n"
238	gezelter	507	"\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
239			"\tvariable ( currently set to %G ) larger.\n",
240	gezelter	1021	orthoTolerance_);
241	gezelter	1390	painCave.severity = OPENMD_WARNING;
242	gezelter	507	simError();
243			}
244	gezelter	246	}
245	gezelter	507	}
246	gezelter	1764
247			/** Returns the inverse H-Matrix */
248			Mat3x3d Snapshot::getInvHmat() {
249			return frameData.invHmat;
250			}
251	gezelter	246
252	gezelter	1764	RealType Snapshot::getVolume() {
253			if (!hasVolume) {
254			frameData.volume = frameData.hmat.determinant();
255			hasVolume = true;
256			}
257			return frameData.volume;
258			}
259	gezelter	246
260	gezelter	1764	void Snapshot::setVolume(RealType vol) {
261			hasVolume = true;
262			frameData.volume = vol;
263			}
264
265			/** Wrap a vector according to periodic boundary conditions */
266	gezelter	1568	void Snapshot::wrapVector(Vector3d& pos) {
267	gezelter	1562
268			Vector3d scaled = scaleVector(pos);
269	gezelter	1568
270	gezelter	1562	for (int i = 0; i < 3; i++)
271			scaled[i] -= roundMe(scaled[i]);
272	gezelter	1764
273	gezelter	1715	if( !frameData.orthoRhombic )
274			pos = frameData.hmat * scaled;
275	gezelter	1562	else {
276	gezelter	1764
277	gezelter	507	// calc the wrapped real coordinates from the wrapped scaled coordinates
278	gezelter	1562	for (int i=0; i<3; i++) {
279	gezelter	1715	pos[i] = scaled[i] * frameData.hmat(i, i);
280	gezelter	1540	}
281	gezelter	246	}
282	gezelter	507	}
283	gezelter	246
284	gezelter	1764	/** Scaling a vector to multiples of the periodic box */
285	gezelter	1562	inline Vector3d Snapshot::scaleVector(Vector3d& pos) {
286
287			Vector3d scaled;
288
289	gezelter	1715	if( !frameData.orthoRhombic )
290			scaled = frameData.invHmat * pos;
291	gezelter	1562	else {
292			// calc the scaled coordinates.
293			for (int i=0; i<3; i++)
294	gezelter	1715	scaled[i] = pos[i] * frameData.invHmat(i, i);
295	gezelter	1562	}
296
297			return scaled;
298			}
299	gezelter	1764
300			void Snapshot::setCOM(const Vector3d& com) {
301			frameData.COM = com;
302			hasCOM = true;
303			}
304	gezelter	1562
305	gezelter	1764	void Snapshot::setCOMvel(const Vector3d& comVel) {
306			frameData.COMvel = comVel;
307			hasCOMvel = true;
308			}
309
310			void Snapshot::setCOMw(const Vector3d& comw) {
311			frameData.COMw = comw;
312			hasCOMw = true;
313			}
314
315	gezelter	1104	Vector3d Snapshot::getCOM() {
316	gezelter	1715	return frameData.COM;
317	gezelter	1104	}
318
319			Vector3d Snapshot::getCOMvel() {
320	gezelter	1715	return frameData.COMvel;
321	gezelter	1104	}
322
323			Vector3d Snapshot::getCOMw() {
324	gezelter	1715	return frameData.COMw;
325	gezelter	1104	}
326	gezelter	1764
327			RealType Snapshot::getTime() {
328			return frameData.currentTime;
329			}
330
331			void Snapshot::increaseTime(RealType dt) {
332			setTime(getTime() + dt);
333			}
334
335			void Snapshot::setTime(RealType time) {
336			frameData.currentTime = time;
337			}
338
339			void Snapshot::setBondPotential(RealType bp) {
340			frameData.bondPotential = bp;
341			}
342
343			void Snapshot::setBendPotential(RealType bp) {
344			frameData.bendPotential = bp;
345			}
346
347			void Snapshot::setTorsionPotential(RealType tp) {
348			frameData.torsionPotential = tp;
349			}
350
351			void Snapshot::setInversionPotential(RealType ip) {
352			frameData.inversionPotential = ip;
353			}
354
355
356			RealType Snapshot::getBondPotential() {
357			return frameData.bondPotential;
358			}
359			RealType Snapshot::getBendPotential() {
360			return frameData.bendPotential;
361			}
362			RealType Snapshot::getTorsionPotential() {
363			return frameData.torsionPotential;
364			}
365			RealType Snapshot::getInversionPotential() {
366			return frameData.inversionPotential;
367			}
368
369			RealType Snapshot::getShortRangePotential() {
370			if (!hasShortRangePotential) {
371			frameData.shortRangePotential = frameData.bondPotential;
372			frameData.shortRangePotential += frameData.bendPotential;
373			frameData.shortRangePotential += frameData.torsionPotential;
374			frameData.shortRangePotential += frameData.inversionPotential;
375			hasShortRangePotential = true;
376			}
377			return frameData.shortRangePotential;
378			}
379
380			void Snapshot::setLongRangePotential(potVec lrPot) {
381			frameData.lrPotentials = lrPot;
382			}
383
384			RealType Snapshot::getLongRangePotential() {
385			if (!hasLongRangePotential) {
386			for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
387			frameData.longRangePotential += frameData.lrPotentials[i];
388			}
389			hasLongRangePotential = true;
390			}
391			return frameData.longRangePotential;
392			}
393
394			potVec Snapshot::getLongRangePotentials() {
395			return frameData.lrPotentials;
396			}
397
398			RealType Snapshot::getPotentialEnergy() {
399			if (!hasPotentialEnergy) {
400			frameData.potentialEnergy = this->getLongRangePotential();
401			frameData.potentialEnergy += this->getShortRangePotential();
402			hasPotentialEnergy = true;
403			}
404			return frameData.potentialEnergy;
405			}
406
407			void Snapshot::setExcludedPotentials(potVec exPot) {
408			frameData.excludedPotentials = exPot;
409			}
410
411			potVec Snapshot::getExcludedPotentials() {
412			return frameData.excludedPotentials;
413			}
414
415			void Snapshot::setRestraintPotential(RealType rp) {
416			frameData.restraintPotential = rp;
417			}
418
419			RealType Snapshot::getRestraintPotential() {
420			return frameData.restraintPotential;
421			}
422
423			void Snapshot::setRawPotential(RealType rp) {
424			frameData.rawPotential = rp;
425			}
426
427			RealType Snapshot::getRawPotential() {
428			return frameData.rawPotential;
429			}
430
431			RealType Snapshot::getTranslationalKineticEnergy() {
432			return frameData.translationalKinetic;
433			}
434
435			RealType Snapshot::getRotationalKineticEnergy() {
436			return frameData.rotationalKinetic;
437			}
438
439			RealType Snapshot::getKineticEnergy() {
440			return frameData.kineticEnergy;
441			}
442
443			void Snapshot::setTranslationalKineticEnergy(RealType tke) {
444			hasTranslationalKineticEnergy = true;
445			frameData.translationalKinetic = tke;
446			}
447
448			void Snapshot::setRotationalKineticEnergy(RealType rke) {
449			hasRotationalKineticEnergy = true;
450			frameData.rotationalKinetic = rke;
451			}
452
453			void Snapshot::setKineticEnergy(RealType ke) {
454			hasKineticEnergy = true;
455			frameData.kineticEnergy = ke;
456			}
457
458			RealType Snapshot::getTotalEnergy() {
459			return frameData.totalEnergy;
460			}
461
462			void Snapshot::setTotalEnergy(RealType te) {
463			hasTotalEnergy = true;
464			frameData.totalEnergy = te;
465			}
466
467			RealType Snapshot::getConservedQuantity() {
468			return frameData.conservedQuantity;
469			}
470
471			void Snapshot::setConservedQuantity(RealType cq) {
472			hasConservedQuantity = true;
473			frameData.conservedQuantity = cq;
474			}
475
476			RealType Snapshot::getTemperature() {
477			return frameData.temperature;
478			}
479
480			void Snapshot::setTemperature(RealType temp) {
481			hasTemperature = true;
482			frameData.temperature = temp;
483			}
484
485			RealType Snapshot::getElectronicTemperature() {
486			return frameData.electronicTemperature;
487			}
488
489			void Snapshot::setElectronicTemperature(RealType eTemp) {
490			hasElectronicTemperature = true;
491			frameData.electronicTemperature = eTemp;
492			}
493
494			RealType Snapshot::getPressure() {
495			return frameData.pressure;
496			}
497
498			void Snapshot::setPressure(RealType pressure) {
499			hasPressure = true;
500			frameData.pressure = pressure;
501			}
502
503			Mat3x3d Snapshot::getPressureTensor() {
504			return frameData.pressureTensor;
505			}
506
507
508			void Snapshot::setPressureTensor(const Mat3x3d& pressureTensor) {
509			hasPressureTensor = true;
510			frameData.pressureTensor = pressureTensor;
511			}
512
513			void Snapshot::setStressTensor(const Mat3x3d& stressTensor) {
514			frameData.stressTensor = stressTensor;
515			}
516
517			Mat3x3d Snapshot::getStressTensor() {
518			return frameData.stressTensor;
519			}
520
521			void Snapshot::setConductiveHeatFlux(const Vector3d& chf) {
522			frameData.conductiveHeatFlux = chf;
523			}
524
525			Vector3d Snapshot::getConductiveHeatFlux() {
526			return frameData.conductiveHeatFlux;
527			}
528
529			Vector3d Snapshot::getConvectiveHeatFlux() {
530			return frameData.convectiveHeatFlux;
531			}
532
533			void Snapshot::setConvectiveHeatFlux(const Vector3d& chf) {
534			hasConvectiveHeatFlux = true;
535			frameData.convectiveHeatFlux = chf;
536			}
537
538			Vector3d Snapshot::getHeatFlux() {
539			// BE CAREFUL WITH UNITS
540			return getConductiveHeatFlux() + getConvectiveHeatFlux();
541			}
542
543			Vector3d Snapshot::getSystemDipole() {
544			return frameData.systemDipole;
545			}
546
547			void Snapshot::setSystemDipole(const Vector3d& bd) {
548			hasSystemDipole = true;
549			frameData.systemDipole = bd;
550			}
551
552			void Snapshot::setThermostat(const pair<RealType, RealType>& thermostat) {
553			frameData.thermostat = thermostat;
554			}
555
556			pair<RealType, RealType> Snapshot::getThermostat() {
557			return frameData.thermostat;
558			}
559
560			void Snapshot::setElectronicThermostat(const pair<RealType, RealType>& eTherm) {
561			frameData.electronicThermostat = eTherm;
562			}
563
564			pair<RealType, RealType> Snapshot::getElectronicThermostat() {
565			return frameData.electronicThermostat;
566			}
567
568			void Snapshot::setBarostat(const Mat3x3d& barostat) {
569			frameData.barostat = barostat;
570			}
571
572			Mat3x3d Snapshot::getBarostat() {
573			return frameData.barostat;
574			}
575
576			void Snapshot::setInertiaTensor(const Mat3x3d& inertiaTensor) {
577			frameData.inertiaTensor = inertiaTensor;
578			hasInertiaTensor = true;
579			}
580
581			Mat3x3d Snapshot::getInertiaTensor() {
582			return frameData.inertiaTensor;
583			}
584
585			void Snapshot::setGyrationalVolume(const RealType gyrationalVolume) {
586			frameData.gyrationalVolume = gyrationalVolume;
587			hasGyrationalVolume = true;
588			}
589
590			RealType Snapshot::getGyrationalVolume() {
591			return frameData.gyrationalVolume;
592			}
593
594			void Snapshot::setHullVolume(const RealType hullVolume) {
595			frameData.hullVolume = hullVolume;
596			hasHullVolume = true;
597			}
598
599			RealType Snapshot::getHullVolume() {
600			return frameData.hullVolume;
601			}
602
603			void Snapshot::setOrthoTolerance(RealType ot) {
604			orthoTolerance_ = ot;
605			}
606	gezelter	246	}