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