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
#	User	Rev	Content
1	gezelter	507	/*
2	gezelter	246	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	gezelter	246	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43	gezelter	507	/**
44			* @file Snapshot.cpp
45			* @author tlin
46			* @date 11/11/2004
47			* @time 10:56am
48			* @version 1.0
49			*/
50	gezelter	246
51			#include "brains/Snapshot.hpp"
52			#include "utils/NumericConstant.hpp"
53			#include "utils/simError.h"
54			#include "utils/Utility.hpp"
55	jmarr	1401	#include <cstdio>
56
57	gezelter	1390	namespace OpenMD {
58	gezelter	246
59	gezelter	1782	Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups) :
60			atomData(nAtoms), rigidbodyData(nRigidbodies),
61			cgData(nCutoffGroups, DataStorage::dslPosition),
62			orthoTolerance_(1e-6) {
63	gezelter	246
64	gezelter	1782	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	chuckv	1112
85	gezelter	1782	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	gezelter	246
116	gezelter	1782	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	gezelter	246	//determine whether the box is orthoTolerance or not
202	gezelter	1782	bool oldOrthoRhombic = frameData.orthoRhombic;
203	gezelter	246
204	gezelter	1782	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	gezelter	1021	RealType tol = smallDiag * orthoTolerance_;
208	gezelter	246
209	gezelter	1782	frameData.orthoRhombic = true;
210	gezelter	246
211			for (int i = 0; i < 3; i++ ) {
212	gezelter	507	for (int j = 0 ; j < 3; j++) {
213			if (i != j) {
214	gezelter	1782	if (frameData.orthoRhombic) {
215			if ( fabs(frameData.hmat(i, j)) >= tol)
216			frameData.orthoRhombic = false;
217	gezelter	507	}
218			}
219			}
220	gezelter	246	}
221	gezelter	1782
222			if( oldOrthoRhombic != frameData.orthoRhombic){
223
224			if( frameData.orthoRhombic ) {
225	gezelter	507	sprintf( painCave.errMsg,
226	gezelter	1390	"OpenMD is switching from the default Non-Orthorhombic\n"
227	gezelter	507	"\tto the faster Orthorhombic periodic boundary computations.\n"
228	gezelter	890	"\tThis is usually a good thing, but if you want the\n"
229	gezelter	507	"\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
230			"\tvariable ( currently set to %G ) smaller.\n",
231	gezelter	1021	orthoTolerance_);
232	gezelter	1390	painCave.severity = OPENMD_INFO;
233	gezelter	507	simError();
234			}
235			else {
236			sprintf( painCave.errMsg,
237	gezelter	1390	"OpenMD is switching from the faster Orthorhombic to the more\n"
238	gezelter	507	"\tflexible Non-Orthorhombic periodic boundary computations.\n"
239			"\tThis is usually because the box has deformed under\n"
240	gezelter	890	"\tNPTf integration. If you want to live on the edge with\n"
241	gezelter	507	"\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
242			"\tvariable ( currently set to %G ) larger.\n",
243	gezelter	1021	orthoTolerance_);
244	gezelter	1390	painCave.severity = OPENMD_WARNING;
245	gezelter	507	simError();
246			}
247	gezelter	246	}
248	gezelter	1782	}
249
250			/** Returns the inverse H-Matrix */
251			Mat3x3d Snapshot::getInvHmat() {
252			return frameData.invHmat;
253			}
254	gezelter	246
255	gezelter	1782	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	gezelter	507	}
264	gezelter	246
265	gezelter	1782	RealType Snapshot::getVolume() {
266			if (!hasVolume) {
267			frameData.volume = frameData.hmat.determinant();
268			hasVolume = true;
269			}
270			return frameData.volume;
271			}
272	gezelter	246
273	gezelter	1782	void Snapshot::setVolume(RealType vol) {
274			hasVolume = true;
275			frameData.volume = vol;
276			}
277
278			/** Wrap a vector according to periodic boundary conditions */
279	gezelter	507	void Snapshot::wrapVector(Vector3d& pos) {
280	gezelter	1782
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	gezelter	246
297	gezelter	1782	/** Scaling a vector to multiples of the periodic box */
298			inline Vector3d Snapshot::scaleVector(Vector3d& pos) {
299
300	gezelter	246	Vector3d scaled;
301
302	gezelter	1782	if( !frameData.orthoRhombic )
303			scaled = frameData.invHmat * pos;
304			else {
305	gezelter	507	// calc the scaled coordinates.
306	gezelter	1782	for (int i=0; i<3; i++)
307			scaled[i] = pos[i] * frameData.invHmat(i, i);
308			}
309	gezelter	246
310	gezelter	1782	return scaled;
311			}
312	gezelter	246
313	gezelter	1782	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	gezelter	246
368
369	gezelter	1782	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	gezelter	246	}
390	gezelter	1782	return frameData.shortRangePotential;
391			}
392	gezelter	246
393	gezelter	1782	void Snapshot::setLongRangePotential(potVec lrPot) {
394			frameData.lrPotentials = lrPot;
395	gezelter	507	}
396	gezelter	1782
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	gezelter	246
407	gezelter	1782	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	gezelter	1104	}
417	gezelter	1782	return frameData.potentialEnergy;
418	gezelter	1104	}
419	gezelter	1782
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	gezelter	1104
432	gezelter	1782	RealType Snapshot::getRestraintPotential() {
433			return frameData.restraintPotential;
434	gezelter	1104	}
435
436	gezelter	1782	void Snapshot::setRawPotential(RealType rp) {
437			frameData.rawPotential = rp;
438	gezelter	1104	}
439	gezelter	1782
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	gezelter	1104
581	gezelter	1782	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	gezelter	246	}