| 43 |
|
#include <algorithm> |
| 44 |
|
#include <fstream> |
| 45 |
|
#include "applications/staticProps/RNEMDStats.hpp" |
| 46 |
+ |
#include "primitives/Molecule.hpp" |
| 47 |
|
#include "utils/PhysicalConstants.hpp" |
| 48 |
|
|
| 49 |
|
namespace OpenMD { |
| 85 |
|
data_.push_back(density); |
| 86 |
|
} |
| 87 |
|
|
| 88 |
< |
void RNEMDZ::processStuntDouble(StuntDouble* sd, int bin) { |
| 89 |
< |
RealType mass = sd->getMass(); |
| 90 |
< |
Vector3d pos = sd->getPos(); |
| 91 |
< |
Vector3d vel = sd->getVel(); |
| 92 |
< |
RealType KE = 0.5 * (mass * vel.lengthSquare()); |
| 93 |
< |
int dof = 3; |
| 88 |
> |
void RNEMDZ::processFrame(int istep) { |
| 89 |
> |
Molecule* mol; |
| 90 |
> |
RigidBody* rb; |
| 91 |
> |
StuntDouble* sd; |
| 92 |
> |
SimInfo::MoleculeIterator mi; |
| 93 |
> |
Molecule::RigidBodyIterator rbIter; |
| 94 |
> |
int i; |
| 95 |
|
|
| 96 |
< |
if (sd->isDirectional()) { |
| 97 |
< |
Vector3d angMom = sd->getJ(); |
| 98 |
< |
Mat3x3d I = sd->getI(); |
| 99 |
< |
if (sd->isLinear()) { |
| 100 |
< |
int i = sd->linearAxis(); |
| 101 |
< |
int j = (i + 1) % 3; |
| 102 |
< |
int k = (i + 2) % 3; |
| 103 |
< |
KE += 0.5 * (angMom[j] * angMom[j] / I(j, j) + |
| 104 |
< |
angMom[k] * angMom[k] / I(k, k)); |
| 105 |
< |
dof += 2; |
| 106 |
< |
} else { |
| 107 |
< |
KE += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + |
| 108 |
< |
angMom[1] * angMom[1] / I(1, 1) + |
| 109 |
< |
angMom[2] * angMom[2] / I(2, 2)); |
| 110 |
< |
dof += 3; |
| 96 |
> |
vector<RealType> binMass(nBins_, 0.0); |
| 97 |
> |
vector<Vector3d> binVel(nBins_, V3Zero); |
| 98 |
> |
vector<RealType> binKE(nBins_, 0.0); |
| 99 |
> |
vector<int> binDof(nBins_, 0); |
| 100 |
> |
vector<int> binCount(nBins_, 0); |
| 101 |
> |
|
| 102 |
> |
|
| 103 |
> |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 104 |
> |
mol = info_->nextMolecule(mi)) { |
| 105 |
> |
|
| 106 |
> |
// change the positions of atoms which belong to the rigidbodies |
| 107 |
> |
|
| 108 |
> |
for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
| 109 |
> |
rb = mol->nextRigidBody(rbIter)) { |
| 110 |
> |
rb->updateAtoms(); |
| 111 |
|
} |
| 112 |
|
} |
| 113 |
|
|
| 114 |
< |
RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb * |
| 115 |
< |
PhysicalConstants::energyConvert); |
| 116 |
< |
RealType den = mass * nBins_ * PhysicalConstants::densityConvert / volume_; |
| 114 |
> |
if (evaluator_.isDynamic()) { |
| 115 |
> |
seleMan_.setSelectionSet(evaluator_.evaluate()); |
| 116 |
> |
} |
| 117 |
|
|
| 118 |
< |
dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp); |
| 119 |
< |
dynamic_cast<VectorAccumulator *>(velocity->accumulator[bin])->add(vel); |
| 120 |
< |
dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den); |
| 118 |
> |
// loop over the selected atoms: |
| 119 |
> |
|
| 120 |
> |
for (sd = seleMan_.beginSelected(i); sd != NULL; |
| 121 |
> |
sd = seleMan_.nextSelected(i)) { |
| 122 |
> |
|
| 123 |
> |
// figure out where that object is: |
| 124 |
> |
Vector3d pos = sd->getPos(); |
| 125 |
> |
currentSnapshot_->wrapVector(pos); |
| 126 |
|
|
| 127 |
+ |
int bin = getBin(pos); |
| 128 |
+ |
binCount[bin]++; |
| 129 |
+ |
|
| 130 |
+ |
RealType m = sd->getMass(); |
| 131 |
+ |
binMass[bin] += m; |
| 132 |
+ |
Vector3d vel = sd->getVel(); |
| 133 |
+ |
binVel[bin] += vel; |
| 134 |
+ |
binKE[bin] += 0.5 * (m * vel.lengthSquare()); |
| 135 |
+ |
binDof[bin] += 3; |
| 136 |
+ |
|
| 137 |
+ |
if (sd->isDirectional()) { |
| 138 |
+ |
Vector3d angMom = sd->getJ(); |
| 139 |
+ |
Mat3x3d I = sd->getI(); |
| 140 |
+ |
if (sd->isLinear()) { |
| 141 |
+ |
int i = sd->linearAxis(); |
| 142 |
+ |
int j = (i + 1) % 3; |
| 143 |
+ |
int k = (i + 2) % 3; |
| 144 |
+ |
binKE[bin] += 0.5 * (angMom[j] * angMom[j] / I(j, j) + |
| 145 |
+ |
angMom[k] * angMom[k] / I(k, k)); |
| 146 |
+ |
binDof[bin] += 2; |
| 147 |
+ |
} else { |
| 148 |
+ |
binKE[bin] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + |
| 149 |
+ |
angMom[1] * angMom[1] / I(1, 1) + |
| 150 |
+ |
angMom[2] * angMom[2] / I(2, 2)); |
| 151 |
+ |
binDof[bin] += 3; |
| 152 |
+ |
} |
| 153 |
+ |
} |
| 154 |
+ |
} |
| 155 |
+ |
|
| 156 |
+ |
for (int i = 0; i < nBins_; i++) { |
| 157 |
+ |
RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb * |
| 158 |
+ |
PhysicalConstants::energyConvert); |
| 159 |
+ |
RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert |
| 160 |
+ |
/ volume_; |
| 161 |
+ |
Vector3d vel = binVel[i] / RealType(binCount[i]); |
| 162 |
+ |
dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp); |
| 163 |
+ |
dynamic_cast<VectorAccumulator *>(velocity->accumulator[i])->add(vel); |
| 164 |
+ |
dynamic_cast<Accumulator *>(density->accumulator[i])->add(den); |
| 165 |
+ |
dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1); |
| 166 |
+ |
} |
| 167 |
|
} |
| 168 |
+ |
|
| 169 |
+ |
void RNEMDZ::processStuntDouble(StuntDouble* sd, int bin) { |
| 170 |
+ |
} |
| 171 |
|
|
| 172 |
|
RNEMDR::RNEMDR(SimInfo* info, const std::string& filename, |
| 173 |
|
const std::string& sele, int nrbins) |
