src/brains/Snapshot.cpp

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

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

namespace OpenMD {

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

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

    clearDerivedProperties();
  }

  void Snapshot::clearDerivedProperties() {
    frameData.totalEnergy = 0.0;     
    frameData.translationalKinetic = 0.0;   
    frameData.rotationalKinetic = 0.0;   
    frameData.kineticEnergy = 0.0;   
    frameData.potentialEnergy = 0.0; 
    frameData.shortRangePotential = 0.0;
    frameData.longRangePotential = 0.0; 
    frameData.pressure = 0.0;        
    frameData.temperature = 0.0;
    frameData.pressureTensor = Mat3x3d(0.0);   
    frameData.systemDipole = Vector3d(0.0);            
    frameData.convectiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
    frameData.electronicTemperature = 0.0;
    frameData.COM = V3Zero;             
    frameData.COMvel = V3Zero;          
    frameData.COMw = V3Zero;            

    hasTotalEnergy = false;         
    hasTranslationalKineticEnergy = false;       
    hasRotationalKineticEnergy = false;       
    hasKineticEnergy = false;       
    hasShortRangePotential = false;
    hasLongRangePotential = false;
    hasPotentialEnergy = false;   
    hasXYarea = false;
    hasVolume = false;         
    hasPressure = false;       
    hasTemperature = false;    
    hasElectronicTemperature = false;
    hasCOM = false;
    hasCOMvel = false;
    hasCOMw = false;
    hasPressureTensor = false;    
    hasSystemDipole = false;      
    hasConvectiveHeatFlux = false;  
    hasInertiaTensor = false;
    hasGyrationalVolume = false;   
    hasHullVolume = false;
    hasConservedQuantity = false; 
  }

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

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

    frameData.orthoRhombic = true;

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

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

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

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

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

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

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

    return scaled;
  }

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

  void Snapshot::setThermostat(const pair<RealType, RealType>& thermostat) {
    frameData.thermostat = thermostat;
  }

  pair<RealType, RealType> Snapshot::getThermostat() {
    return frameData.thermostat;
  }

  void Snapshot::setElectronicThermostat(const pair<RealType, RealType>& eTherm) {
    frameData.electronicThermostat = eTherm;
  }

  pair<RealType, RealType> Snapshot::getElectronicThermostat() {
    return frameData.electronicThermostat;
  }
 
  void Snapshot::setBarostat(const Mat3x3d& barostat) {
    frameData.barostat = barostat;
  }

  Mat3x3d Snapshot::getBarostat() {
    return frameData.barostat;
  }

  void Snapshot::setInertiaTensor(const Mat3x3d& inertiaTensor) {
    frameData.inertiaTensor = inertiaTensor;
    hasInertiaTensor = true;
  }

  Mat3x3d Snapshot::getInertiaTensor() {
    return frameData.inertiaTensor;
  }

  void Snapshot::setGyrationalVolume(const RealType gyrationalVolume) {
    frameData.gyrationalVolume = gyrationalVolume;
    hasGyrationalVolume = true;
  }

  RealType Snapshot::getGyrationalVolume() {
    return frameData.gyrationalVolume;
  }

  void Snapshot::setHullVolume(const RealType hullVolume) {
    frameData.hullVolume = hullVolume;
    hasHullVolume = true;
  }

  RealType Snapshot::getHullVolume() {
    return frameData.hullVolume;
  }

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