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

Comparing trunk/src/applications/staticProps/RNEMDStats.cpp (file contents):
Revision 1882 by gezelter, Tue Jun 18 16:10:07 2013 UTC vs.
Revision 1888 by gezelter, Tue Jun 18 17:52:37 2013 UTC

# Line 86 | Line 86 | namespace OpenMD {
86    }
87  
88    void RNEMDZ::processFrame(int istep) {
89 +    RealType z;
90 +
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;
# Line 96 | Line 106 | namespace OpenMD {
106      vector<RealType> binMass(nBins_, 0.0);
107      vector<Vector3d> binVel(nBins_, V3Zero);
108      vector<RealType> binKE(nBins_, 0.0);
109 <    vector<int> binDof(nBins_, 0);
110 <    vector<int> binCount(nBins_, 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;
# Line 110 | Line 120 | namespace OpenMD {
120          rb->updateAtoms();
121        }
122      }
123 <    
123 >  
124      if (evaluator_.isDynamic()) {
125        seleMan_.setSelectionSet(evaluator_.evaluate());
126      }
# Line 122 | Line 132 | namespace OpenMD {
132        
133        // figure out where that object is:
134        Vector3d pos = sd->getPos();
135 <      currentSnapshot_->wrapVector(pos);
135 >      Vector3d vel = sd->getVel();
136 >      RealType m = sd->getMass();
137  
138        int bin = getBin(pos);
128 –      binCount[bin]++;
139  
140 <      RealType m = sd->getMass();
140 >      binCount[bin] += 1;
141 >
142        binMass[bin] += m;
132 –      Vector3d vel = sd->getVel();
143        binVel[bin] += vel;
144        binKE[bin] += 0.5 * (m * vel.lengthSquare());
145        binDof[bin] += 3;
# Line 153 | Line 163 | namespace OpenMD {
163        }
164      }
165      
166 <    for (int i = 0; i < nBins_; i++) {
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);
# Line 209 | Line 220 | namespace OpenMD {
220      data_.push_back(density);
221    }
222  
212 –  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
213 –    RealType mass = sd->getMass();
214 –    Vector3d vel = sd->getVel();
215 –    Vector3d rPos = sd->getPos() - coordinateOrigin_;
216 –    Vector3d aVel = cross(rPos, vel);
223  
224 <    RealType KE = 0.5 * (mass * vel.lengthSquare());
219 <    int dof = 3;
224 >  void RNEMDR::processFrame(int istep) {
225  
226 <    if (sd->isDirectional()) {
227 <      Vector3d angMom = sd->getJ();
228 <      Mat3x3d I = sd->getI();
229 <      if (sd->isLinear()) {
230 <        int i = sd->linearAxis();
231 <        int j = (i + 1) % 3;
227 <        int k = (i + 2) % 3;
228 <        KE += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
229 <                     angMom[k] * angMom[k] / I(k, k));
230 <        dof += 2;
231 <      } else {
232 <        KE += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
233 <                     angMom[1] * angMom[1] / I(1, 1) +
234 <                     angMom[2] * angMom[2] / I(2, 2));
235 <        dof += 3;
236 <      }
237 <    }
238 <    
239 <    RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb *
240 <                                PhysicalConstants::energyConvert);
226 >    Molecule* mol;
227 >    RigidBody* rb;
228 >    StuntDouble* sd;
229 >    SimInfo::MoleculeIterator mi;
230 >    Molecule::RigidBodyIterator rbIter;
231 >    int i;
232  
233 <    RealType rinner = (RealType)bin * binWidth_;
234 <    RealType router = (RealType)(bin+1) * binWidth_;
235 <    RealType den = mass * 3.0 * PhysicalConstants::densityConvert
236 <      / (4.0 * M_PI * (pow(router,3) - pow(rinner,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 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
240 <    dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[bin])->add(aVel);
241 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
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 >    // 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 +      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(rPos);
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 +    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