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 1881 by gezelter, Tue Jun 18 16:07:27 2013 UTC vs.
Revision 1892 by gezelter, Wed Jun 19 14:14:56 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 94 | Line 104 | namespace OpenMD {
104      int i;
105  
106      vector<RealType> binMass(nBins_, 0.0);
107 <    vector<Vector3d> binVel(nBins_, V3Zero);
107 >    vector<RealType> binPx(nBins_, 0.0);
108 >    vector<RealType> binPy(nBins_, 0.0);
109 >    vector<RealType> binPz(nBins_, 0.0);
110      vector<RealType> binKE(nBins_, 0.0);
111 <    vector<int> binDof(nBins_, 0);
112 <    vector<int> binCount(nBins_, 0);
111 >    vector<unsigned int> binDof(nBins_, 0);
112 >    vector<unsigned int> binCount(nBins_, 0);
113      
114      
115      for (mol = info_->beginMolecule(mi); mol != NULL;
# Line 107 | Line 119 | namespace OpenMD {
119        
120        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
121             rb = mol->nextRigidBody(rbIter)) {
122 <        rb->updateAtoms();
122 >        rb->updateAtomVel();
123        }
124      }
125 <    
125 >  
126      if (evaluator_.isDynamic()) {
127        seleMan_.setSelectionSet(evaluator_.evaluate());
128      }
# Line 122 | Line 134 | namespace OpenMD {
134        
135        // figure out where that object is:
136        Vector3d pos = sd->getPos();
137 <      currentSnapshot_->wrapVector(pos);
137 >      Vector3d vel = sd->getVel();
138 >      RealType m = sd->getMass();
139  
140        int bin = getBin(pos);
128 –      binCount[bin]++;
141  
142 <      RealType m = sd->getMass();
142 >      binCount[bin] += 1;
143 >
144        binMass[bin] += m;
145 <      Vector3d vel = sd->getVel();
146 <      binVel[bin] += vel;
145 >      binPx[bin] += m * vel.x();
146 >      binPy[bin] += m * vel.y();
147 >      binPz[bin] += m * vel.z();
148        binKE[bin] += 0.5 * (m * vel.lengthSquare());
149        binDof[bin] += 3;
150        
# Line 153 | Line 167 | namespace OpenMD {
167        }
168      }
169      
170 <    for (int i = 0; i < nBins_; i++) {
171 <      RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
172 <                                        PhysicalConstants::energyConvert);
173 <      RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
174 <        / volume_;
175 <      Vector3d vel = binVel[i] / RealType(binCount[i]);
176 <      dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
177 <      dynamic_cast<VectorAccumulator *>(velocity->accumulator[i])->add(vel);
178 <      dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
179 <      dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
170 >    for (unsigned int i = 0; i < nBins_; i++) {
171 >
172 >      if (binDof[i] > 0) {
173 >        RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
174 >                                          PhysicalConstants::energyConvert);
175 >        RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
176 >          / volume_;
177 >        Vector3d vel;
178 >        vel.x() = binPx[i] / binMass[i];
179 >        vel.y() = binPy[i] / binMass[i];
180 >        vel.z() = binPz[i] / binMass[i];
181 >
182 >        dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
183 >        dynamic_cast<VectorAccumulator *>(velocity->accumulator[i])->add(vel);
184 >        dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
185 >        dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
186 >      }
187      }
188    }
189    
# Line 207 | Line 228 | namespace OpenMD {
228      data_.push_back(density);
229    }
230  
210 –  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
211 –    RealType mass = sd->getMass();
212 –    Vector3d vel = sd->getVel();
213 –    Vector3d rPos = sd->getPos() - coordinateOrigin_;
214 –    Vector3d aVel = cross(rPos, vel);
231  
232 <    RealType KE = 0.5 * (mass * vel.lengthSquare());
217 <    int dof = 3;
232 >  void RNEMDR::processFrame(int istep) {
233  
234 <    if (sd->isDirectional()) {
235 <      Vector3d angMom = sd->getJ();
236 <      Mat3x3d I = sd->getI();
237 <      if (sd->isLinear()) {
238 <        int i = sd->linearAxis();
239 <        int j = (i + 1) % 3;
240 <        int k = (i + 2) % 3;
241 <        KE += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
242 <                     angMom[k] * angMom[k] / I(k, k));
243 <        dof += 2;
244 <      } else {
245 <        KE += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
246 <                     angMom[1] * angMom[1] / I(1, 1) +
247 <                     angMom[2] * angMom[2] / I(2, 2));
248 <        dof += 3;
234 >    Molecule* mol;
235 >    RigidBody* rb;
236 >    StuntDouble* sd;
237 >    SimInfo::MoleculeIterator mi;
238 >    Molecule::RigidBodyIterator rbIter;
239 >    int i;
240 >
241 >    vector<RealType> binMass(nBins_, 0.0);
242 >    vector<Vector3d> binaVel(nBins_, V3Zero);
243 >    vector<RealType> binKE(nBins_, 0.0);
244 >    vector<unsigned int> binDof(nBins_, 0);
245 >    vector<unsigned int> binCount(nBins_, 0);
246 >    
247 >    for (mol = info_->beginMolecule(mi); mol != NULL;
248 >         mol = info_->nextMolecule(mi)) {
249 >      
250 >      // change the positions of atoms which belong to the rigidbodies
251 >      
252 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
253 >           rb = mol->nextRigidBody(rbIter)) {
254 >        rb->updateAtomVel();
255        }
256      }
257 +  
258 +    if (evaluator_.isDynamic()) {
259 +      seleMan_.setSelectionSet(evaluator_.evaluate());
260 +    }
261      
262 <    RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb *
263 <                                PhysicalConstants::energyConvert);
262 >    // loop over the selected atoms:
263 >    
264 >    for (sd = seleMan_.beginSelected(i); sd != NULL;
265 >         sd = seleMan_.nextSelected(i)) {
266 >      
267 >      // figure out where that object is:
268  
269 <    RealType rinner = (RealType)bin * binWidth_;
270 <    RealType router = (RealType)(bin+1) * binWidth_;
271 <    RealType den = mass * 3.0 * PhysicalConstants::densityConvert
272 <      / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));  
269 >      Vector3d rPos = sd->getPos() - coordinateOrigin_;
270 >      Vector3d vel = sd->getVel();      
271 >      Vector3d aVel = cross(rPos, vel);
272 >      RealType m = sd->getMass();
273 >
274 >      int bin = getBin(rPos);
275 >
276 >      binCount[bin] += 1;
277 >
278 >      binMass[bin] += m;
279 >      binaVel[bin] += aVel;
280 >      binKE[bin] += 0.5 * (m * vel.lengthSquare());
281 >      binDof[bin] += 3;
282 >      
283 >      if (sd->isDirectional()) {
284 >        Vector3d angMom = sd->getJ();
285 >        Mat3x3d I = sd->getI();
286 >        if (sd->isLinear()) {
287 >          int i = sd->linearAxis();
288 >          int j = (i + 1) % 3;
289 >          int k = (i + 2) % 3;
290 >          binKE[bin] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
291 >                               angMom[k] * angMom[k] / I(k, k));
292 >          binDof[bin] += 2;
293 >        } else {
294 >          binKE[bin] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
295 >                               angMom[1] * angMom[1] / I(1, 1) +
296 >                               angMom[2] * angMom[2] / I(2, 2));
297 >          binDof[bin] += 3;
298 >        }
299 >      }
300 >    }
301      
302 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
303 <    dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[bin])->add(aVel);
304 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
302 >    for (unsigned int i = 0; i < nBins_; i++) {
303 >      RealType rinner = (RealType)i * binWidth_;
304 >      RealType router = (RealType)(i+1) * binWidth_;
305 >      if (binDof[i] > 0) {
306 >        RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
307 >                                          PhysicalConstants::energyConvert);
308 >        RealType den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
309 >          / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));  
310 >        Vector3d aVel = binaVel[i] / RealType(binCount[i]);
311 >        dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
312 >        dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[i])->add(aVel);
313 >        dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
314 >        dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
315 >      }
316 >    }
317 >  }
318  
319 +
320 +  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
321    }
322   }
323  

Diff Legend

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