--- trunk/src/brains/Snapshot.cpp 2006/02/16 22:05:48 890 +++ trunk/src/brains/Snapshot.cpp 2012/08/22 02:28:28 1782 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 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. @@ -37,6 +28,16 @@ * 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). */ /** @@ -51,110 +52,568 @@ #include "utils/NumericConstant.hpp" #include "utils/simError.h" #include "utils/Utility.hpp" -namespace oopse { +#include - void Snapshot::setHmat(const Mat3x3d& m) { - const double orthoTolerance = NumericConstant::epsilon; - hmat_ = m; - invHmat_ = hmat_.inverse(); +namespace OpenMD { + + Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups) : + atomData(nAtoms), rigidbodyData(nRigidbodies), + cgData(nCutoffGroups, DataStorage::dslPosition), + orthoTolerance_(1e-6) { - //prepare fortran Hmat - double fortranHmat[9]; - double fortranInvHmat[9]; - hmat_.getArray(fortranHmat); - invHmat_.getArray(fortranInvHmat); + 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 - int oldOrthoRhombic = orthoRhombic_; + bool oldOrthoRhombic = frameData.orthoRhombic; - double smallDiag = fabs(hmat_(0, 0)); - if(smallDiag > fabs(hmat_(1, 1))) smallDiag = fabs(hmat_(1, 1)); - if(smallDiag > fabs(hmat_(2, 2))) smallDiag = fabs(hmat_(2, 2)); - double tol = smallDiag * orthoTolerance; + 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_; - orthoRhombic_ = 1; + frameData.orthoRhombic = true; for (int i = 0; i < 3; i++ ) { for (int j = 0 ; j < 3; j++) { if (i != j) { - if (orthoRhombic_) { - if ( fabs(hmat_(i, j)) >= tol) - orthoRhombic_ = 0; + if (frameData.orthoRhombic) { + if ( fabs(frameData.hmat(i, j)) >= tol) + frameData.orthoRhombic = false; } } } } - - if( oldOrthoRhombic != orthoRhombic_ ){ - - if( orthoRhombic_ ) { + + if( oldOrthoRhombic != frameData.orthoRhombic){ + + if( frameData.orthoRhombic ) { sprintf( painCave.errMsg, - "OOPSE is switching from the default Non-Orthorhombic\n" + "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 = OOPSE_INFO; + orthoTolerance_); + painCave.severity = OPENMD_INFO; simError(); } else { sprintf( painCave.errMsg, - "OOPSE is switching from the faster Orthorhombic to the more\n" + "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 = OOPSE_WARNING; + orthoTolerance_); + painCave.severity = OPENMD_WARNING; simError(); } } + } + + /** Returns the inverse H-Matrix */ + Mat3x3d Snapshot::getInvHmat() { + return frameData.invHmat; + } - //notify fortran simulation box has changed - setFortranBox(fortranHmat, fortranInvHmat, &orthoRhombic_); + 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); + } + } + } - int i; + /** Scaling a vector to multiples of the periodic box */ + inline Vector3d Snapshot::scaleVector(Vector3d& pos) { + Vector3d scaled; - if( !orthoRhombic_ ){ - + if( !frameData.orthoRhombic ) + scaled = frameData.invHmat * pos; + else { // calc the scaled coordinates. - scaled = invHmat_* pos; + for (int i=0; i<3; i++) + scaled[i] = pos[i] * frameData.invHmat(i, i); + } - // wrap the scaled coordinates - for (i = 0; i < 3; ++i) { - scaled[i] -= roundMe(scaled[i]); - } + return scaled; + } - // calc the wrapped real coordinates from the wrapped scaled coordinates - pos = hmat_ * scaled; - - } else {//if it is orthoRhombic, we could improve efficiency by only caculating the diagonal element + 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; + } - // calc the scaled coordinates. - for (i=0; i<3; i++) { - scaled[i] = pos[i] * invHmat_(i, i); + RealType Snapshot::getLongRangePotential() { + if (!hasLongRangePotential) { + for (int i = 0; i < N_INTERACTION_FAMILIES; i++) { + frameData.longRangePotential += frameData.lrPotentials[i]; } - - // wrap the scaled coordinates - for (i = 0; i < 3; ++i) { - scaled[i] -= roundMe(scaled[i]); - } + hasLongRangePotential = true; + } + return frameData.longRangePotential; + } - // calc the wrapped real coordinates from the wrapped scaled coordinates - for (i=0; i<3; i++) { - pos[i] = scaled[i] * hmat_(i, i); - } - + 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& thermostat) { + frameData.thermostat = thermostat; + } + + pair Snapshot::getThermostat() { + return frameData.thermostat; + } + + void Snapshot::setElectronicThermostat(const pair& eTherm) { + frameData.electronicThermostat = eTherm; + } + + pair 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; + } +}