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 1944 by gezelter, Wed Nov 6 18:30:25 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<Vector3d> binP(nBins_, V3Zero);
108      vector<RealType> binKE(nBins_, 0.0);
109 <    vector<int> binDof(nBins_, 0);
100 <    vector<int> binCount(nBins_, 0);
109 >    vector<unsigned int> binDof(nBins_, 0);
110      
102 –    
111      for (mol = info_->beginMolecule(mi); mol != NULL;
112           mol = info_->nextMolecule(mi)) {
113        
# Line 107 | Line 115 | namespace OpenMD {
115        
116        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
117             rb = mol->nextRigidBody(rbIter)) {
118 <        rb->updateAtoms();
118 >        rb->updateAtomVel();
119        }
120      }
121 <    
121 >  
122      if (evaluator_.isDynamic()) {
123        seleMan_.setSelectionSet(evaluator_.evaluate());
124      }
# Line 122 | Line 130 | namespace OpenMD {
130        
131        // figure out where that object is:
132        Vector3d pos = sd->getPos();
133 <      currentSnapshot_->wrapVector(pos);
133 >      Vector3d vel = sd->getVel();
134 >      RealType m = sd->getMass();
135  
136        int bin = getBin(pos);
128 –      binCount[bin]++;
137  
130 –      RealType m = sd->getMass();
138        binMass[bin] += m;
139 <      Vector3d vel = sd->getVel();
133 <      binVel[bin] += vel;
139 >      binP[bin] += m * vel;
140        binKE[bin] += 0.5 * (m * vel.lengthSquare());
141        binDof[bin] += 3;
142        
# Line 153 | Line 159 | namespace OpenMD {
159        }
160      }
161      
162 <    for (int i = 0; i < nBins_; i++) {
162 >    for (unsigned int i = 0; i < nBins_; i++) {
163 >
164        if (binDof[i] > 0) {
165          RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
166                                            PhysicalConstants::energyConvert);
167          RealType den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
168            / volume_;
169 <        Vector3d vel = binVel[i] / RealType(binCount[i]);
169 >        Vector3d vel = binP[i] / binMass[i];
170 >
171          dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
172          dynamic_cast<VectorAccumulator *>(velocity->accumulator[i])->add(vel);
173          dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
# Line 209 | Line 217 | namespace OpenMD {
217      data_.push_back(density);
218    }
219  
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);
220  
221 <    RealType KE = 0.5 * (mass * vel.lengthSquare());
219 <    int dof = 3;
221 >  void RNEMDR::processFrame(int istep) {
222  
223 <    if (sd->isDirectional()) {
224 <      Vector3d angMom = sd->getJ();
225 <      Mat3x3d I = sd->getI();
226 <      if (sd->isLinear()) {
227 <        int i = sd->linearAxis();
228 <        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);
223 >    Molecule* mol;
224 >    RigidBody* rb;
225 >    StuntDouble* sd;
226 >    SimInfo::MoleculeIterator mi;
227 >    Molecule::RigidBodyIterator rbIter;
228 >    int i;
229  
230 <    RealType rinner = (RealType)bin * binWidth_;
231 <    RealType router = (RealType)(bin+1) * binWidth_;
232 <    RealType den = mass * 3.0 * PhysicalConstants::densityConvert
233 <      / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));  
230 >    vector<RealType> binMass(nBins_, 0.0);
231 >    vector<Mat3x3d>  binI(nBins_);
232 >    vector<Vector3d> binL(nBins_, V3Zero);
233 >    vector<RealType> binKE(nBins_, 0.0);
234 >    vector<int> binDof(nBins_, 0);
235      
236 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
237 <    dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[bin])->add(aVel);
238 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
236 >    for (mol = info_->beginMolecule(mi); mol != NULL;
237 >         mol = info_->nextMolecule(mi)) {
238 >      
239 >      // change the positions of atoms which belong to the rigidbodies
240 >      
241 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
242 >           rb = mol->nextRigidBody(rbIter)) {
243 >        rb->updateAtomVel();
244 >      }
245 >    }
246 >  
247 >    if (evaluator_.isDynamic()) {
248 >      seleMan_.setSelectionSet(evaluator_.evaluate());
249 >    }
250 >    
251 >    // loop over the selected atoms:
252 >    
253 >    for (sd = seleMan_.beginSelected(i); sd != NULL;
254 >         sd = seleMan_.nextSelected(i)) {
255  
256 +      // figure out where that object is:
257 +      int bin = getBin(sd->getPos() );      
258 +
259 +      if (bin >= 0 && bin < nBins_)  {
260 +
261 +        Vector3d rPos = sd->getPos() - coordinateOrigin_;
262 +        Vector3d vel = sd->getVel();      
263 +        RealType m = sd->getMass();
264 +        Vector3d L = m * cross(rPos, vel);
265 +        Mat3x3d I(0.0);
266 +        I = outProduct(rPos, rPos) * m;
267 +        RealType r2 = rPos.lengthSquare();
268 +        I(0, 0) += m * r2;
269 +        I(1, 1) += m * r2;
270 +        I(2, 2) += m * r2;      
271 +
272 +        binMass[bin] += m;
273 +        binI[bin] += I;
274 +        binL[bin] += L;
275 +        binKE[bin] += 0.5 * (m * vel.lengthSquare());
276 +        binDof[bin] += 3;
277 +        
278 +        if (sd->isDirectional()) {
279 +          Vector3d angMom = sd->getJ();
280 +          Mat3x3d Ia = sd->getI();
281 +          if (sd->isLinear()) {
282 +            int i = sd->linearAxis();
283 +            int j = (i + 1) % 3;
284 +            int k = (i + 2) % 3;
285 +            binKE[bin] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
286 +                                 angMom[k] * angMom[k] / Ia(k, k));
287 +            binDof[bin] += 2;
288 +          } else {
289 +            binKE[bin] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
290 +                                 angMom[1] * angMom[1] / Ia(1, 1) +
291 +                                 angMom[2] * angMom[2] / Ia(2, 2));
292 +            binDof[bin] += 3;
293 +          }
294 +        }
295 +      }
296 +    }
297 +    
298 +    for (unsigned int i = 0; i < nBins_; i++) {
299 +      RealType rinner = (RealType)i * binWidth_;
300 +      RealType router = (RealType)(i+1) * binWidth_;
301 +      if (binDof[i] > 0) {
302 +        RealType temp = 2.0 * binKE[i] / (binDof[i] * PhysicalConstants::kb *
303 +                                          PhysicalConstants::energyConvert);
304 +        RealType den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
305 +          / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
306 +
307 +        Vector3d omega = binI[i].inverse() * binL[i];
308 +
309 +        dynamic_cast<Accumulator *>(temperature->accumulator[i])->add(temp);
310 +        dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[i])->add(omega);
311 +        dynamic_cast<Accumulator *>(density->accumulator[i])->add(den);
312 +        dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
313 +      }
314 +    }
315    }
316 +
317 +
318 +  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
319 +  }
320   }
321  

Diff Legend

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