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, 234107 (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
 * @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.reciprocalPotential = 0.0;
    frameData.excludedPotentials = potVec(0.0); 
    frameData.restraintPotential = 0.0; 
    frameData.rawPotential = 0.0;   
    frameData.xyArea = 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.bBox = Mat3x3d(0.0);             
    frameData.invBbox = 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.reciprocalPotential = 0.0;
    frameData.excludedPotentials = potVec(0.0); 
    frameData.restraintPotential = 0.0; 
    frameData.rawPotential = 0.0;       
    frameData.xyArea = 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;   
    hasXYarea = 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; 
    hasBoundingBox = 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 number of bytes in a FrameData structure */
  int Snapshot::getFrameDataSize() {
    return sizeof(FrameData);
  }
  
  /** 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){
      
      // It is finally time to suppress these warnings once and for
      // all.  They were annoying and not very informative.

      // 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;
  }

  /** Returns the Bounding Box */
  Mat3x3d Snapshot::getBoundingBox() {
    return frameData.bBox;
  }

  /** Sets the Bounding Box */  
  void Snapshot::setBoundingBox(const Mat3x3d& m) {
    frameData.bBox = m;
    frameData.invBbox = frameData.bBox.inverse();
    hasBoundingBox = true;
  }

  /** Returns the inverse Bounding Box */
  Mat3x3d Snapshot::getInvBoundingBox() {
    return frameData.invBbox;
  }

  RealType Snapshot::getXYarea() {
    if (!hasXYarea) {
      Vector3d x = frameData.hmat.getColumn(0);
      Vector3d y = frameData.hmat.getColumn(1);
      frameData.xyArea = cross(x,y).length();
      hasXYarea = true;
    }
    return frameData.xyArea;
  }

  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) {
    
    if( !frameData.orthoRhombic ) {
      Vector3d scaled = frameData.invHmat * pos;
      for (int i = 0; i < 3; i++) {
        scaled[i] -= roundMe( scaled[i] );        
      }
      // calc the wrapped real coordinates from the wrapped scaled coordinates
      pos = frameData.hmat * scaled;
    } else {
      RealType scaled;
      for (int i=0; i<3; i++) {      
        scaled = pos[i] * frameData.invHmat(i,i);
        scaled -= roundMe( scaled );
        pos[i] = scaled * 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::setReciprocalPotential(RealType rp){
    frameData.reciprocalPotential = rp;
  }

  RealType Snapshot::getReciprocalPotential() {
    return frameData.reciprocalPotential;
  }

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