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 1892 by gezelter, Wed Jun 19 14:14:56 2013 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

# Line 104 | Line 104 | namespace OpenMD {
104      int i;
105  
106      vector<RealType> binMass(nBins_, 0.0);
107 <    vector<RealType> binPx(nBins_, 0.0);
108 <    vector<RealType> binPy(nBins_, 0.0);
109 <    vector<RealType> binPz(nBins_, 0.0);
107 >    vector<Vector3d> binP(nBins_, V3Zero);
108      vector<RealType> binKE(nBins_, 0.0);
109      vector<unsigned int> binDof(nBins_, 0);
112 –    vector<unsigned int> binCount(nBins_, 0);
110      
114 –    
111      for (mol = info_->beginMolecule(mi); mol != NULL;
112           mol = info_->nextMolecule(mi)) {
113        
# Line 139 | Line 135 | namespace OpenMD {
135  
136        int bin = getBin(pos);
137  
142 –      binCount[bin] += 1;
143 –
138        binMass[bin] += m;
139 <      binPx[bin] += m * vel.x();
146 <      binPy[bin] += m * vel.y();
147 <      binPz[bin] += m * vel.z();
139 >      binP[bin] += m * vel;
140        binKE[bin] += 0.5 * (m * vel.lengthSquare());
141        binDof[bin] += 3;
142        
# Line 167 | Line 159 | namespace OpenMD {
159        }
160      }
161      
162 <    for (unsigned int i = 0; i < nBins_; i++) {
162 >    for (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;
178 <        vel.x() = binPx[i] / binMass[i];
179 <        vel.y() = binPy[i] / binMass[i];
180 <        vel.z() = binPz[i] / binMass[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);
# Line 209 | Line 198 | namespace OpenMD {
198      
199      angularVelocity = new OutputData;
200      angularVelocity->units = "angstroms^2/fs";
201 <    angularVelocity->title =  "Velocity";  
201 >    angularVelocity->title =  "Angular Velocity";  
202      angularVelocity->dataType = odtVector3;
203      angularVelocity->dataHandling = odhAverage;
204      angularVelocity->accumulator.reserve(nBins_);
# Line 239 | Line 228 | namespace OpenMD {
228      int i;
229  
230      vector<RealType> binMass(nBins_, 0.0);
231 <    vector<Vector3d> binaVel(nBins_, V3Zero);
231 >    vector<Mat3x3d>  binI(nBins_);
232 >    vector<Vector3d> binL(nBins_, V3Zero);
233      vector<RealType> binKE(nBins_, 0.0);
234 <    vector<unsigned int> binDof(nBins_, 0);
245 <    vector<unsigned int> binCount(nBins_, 0);
234 >    vector<int> binDof(nBins_, 0);
235      
236      for (mol = info_->beginMolecule(mi); mol != NULL;
237           mol = info_->nextMolecule(mi)) {
# Line 263 | Line 252 | namespace OpenMD {
252      
253      for (sd = seleMan_.beginSelected(i); sd != NULL;
254           sd = seleMan_.nextSelected(i)) {
255 <      
255 >
256        // figure out where that object is:
257 +      int bin = getBin( sd->getPos() );  
258  
259 <      Vector3d rPos = sd->getPos() - coordinateOrigin_;
270 <      Vector3d vel = sd->getVel();      
271 <      Vector3d aVel = cross(rPos, vel);
272 <      RealType m = sd->getMass();
259 >      if (bin >= 0 && bin < nBins_)  {
260  
261 <      int bin = getBin(rPos);
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 <      binCount[bin] += 1;
273 <
274 <      binMass[bin] += m;
275 <      binaVel[bin] += aVel;
276 <      binKE[bin] += 0.5 * (m * vel.lengthSquare());
277 <      binDof[bin] += 3;
278 <      
279 <      if (sd->isDirectional()) {
280 <        Vector3d angMom = sd->getJ();
281 <        Mat3x3d I = sd->getI();
282 <        if (sd->isLinear()) {
283 <          int i = sd->linearAxis();
284 <          int j = (i + 1) % 3;
285 <          int k = (i + 2) % 3;
286 <          binKE[bin] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
287 <                               angMom[k] * angMom[k] / I(k, k));
288 <          binDof[bin] += 2;
289 <        } else {
290 <          binKE[bin] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
291 <                               angMom[1] * angMom[1] / I(1, 1) +
292 <                               angMom[2] * angMom[2] / I(2, 2));
293 <          binDof[bin] += 3;
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++) {
298 >    for (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 <        Vector3d aVel = binaVel[i] / RealType(binCount[i]);
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(aVel);
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        }
# Line 319 | Line 317 | namespace OpenMD {
317  
318    void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
319    }
320 +  
321 +  RNEMDRTheta::RNEMDRTheta(SimInfo* info, const std::string& filename,
322 +                           const std::string& sele, int nrbins, int nangleBins)
323 +    : ShellStatistics(info, filename, sele, nrbins), nAngleBins_(nangleBins) {
324 +    
325 +    Globals* simParams = info->getSimParams();
326 +    RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
327 +    bool hasAngularMomentumFluxVector = rnemdParams->haveAngularMomentumFluxVector();
328 +    
329 +    if (hasAngularMomentumFluxVector) {
330 +      fluxVector_ = rnemdParams->getAngularMomentumFluxVector();
331 +    } else {
332 +      
333 +      std::string fluxStr = rnemdParams->getFluxType();
334 +      if (fluxStr.find("Lx") != std::string::npos) {
335 +        fluxVector_ = V3X;
336 +      } else if (fluxStr.find("Ly") != std::string::npos) {
337 +        fluxVector_ = V3Y;
338 +      } else {
339 +        fluxVector_ = V3Z;
340 +      }
341 +    }
342 +    
343 +    fluxVector_.normalize();
344 +
345 +    setOutputName(getPrefix(filename) + ".rnemdRTheta");
346 +
347 +    angularVelocity = new OutputData;
348 +    angularVelocity->units = "angstroms^2/fs";
349 +    angularVelocity->title =  "Angular Velocity";  
350 +    angularVelocity->dataType = odtArray2d;
351 +    angularVelocity->dataHandling = odhAverage;
352 +    angularVelocity->accumulatorArray2d.reserve(nBins_);
353 +    for (int i = 0; i < nBins_; i++) {
354 +      angularVelocity->accumulatorArray2d[i].reserve(nAngleBins_);
355 +      for (int j = 0 ; j < nAngleBins_; j++) {      
356 +        angularVelocity->accumulatorArray2d[i][j] = new Accumulator();
357 +      }
358 +    }
359 +    data_.push_back(angularVelocity);
360 +
361 +  }
362 +
363 +
364 +  std::pair<int,int> RNEMDRTheta::getBins(Vector3d pos) {
365 +    std::pair<int,int> result;
366 +
367 +    Vector3d rPos = pos - coordinateOrigin_;
368 +    RealType cosAngle= dot(rPos, fluxVector_) / rPos.length();
369 +
370 +    result.first = int(rPos.length() / binWidth_);
371 +    result.second = int( (nAngleBins_ - 1) * 0.5 * (cosAngle + 1.0) );
372 +    return result;
373 +  }
374 +
375 +  void RNEMDRTheta::processStuntDouble(StuntDouble* sd, int bin) {
376 +  }
377 +
378 +  void RNEMDRTheta::processFrame(int istep) {
379 +
380 +    Molecule* mol;
381 +    RigidBody* rb;
382 +    StuntDouble* sd;
383 +    SimInfo::MoleculeIterator mi;
384 +    Molecule::RigidBodyIterator rbIter;
385 +    int i;
386 +
387 +    vector<vector<Mat3x3d> >  binI;
388 +    vector<vector<Vector3d> > binL;
389 +    vector<vector<int> > binCount;
390 +    
391 +    for (mol = info_->beginMolecule(mi); mol != NULL;
392 +         mol = info_->nextMolecule(mi)) {
393 +      
394 +      // change the positions of atoms which belong to the rigidbodies
395 +      
396 +      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
397 +           rb = mol->nextRigidBody(rbIter)) {
398 +        rb->updateAtomVel();
399 +      }
400 +    }
401 +  
402 +    if (evaluator_.isDynamic()) {
403 +      seleMan_.setSelectionSet(evaluator_.evaluate());
404 +    }
405 +    
406 +    // loop over the selected atoms:
407 +    
408 +    for (sd = seleMan_.beginSelected(i); sd != NULL;
409 +         sd = seleMan_.nextSelected(i)) {
410 +
411 +      // figure out where that object is:
412 +      std::pair<int,int> bins = getBins( sd->getPos() );  
413 +
414 +      if (bins.first >= 0 && bins.first < nBins_)  {
415 +        if (bins.second >= 0 && bins.second < nAngleBins_) {
416 +
417 +          Vector3d rPos = sd->getPos() - coordinateOrigin_;
418 +          Vector3d vel = sd->getVel();      
419 +          RealType m = sd->getMass();
420 +          Vector3d L = m * cross(rPos, vel);
421 +          Mat3x3d I(0.0);
422 +          I = outProduct(rPos, rPos) * m;
423 +          RealType r2 = rPos.lengthSquare();
424 +          I(0, 0) += m * r2;
425 +          I(1, 1) += m * r2;
426 +          I(2, 2) += m * r2;      
427 +
428 +          binI[bins.first][bins.second] += I;
429 +          binL[bins.first][bins.second] += L;
430 +          binCount[bins.first][bins.second]++;
431 +        }
432 +      }
433 +    }
434 +  
435 +    
436 +    for (int i = 0; i < nBins_; i++) {
437 +      for (int j = 0; j < nAngleBins_; j++) {
438 +
439 +        if (binCount[i][j] > 0) {
440 +          Vector3d omega = binI[i][j].inverse() * binL[i][j];
441 +          RealType omegaProj = dot(omega, fluxVector_);
442 +
443 +          dynamic_cast<Accumulator *>(angularVelocity->accumulatorArray2d[i][j])->add(omegaProj);
444 +        }
445 +      }
446 +    }
447 +  }
448 +
449 +  void RNEMDRTheta::writeOutput() {
450 +    
451 +    vector<OutputData*>::iterator i;
452 +    OutputData* outputData;
453 +    
454 +    ofstream outStream(outputFilename_.c_str());
455 +    if (outStream.is_open()) {
456 +      
457 +      //write title
458 +      outStream << "# SPATIAL STATISTICS\n";
459 +      outStream << "#";
460 +      
461 +      for(outputData = beginOutputData(i); outputData;
462 +          outputData = nextOutputData(i)) {
463 +        outStream << "\t" << outputData->title <<
464 +          "(" << outputData->units << ")";
465 +        // add some extra tabs for column alignment
466 +        if (outputData->dataType == odtVector3) outStream << "\t\t";
467 +      }
468 +      
469 +      outStream << std::endl;
470 +      
471 +      outStream.precision(8);
472 +      
473 +      for (int j = 0; j < nBins_; j++) {        
474 +        
475 +        int counts = counts_->accumulator[j]->count();
476 +
477 +        if (counts > 0) {
478 +          for(outputData = beginOutputData(i); outputData;
479 +              outputData = nextOutputData(i)) {
480 +            
481 +            int n = outputData->accumulator[j]->count();
482 +            if (n != 0) {
483 +              writeData( outStream, outputData, j );
484 +            }
485 +          }
486 +          outStream << std::endl;
487 +        }
488 +      }
489 +    }
490 +  }
491   }
492  

Diff Legend

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