src/brains/Snapshot.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
  
/**
 * @file Snapshot.cpp
 * @author tlin
 * @date 11/11/2004
 * @version 1.0
 */

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

namespace OpenMD {

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

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

    clearDerivedProperties();
  }

  void Snapshot::clearDerivedProperties() {
    frameData.totalEnergy = 0.0;     
    frameData.translationalKinetic = 0.0;   
    frameData.rotationalKinetic = 0.0;   
    frameData.kineticEnergy = 0.0;   
    frameData.potentialEnergy = 0.0; 
    frameData.shortRangePotential = 0.0;
    frameData.longRangePotential = 0.0; 
    frameData.pressure = 0.0;        
    frameData.temperature = 0.0;
    frameData.pressureTensor = Mat3x3d(0.0);   
    frameData.systemDipole = Vector3d(0.0);
    frameData.systemQuadrupole = Mat3x3d(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;      
    hasSystemQuadrupole = false;
    hasConvectiveHeatFlux = false;  
    hasInertiaTensor = false;
    hasGyrationalVolume = false;   
    hasHullVolume = false;
    hasConservedQuantity = false; 
    hasBoundingBox = false;
  }

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

    /** Returns the number of bytes in a FrameData structure */
  int Snapshot::getFrameDataSize() {
    return sizeof(FrameData);
  }
  
  /** Returns the H-Matrix */
  Mat3x3d Snapshot::getHmat() {
    return frameData.hmat;
  }

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

    frameData.orthoRhombic = true;

    for (int i = 0; i < 3; i++ ) {
      for (int j = 0 ; j < 3; j++) {
        if (i != j) {
          if (frameData.orthoRhombic) {
            if ( fabs(frameData.hmat(i, j)) >= tol)
              frameData.orthoRhombic = false;
          }        
        }
      }
    }
    
    if( oldOrthoRhombic != frameData.orthoRhombic){
      
      // It is finally time to suppress these warnings once and for
      // all.  They were annoying and not very informative.

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

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

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

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

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

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

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


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

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

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

    return scaled;
  }

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


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

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

  void Snapshot::setReciprocalPotential(RealType rp){
    frameData.reciprocalPotential = rp;
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

  Mat3x3d Snapshot::getSystemQuadrupole() {
    return frameData.systemQuadrupole;
  }

  void Snapshot::setSystemQuadrupole(const Mat3x3d& bq) {    
    hasSystemQuadrupole = true;
    frameData.systemQuadrupole = bq;
  }

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