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
 * @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.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.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.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; 
  }

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