ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/RNEMDStats.cpp
(Generate patch)

Comparing:
branches/development/src/applications/staticProps/RNEMDStats.cpp (file contents), Revision 1865 by gezelter, Wed Apr 17 18:24:08 2013 UTC vs.
trunk/src/applications/staticProps/RNEMDStats.cpp (file contents), Revision 1887 by gezelter, Tue Jun 18 17:48:52 2013 UTC

# Line 43 | Line 43
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 {
# Line 84 | Line 85 | namespace OpenMD {
85      data_.push_back(density);
86    }
87  
88 <  void RNEMDZ::processStuntDouble(StuntDouble* sd, int bin) {
89 <    RealType mass = sd->getMass();
89 <    Vector3d pos = sd->getPos();    
90 <    Vector3d vel = sd->getVel();
91 <    RealType KE = 0.5 * (mass * vel.lengthSquare());
92 <    int dof = 3;
88 >  void RNEMDZ::processFrame(int istep) {
89 >    RealType z;
90  
91 <    if (sd->isDirectional()) {
92 <      Vector3d angMom = sd->getJ();
93 <      Mat3x3d I = sd->getI();
94 <      if (sd->isLinear()) {
95 <        int i = sd->linearAxis();
96 <        int j = (i + 1) % 3;
97 <        int k = (i + 2) % 3;
98 <        KE += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
99 <                     angMom[k] * angMom[k] / I(k, k));
100 <        dof += 2;
101 <      } else {
102 <        KE += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
103 <                     angMom[1] * angMom[1] / I(1, 1) +
104 <                     angMom[2] * angMom[2] / I(2, 2));
105 <        dof += 3;
91 >    hmat_ = currentSnapshot_->getHmat();
92 >    for (int i = 0; i < nBins_; i++) {
93 >      z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat_(2,2);
94 >      dynamic_cast<Accumulator*>(z_->accumulator[i])->add(z);
95 >    }
96 >    volume_ = currentSnapshot_->getVolume();
97 >
98 >
99 >    Molecule* mol;
100 >    RigidBody* rb;
101 >    StuntDouble* sd;
102 >    SimInfo::MoleculeIterator mi;
103 >    Molecule::RigidBodyIterator rbIter;
104 >    int i;
105 >
106 >    vector<RealType> binMass(nBins_, 0.0);
107 >    vector<Vector3d> binVel(nBins_, V3Zero);
108 >    vector<RealType> binKE(nBins_, 0.0);
109 >    vector<unsigned int> binDof(nBins_, 0);
110 >    vector<unsigned int> binCount(nBins_, 0);
111 >    
112 >    
113 >    for (mol = info_->beginMolecule(mi); mol != NULL;
114 >         mol = info_->nextMolecule(mi)) {
115 >      
116 >      // change the positions of atoms which belong to the rigidbodies
117 >      
118 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
119 >           rb = mol->nextRigidBody(rbIter)) {
120 >        rb->updateAtoms();
121        }
122      }
123 +  
124 +    if (evaluator_.isDynamic()) {
125 +      seleMan_.setSelectionSet(evaluator_.evaluate());
126 +    }
127      
128 <    RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb *
113 <                                PhysicalConstants::energyConvert);
114 <    RealType den = mass * nBins_ * PhysicalConstants::densityConvert / volume_;
128 >    // loop over the selected atoms:
129      
130 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
131 <    dynamic_cast<VectorAccumulator *>(velocity->accumulator[bin])->add(vel);
132 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
130 >    for (sd = seleMan_.beginSelected(i); sd != NULL;
131 >         sd = seleMan_.nextSelected(i)) {
132 >      
133 >      // figure out where that object is:
134 >      Vector3d pos = sd->getPos();
135 >      Vector3d vel = sd->getVel();
136 >      RealType m = sd->getMass();
137  
138 +      int bin = getBin(pos);
139 +
140 +      binCount[bin] += 1;
141 +
142 +      binMass[bin] += m;
143 +      binVel[bin] += vel;
144 +      binKE[bin] += 0.5 * (m * vel.lengthSquare());
145 +      binDof[bin] += 3;
146 +      
147 +      if (sd->isDirectional()) {
148 +        Vector3d angMom = sd->getJ();
149 +        Mat3x3d I = sd->getI();
150 +        if (sd->isLinear()) {
151 +          int i = sd->linearAxis();
152 +          int j = (i + 1) % 3;
153 +          int k = (i + 2) % 3;
154 +          binKE[bin] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
155 +                               angMom[k] * angMom[k] / I(k, k));
156 +          binDof[bin] += 2;
157 +        } else {
158 +          binKE[bin] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
159 +                               angMom[1] * angMom[1] / I(1, 1) +
160 +                               angMom[2] * angMom[2] / I(2, 2));
161 +          binDof[bin] += 3;
162 +        }
163 +      }
164 +    }
165 +    
166 +    for (unsigned int i = 0; i < nBins_; i++) {
167 +
168 +      if (binDof[i] > 0) {
169 +        RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
170 +                                          PhysicalConstants::energyConvert);
171 +        RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
172 +          / volume_;
173 +        Vector3d vel = binVel[i] / RealType(binCount[i]);
174 +        dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
175 +        dynamic_cast<VectorAccumulator *>(velocity->accumulator[i])->add(vel);
176 +        dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
177 +        dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
178 +      }
179 +    }
180    }
181 +  
182 +  void RNEMDZ::processStuntDouble(StuntDouble* sd, int bin) {
183 +  }
184  
185    RNEMDR::RNEMDR(SimInfo* info, const std::string& filename,
186                   const std::string& sele, int nrbins)
# Line 157 | Line 220 | namespace OpenMD {
220      data_.push_back(density);
221    }
222  
223 <  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
224 <    RealType mass = sd->getMass();
162 <    Vector3d vel = sd->getVel();
163 <    Vector3d rPos = sd->getPos() - coordinateOrigin_;
164 <    Vector3d aVel = cross(rPos, vel);
223 >
224 >  void RNEMDR::processFrame(int istep) {
225  
226 <    RealType KE = 0.5 * (mass * vel.lengthSquare());
227 <    int dof = 3;
226 >    Molecule* mol;
227 >    RigidBody* rb;
228 >    StuntDouble* sd;
229 >    SimInfo::MoleculeIterator mi;
230 >    Molecule::RigidBodyIterator rbIter;
231 >    int i;
232  
233 <    if (sd->isDirectional()) {
234 <      Vector3d angMom = sd->getJ();
235 <      Mat3x3d I = sd->getI();
236 <      if (sd->isLinear()) {
237 <        int i = sd->linearAxis();
238 <        int j = (i + 1) % 3;
239 <        int k = (i + 2) % 3;
240 <        KE += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
241 <                     angMom[k] * angMom[k] / I(k, k));
242 <        dof += 2;
243 <      } else {
244 <        KE += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
245 <                     angMom[1] * angMom[1] / I(1, 1) +
246 <                     angMom[2] * angMom[2] / I(2, 2));
183 <        dof += 3;
233 >    vector<RealType> binMass(nBins_, 0.0);
234 >    vector<Vector3d> binaVel(nBins_, V3Zero);
235 >    vector<RealType> binKE(nBins_, 0.0);
236 >    vector<unsigned int> binDof(nBins_, 0);
237 >    vector<unsigned int> binCount(nBins_, 0);
238 >    
239 >    for (mol = info_->beginMolecule(mi); mol != NULL;
240 >         mol = info_->nextMolecule(mi)) {
241 >      
242 >      // change the positions of atoms which belong to the rigidbodies
243 >      
244 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
245 >           rb = mol->nextRigidBody(rbIter)) {
246 >        rb->updateAtoms();
247        }
248      }
249 +  
250 +    if (evaluator_.isDynamic()) {
251 +      seleMan_.setSelectionSet(evaluator_.evaluate());
252 +    }
253      
254 <    RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb *
255 <                                PhysicalConstants::energyConvert);
254 >    // loop over the selected atoms:
255 >    
256 >    for (sd = seleMan_.beginSelected(i); sd != NULL;
257 >         sd = seleMan_.nextSelected(i)) {
258 >      
259 >      // figure out where that object is:
260  
261 <    RealType rinner = (RealType)bin * binWidth_;
262 <    RealType router = (RealType)(bin+1) * binWidth_;
263 <    RealType den = mass * 3.0 * PhysicalConstants::densityConvert
264 <      / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));  
261 >      Vector3d rPos = sd->getPos() - coordinateOrigin_;
262 >      Vector3d vel = sd->getVel();      
263 >      Vector3d aVel = cross(rPos, vel);
264 >      RealType m = sd->getMass();
265 >
266 >      int bin = getBin(pos);
267 >
268 >      binCount[bin] += 1;
269 >
270 >      binMass[bin] += m;
271 >      binaVel[bin] += aVel;
272 >      binKE[bin] += 0.5 * (m * vel.lengthSquare());
273 >      binDof[bin] += 3;
274 >      
275 >      if (sd->isDirectional()) {
276 >        Vector3d angMom = sd->getJ();
277 >        Mat3x3d I = sd->getI();
278 >        if (sd->isLinear()) {
279 >          int i = sd->linearAxis();
280 >          int j = (i + 1) % 3;
281 >          int k = (i + 2) % 3;
282 >          binKE[bin] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
283 >                               angMom[k] * angMom[k] / I(k, k));
284 >          binDof[bin] += 2;
285 >        } else {
286 >          binKE[bin] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
287 >                               angMom[1] * angMom[1] / I(1, 1) +
288 >                               angMom[2] * angMom[2] / I(2, 2));
289 >          binDof[bin] += 3;
290 >        }
291 >      }
292 >    }
293      
294 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
295 <    dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[bin])->add(aVel);
296 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
294 >    for (unsigned int i = 0; i < nBins_; i++) {
295 >      RealType rinner = (RealType)i * binWidth_;
296 >      RealType router = (RealType)(i+1) * binWidth_;
297 >      if (binDof[i] > 0) {
298 >        RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
299 >                                          PhysicalConstants::energyConvert);
300 >        RealType den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
301 >          / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));  
302 >        Vector3d aVel = binaVel[i] / RealType(binCount[i]);
303 >        dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
304 >        dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[i])->add(aVel);
305 >        dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
306 >        dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
307 >      }
308 >    }
309 >  }
310  
311 +
312 +  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
313    }
314   }
315  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines