55#include "brains/Snapshot.hpp"
59#include "utils/Utility.hpp"
60#include "utils/simError.h"
64 Snapshot::Snapshot(
int nAtoms,
int nRigidbodies,
int nCutoffGroups,
67 rigidbodyData(nRigidbodies),
68 cgData(nCutoffGroups,
DataStorage::dslPosition), orthoTolerance_(1e-6) {
70 frameData.currentTime = 0;
71 frameData.hmat = Mat3x3d(0.0);
72 frameData.invHmat = Mat3x3d(0.0);
73 frameData.orthoRhombic =
false;
74 frameData.usePBC = usePBC;
75 frameData.bondPotential = 0.0;
76 frameData.bendPotential = 0.0;
77 frameData.torsionPotential = 0.0;
78 frameData.inversionPotential = 0.0;
79 frameData.lrPotentials = potVec(0.0);
80 frameData.surfacePotential = 0.0;
81 frameData.reciprocalPotential = 0.0;
82 frameData.selfPotentials = potVec(0.0);
83 frameData.excludedPotentials = potVec(0.0);
84 frameData.restraintPotential = 0.0;
85 frameData.rawPotential = 0.0;
86 frameData.xyArea = 0.0;
87 frameData.xzArea = 0.0;
88 frameData.yzArea = 0.0;
89 frameData.volume = 0.0;
90 frameData.thermostat = make_pair(0.0, 0.0);
91 frameData.electronicThermostat = make_pair(0.0, 0.0);
92 frameData.barostat = Mat3x3d(0.0);
93 frameData.virialTensor = Mat3x3d(0.0);
94 frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
95 frameData.spfData = std::make_shared<SPFData>();
96 clearDerivedProperties();
99 Snapshot::Snapshot(
int nAtoms,
int nRigidbodies,
int nCutoffGroups,
100 int atomStorageLayout,
int rigidBodyStorageLayout,
101 int cutoffGroupStorageLayout,
bool usePBC) :
102 atomData(nAtoms, atomStorageLayout),
103 rigidbodyData(nRigidbodies, rigidBodyStorageLayout),
104 cgData(nCutoffGroups, cutoffGroupStorageLayout), orthoTolerance_(1e-6) {
106 frameData.currentTime = 0;
107 frameData.hmat = Mat3x3d(0.0);
108 frameData.invHmat = Mat3x3d(0.0);
109 frameData.bBox = Mat3x3d(0.0);
110 frameData.invBbox = Mat3x3d(0.0);
111 frameData.orthoRhombic =
false;
112 frameData.usePBC = usePBC;
113 frameData.bondPotential = 0.0;
114 frameData.bendPotential = 0.0;
115 frameData.torsionPotential = 0.0;
116 frameData.inversionPotential = 0.0;
117 frameData.lrPotentials = potVec(0.0);
118 frameData.surfacePotential = 0.0;
119 frameData.reciprocalPotential = 0.0;
120 frameData.selfPotentials = potVec(0.0);
121 frameData.excludedPotentials = potVec(0.0);
122 frameData.restraintPotential = 0.0;
123 frameData.rawPotential = 0.0;
124 frameData.xyArea = 0.0;
125 frameData.xzArea = 0.0;
126 frameData.yzArea = 0.0;
127 frameData.volume = 0.0;
128 frameData.thermostat = make_pair(0.0, 0.0);
129 frameData.electronicThermostat = make_pair(0.0, 0.0);
130 frameData.barostat = Mat3x3d(0.0);
131 frameData.virialTensor = Mat3x3d(0.0);
132 frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
133 frameData.spfData = std::make_shared<SPFData>();
135 clearDerivedProperties();
139 frameData.totalEnergy = 0.0;
140 frameData.translationalKinetic = 0.0;
141 frameData.rotationalKinetic = 0.0;
142 frameData.electronicKinetic = 0.0;
143 frameData.kineticEnergy = 0.0;
144 frameData.potentialEnergy = 0.0;
145 frameData.shortRangePotential = 0.0;
146 frameData.longRangePotential = 0.0;
147 frameData.excludedPotential = 0.0;
148 frameData.selfPotential = 0.0;
149 frameData.pressure = 0.0;
150 frameData.temperature = 0.0;
151 frameData.pressureTensor = Mat3x3d(0.0);
152 frameData.systemDipole = Vector3d(0.0);
153 frameData.systemQuadrupole = Mat3x3d(0.0);
154 frameData.convectiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
155 frameData.electronicTemperature = 0.0;
156 frameData.netCharge = 0.0;
157 frameData.chargeMomentum = 0.0;
158 frameData.COM = V3Zero;
159 frameData.COMvel = V3Zero;
160 frameData.COMw = V3Zero;
162 hasTotalEnergy =
false;
163 hasTranslationalKineticEnergy =
false;
164 hasRotationalKineticEnergy =
false;
165 hasElectronicKineticEnergy =
false;
166 hasKineticEnergy =
false;
167 hasShortRangePotential =
false;
168 hasLongRangePotential =
false;
169 hasExcludedPotential =
false;
170 hasSelfPotential =
false;
171 hasPotentialEnergy =
false;
177 hasTemperature =
false;
178 hasElectronicTemperature =
false;
179 hasNetCharge =
false;
180 hasChargeMomentum =
false;
184 hasPressureTensor =
false;
185 hasSystemDipole =
false;
186 hasSystemQuadrupole =
false;
187 hasConvectiveHeatFlux =
false;
188 hasInertiaTensor =
false;
189 hasGyrationalVolume =
false;
190 hasHullVolume =
false;
191 hasBoundingBox =
false;
200 int Snapshot::getSize() {
201 return atomData.getSize() + rigidbodyData.getSize();
223 frameData.invHmat = frameData.hmat.
inverse();
226 bool oldOrthoRhombic = frameData.orthoRhombic;
228 RealType smallDiag = fabs(frameData.hmat(0, 0));
229 if (smallDiag > fabs(frameData.hmat(1, 1)))
230 smallDiag = fabs(frameData.hmat(1, 1));
231 if (smallDiag > fabs(frameData.hmat(2, 2)))
232 smallDiag = fabs(frameData.hmat(2, 2));
233 RealType tol = smallDiag * orthoTolerance_;
235 frameData.orthoRhombic =
true;
237 for (
int i = 0; i < 3; i++) {
238 for (
int j = 0; j < 3; j++) {
240 if (frameData.orthoRhombic) {
241 if (fabs(frameData.hmat(i, j)) >= tol)
242 frameData.orthoRhombic =
false;
248 if (oldOrthoRhombic != frameData.orthoRhombic) {
289 frameData.invBbox = frameData.bBox.
inverse();
290 hasBoundingBox =
true;
296 RealType Snapshot::getXYarea() {
298 Vector3d x = frameData.hmat.getColumn(0);
299 Vector3d y = frameData.hmat.getColumn(1);
300 frameData.xyArea =
cross(x, y).length();
303 return frameData.xyArea;
306 RealType Snapshot::getXZarea() {
308 Vector3d x = frameData.hmat.getColumn(0);
309 Vector3d z = frameData.hmat.getColumn(2);
310 frameData.xzArea =
cross(x, z).length();
313 return frameData.xzArea;
316 RealType Snapshot::getYZarea() {
318 Vector3d y = frameData.hmat.getColumn(1);
319 Vector3d z = frameData.hmat.getColumn(2);
320 frameData.yzArea =
cross(y, z).length();
323 return frameData.yzArea;
326 RealType Snapshot::getVolume() {
328 frameData.volume = frameData.hmat.determinant();
331 return frameData.volume;
334 void Snapshot::setVolume(RealType vol) {
336 frameData.volume = vol;
341 if (!frameData.usePBC)
return;
343 if (!frameData.orthoRhombic) {
344 Vector3d scaled = frameData.invHmat * pos;
345 for (
int i = 0; i < 3; i++) {
346 scaled[i] -= roundMe(scaled[i]);
349 pos = frameData.hmat * scaled;
352 for (
int i = 0; i < 3; i++) {
353 scaled = pos[i] * frameData.invHmat(i, i);
354 scaled -= roundMe(scaled);
355 pos[i] = scaled * frameData.hmat(i, i);
364 if (!frameData.orthoRhombic)
365 scaled = frameData.invHmat * pos;
368 for (
int i = 0; i < 3; i++)
369 scaled[i] = pos[i] * frameData.invHmat(i, i);
375 void Snapshot::setCOM(
const Vector3d& com) {
380 void Snapshot::setCOMvel(
const Vector3d& comVel) {
381 frameData.COMvel = comVel;
385 void Snapshot::setCOMw(
const Vector3d& comw) {
386 frameData.COMw = comw;
390 Vector3d Snapshot::getCOM() {
return frameData.COM; }
392 Vector3d Snapshot::getCOMvel() {
return frameData.COMvel; }
394 Vector3d Snapshot::getCOMw() {
return frameData.COMw; }
396 RealType Snapshot::getTime() {
return frameData.currentTime; }
398 void Snapshot::increaseTime(RealType dt) { setTime(getTime() + dt); }
400 void Snapshot::setTime(RealType time) { frameData.currentTime = time; }
402 void Snapshot::setBondPotential(RealType bp) {
403 frameData.bondPotential = bp;
404 hasShortRangePotential =
false;
405 hasPotentialEnergy =
false;
406 hasTotalEnergy =
false;
409 void Snapshot::setBendPotential(RealType bp) {
410 frameData.bendPotential = bp;
411 hasShortRangePotential =
false;
412 hasPotentialEnergy =
false;
413 hasTotalEnergy =
false;
416 void Snapshot::setTorsionPotential(RealType tp) {
417 frameData.torsionPotential = tp;
418 hasShortRangePotential =
false;
419 hasPotentialEnergy =
false;
420 hasTotalEnergy =
false;
423 void Snapshot::setInversionPotential(RealType ip) {
424 frameData.inversionPotential = ip;
425 hasShortRangePotential =
false;
426 hasPotentialEnergy =
false;
427 hasTotalEnergy =
false;
430 RealType Snapshot::getBondPotential() {
return frameData.bondPotential; }
432 RealType Snapshot::getBendPotential() {
return frameData.bendPotential; }
434 RealType Snapshot::getTorsionPotential() {
435 return frameData.torsionPotential;
438 RealType Snapshot::getInversionPotential() {
439 return frameData.inversionPotential;
442 RealType Snapshot::getShortRangePotential() {
443 if (!hasShortRangePotential) {
444 frameData.shortRangePotential = frameData.bondPotential;
445 frameData.shortRangePotential += frameData.bendPotential;
446 frameData.shortRangePotential += frameData.torsionPotential;
447 frameData.shortRangePotential += frameData.inversionPotential;
448 hasShortRangePotential =
true;
449 hasPotentialEnergy =
false;
450 hasTotalEnergy =
false;
452 return frameData.shortRangePotential;
455 void Snapshot::setSurfacePotential(RealType sp) {
456 frameData.surfacePotential = sp;
457 hasLongRangePotential =
false;
458 hasPotentialEnergy =
false;
459 hasTotalEnergy =
false;
462 RealType Snapshot::getSurfacePotential() {
463 return frameData.surfacePotential;
466 void Snapshot::setReciprocalPotential(RealType rp) {
467 frameData.reciprocalPotential = rp;
468 hasLongRangePotential =
false;
469 hasPotentialEnergy =
false;
472 RealType Snapshot::getReciprocalPotential() {
473 return frameData.reciprocalPotential;
476 void Snapshot::setSelfPotentials(potVec sp) {
477 frameData.selfPotentials = sp;
478 hasSelfPotential =
false;
479 hasPotentialEnergy =
false;
480 hasTotalEnergy =
false;
483 potVec Snapshot::getSelfPotentials() {
return frameData.selfPotentials; }
485 RealType Snapshot::getSelfPotential() {
486 if (!hasSelfPotential) {
487 frameData.selfPotential = 0.0;
488 for (
int i = 0; i < N_INTERACTION_FAMILIES; i++) {
489 frameData.selfPotential += frameData.selfPotentials[i];
491 hasSelfPotential =
true;
492 hasPotentialEnergy =
false;
493 hasTotalEnergy =
false;
495 return frameData.selfPotential;
498 void Snapshot::setLongRangePotentials(potVec lrPot) {
499 frameData.lrPotentials = lrPot;
500 hasLongRangePotential =
false;
501 hasPotentialEnergy =
false;
502 hasTotalEnergy =
false;
505 RealType Snapshot::getLongRangePotential() {
506 if (!hasLongRangePotential) {
507 frameData.longRangePotential = 0.0;
508 for (
int i = 0; i < N_INTERACTION_FAMILIES; i++) {
509 frameData.longRangePotential += frameData.lrPotentials[i];
511 frameData.longRangePotential += frameData.reciprocalPotential;
512 frameData.longRangePotential += frameData.surfacePotential;
513 hasLongRangePotential =
true;
514 hasPotentialEnergy =
false;
515 hasTotalEnergy =
false;
517 return frameData.longRangePotential;
520 potVec Snapshot::getLongRangePotentials() {
return frameData.lrPotentials; }
522 RealType Snapshot::getPotentialEnergy() {
523 if (!hasPotentialEnergy) {
524 frameData.potentialEnergy = this->getLongRangePotential();
525 frameData.potentialEnergy += this->getShortRangePotential();
526 frameData.potentialEnergy += this->getSelfPotential();
527 frameData.potentialEnergy += this->getExcludedPotential();
528 hasPotentialEnergy =
true;
529 hasTotalEnergy =
false;
531 return frameData.potentialEnergy;
534 void Snapshot::setPotentialEnergy(
const RealType pe) {
535 frameData.potentialEnergy = pe;
536 hasPotentialEnergy =
true;
537 hasTotalEnergy =
false;
540 void Snapshot::setExcludedPotentials(potVec exPot) {
541 frameData.excludedPotentials = exPot;
542 hasExcludedPotential =
false;
543 hasPotentialEnergy =
false;
544 hasTotalEnergy =
false;
547 potVec Snapshot::getExcludedPotentials() {
548 return frameData.excludedPotentials;
551 RealType Snapshot::getExcludedPotential() {
552 if (!hasExcludedPotential) {
553 frameData.excludedPotential = 0.0;
554 for (
int i = 0; i < N_INTERACTION_FAMILIES; i++) {
555 frameData.excludedPotential += frameData.excludedPotentials[i];
557 hasExcludedPotential =
true;
558 hasPotentialEnergy =
false;
559 hasTotalEnergy =
false;
561 return frameData.excludedPotential;
564 void Snapshot::setRestraintPotential(RealType rp) {
565 frameData.restraintPotential = rp;
568 RealType Snapshot::getRestraintPotential() {
569 return frameData.restraintPotential;
572 void Snapshot::setRawPotential(RealType rp) { frameData.rawPotential = rp; }
574 RealType Snapshot::getRawPotential() {
return frameData.rawPotential; }
576 void Snapshot::setSelectionPotentials(potVec selPot) {
577 frameData.selectionPotentials = selPot;
580 potVec Snapshot::getSelectionPotentials() {
581 return frameData.selectionPotentials;
584 RealType Snapshot::getTranslationalKineticEnergy() {
585 return frameData.translationalKinetic;
588 RealType Snapshot::getRotationalKineticEnergy() {
589 return frameData.rotationalKinetic;
592 RealType Snapshot::getElectronicKineticEnergy() {
593 return frameData.electronicKinetic;
596 RealType Snapshot::getKineticEnergy() {
return frameData.kineticEnergy; }
598 void Snapshot::setTranslationalKineticEnergy(RealType tke) {
599 frameData.translationalKinetic = tke;
600 hasTranslationalKineticEnergy =
true;
601 hasKineticEnergy =
false;
602 hasTotalEnergy =
false;
605 void Snapshot::setRotationalKineticEnergy(RealType rke) {
606 frameData.rotationalKinetic = rke;
607 hasRotationalKineticEnergy =
true;
608 hasKineticEnergy =
false;
609 hasTotalEnergy =
false;
612 void Snapshot::setElectronicKineticEnergy(RealType eke) {
613 frameData.electronicKinetic = eke;
614 hasElectronicKineticEnergy =
true;
615 hasKineticEnergy =
false;
616 hasTotalEnergy =
false;
619 void Snapshot::setKineticEnergy(RealType ke) {
620 frameData.kineticEnergy = ke;
621 hasKineticEnergy =
true;
622 hasTotalEnergy =
false;
625 RealType Snapshot::getTotalEnergy() {
return frameData.totalEnergy; }
627 void Snapshot::setTotalEnergy(RealType te) {
628 frameData.totalEnergy = te;
629 hasTotalEnergy =
true;
632 RealType Snapshot::getConservedQuantity() {
633 return frameData.conservedQuantity;
636 void Snapshot::setConservedQuantity(RealType cq) {
637 frameData.conservedQuantity = cq;
640 RealType Snapshot::getTemperature() {
return frameData.temperature; }
642 void Snapshot::setTemperature(RealType temp) {
643 hasTemperature =
true;
644 frameData.temperature = temp;
647 RealType Snapshot::getElectronicTemperature() {
648 return frameData.electronicTemperature;
651 void Snapshot::setElectronicTemperature(RealType eTemp) {
652 hasElectronicTemperature =
true;
653 frameData.electronicTemperature = eTemp;
656 RealType Snapshot::getNetCharge() {
return frameData.netCharge; }
658 void Snapshot::setNetCharge(RealType nChg) {
660 frameData.netCharge = nChg;
663 RealType Snapshot::getChargeMomentum() {
return frameData.chargeMomentum; }
665 void Snapshot::setChargeMomentum(RealType cMom) {
666 hasChargeMomentum =
true;
667 frameData.chargeMomentum = cMom;
670 RealType Snapshot::getPressure() {
return frameData.pressure; }
672 void Snapshot::setPressure(RealType pressure) {
674 frameData.pressure = pressure;
677 Mat3x3d Snapshot::getPressureTensor() {
return frameData.pressureTensor; }
679 void Snapshot::setPressureTensor(
const Mat3x3d& pressureTensor) {
680 hasPressureTensor =
true;
681 frameData.pressureTensor = pressureTensor;
684 void Snapshot::setVirialTensor(
const Mat3x3d& virialTensor) {
685 frameData.virialTensor = virialTensor;
688 Mat3x3d Snapshot::getVirialTensor() {
return frameData.virialTensor; }
690 void Snapshot::setConductiveHeatFlux(
const Vector3d& chf) {
691 frameData.conductiveHeatFlux = chf;
694 Vector3d Snapshot::getConductiveHeatFlux() {
695 return frameData.conductiveHeatFlux;
698 Vector3d Snapshot::getConvectiveHeatFlux() {
699 return frameData.convectiveHeatFlux;
702 void Snapshot::setConvectiveHeatFlux(
const Vector3d& chf) {
703 hasConvectiveHeatFlux =
true;
704 frameData.convectiveHeatFlux = chf;
707 Vector3d Snapshot::getHeatFlux() {
709 return getConductiveHeatFlux() + getConvectiveHeatFlux();
712 Vector3d Snapshot::getSystemDipole() {
return frameData.systemDipole; }
714 void Snapshot::setSystemDipole(
const Vector3d& bd) {
715 hasSystemDipole =
true;
716 frameData.systemDipole = bd;
719 Mat3x3d Snapshot::getSystemQuadrupole() {
return frameData.systemQuadrupole; }
721 void Snapshot::setSystemQuadrupole(
const Mat3x3d& bq) {
722 hasSystemQuadrupole =
true;
723 frameData.systemQuadrupole = bq;
726 void Snapshot::setThermostat(
const pair<RealType, RealType>& thermostat) {
727 frameData.thermostat = thermostat;
730 pair<RealType, RealType> Snapshot::getThermostat() {
731 return frameData.thermostat;
734 void Snapshot::setElectronicThermostat(
735 const pair<RealType, RealType>& eTherm) {
736 frameData.electronicThermostat = eTherm;
739 pair<RealType, RealType> Snapshot::getElectronicThermostat() {
740 return frameData.electronicThermostat;
743 void Snapshot::setBarostat(
const Mat3x3d& barostat) {
744 frameData.barostat = barostat;
747 Mat3x3d Snapshot::getBarostat() {
return frameData.barostat; }
749 void Snapshot::setSPFData(std::shared_ptr<SPFData> data) {
750 frameData.spfData = data;
753 std::shared_ptr<SPFData> Snapshot::getSPFData() {
return frameData.spfData; }
755 void Snapshot::setInertiaTensor(
const Mat3x3d& inertiaTensor) {
756 frameData.inertiaTensor = inertiaTensor;
757 hasInertiaTensor =
true;
760 Mat3x3d Snapshot::getInertiaTensor() {
return frameData.inertiaTensor; }
762 void Snapshot::setGyrationalVolume(
const RealType gyrationalVolume) {
763 frameData.gyrationalVolume = gyrationalVolume;
764 hasGyrationalVolume =
true;
767 RealType Snapshot::getGyrationalVolume() {
768 return frameData.gyrationalVolume;
771 void Snapshot::setHullVolume(
const RealType hullVolume) {
772 frameData.hullVolume = hullVolume;
773 hasHullVolume =
true;
776 RealType Snapshot::getHullVolume() {
return frameData.hullVolume; }
778 void Snapshot::setOrthoTolerance(RealType ot) { orthoTolerance_ = ot; }
Vector3d scaleVector(Vector3d &v)
Scaling a vector to multiples of the periodic box.
Mat3x3d getHmat()
Returns the H-Matrix.
Mat3x3d getInvHmat()
Returns the inverse H-Matrix.
void clearDerivedProperties()
sets the state of the computed properties to false
void setHmat(const Mat3x3d &m)
Sets the H-Matrix.
Mat3x3d getBoundingBox()
Returns the Bounding Box.
int getNumberOfAtoms()
Returns the number of atoms.
int getNumberOfRigidBodies()
Returns the number of rigid bodies.
void setID(int id)
Sets the id of this Snapshot.
Mat3x3d getInvBoundingBox()
Returns the inverse Bounding Box.
void setBoundingBox(const Mat3x3d &m)
Sets the Bounding Box.
static int getFrameDataSize()
Returns the number of bytes in a FrameData structure.
void wrapVector(Vector3d &v)
Wrapping the vector according to periodic boundary condition.
int getNumberOfCutoffGroups()
Returns the number of rigid bodies.
int getID()
Returns the id of this Snapshot.
SquareMatrix3< Real > inverse() const
Sets the value of this matrix to the inverse of itself.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
FrameData is a structure for holding system-wide dynamic data about the simulation.