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 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2026 by gezelter, Wed Oct 22 12:23:59 2014 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> binP(nBins_, V3Zero);
108 >    vector<RealType> binKE(nBins_, 0.0);
109 >    vector<unsigned int> binDof(nBins_, 0);
110 >    
111 >    for (mol = info_->beginMolecule(mi); mol != NULL;
112 >         mol = info_->nextMolecule(mi)) {
113 >      
114 >      // change the positions of atoms which belong to the rigidbodies
115 >      
116 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
117 >           rb = mol->nextRigidBody(rbIter)) {
118 >        rb->updateAtomVel();
119        }
120      }
121 +  
122 +    if (evaluator_.isDynamic()) {
123 +      seleMan_.setSelectionSet(evaluator_.evaluate());
124 +    }
125      
126 <    RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb *
113 <                                PhysicalConstants::energyConvert);
114 <    RealType den = mass * nBins_ * PhysicalConstants::densityConvert / volume_;
126 >    // loop over the selected atoms:
127      
128 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
129 <    dynamic_cast<VectorAccumulator *>(velocity->accumulator[bin])->add(vel);
130 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
128 >    for (sd = seleMan_.beginSelected(i); sd != NULL;
129 >         sd = seleMan_.nextSelected(i)) {
130 >      
131 >      // figure out where that object is:
132 >      Vector3d pos = sd->getPos();
133 >      Vector3d vel = sd->getVel();
134 >      RealType m = sd->getMass();
135  
136 +      int bin = getBin(pos);
137 +
138 +      binMass[bin] += m;
139 +      binP[bin] += m * vel;
140 +      binKE[bin] += 0.5 * (m * vel.lengthSquare());
141 +      binDof[bin] += 3;
142 +      
143 +      if (sd->isDirectional()) {
144 +        Vector3d angMom = sd->getJ();
145 +        Mat3x3d I = sd->getI();
146 +        if (sd->isLinear()) {
147 +          int i = sd->linearAxis();
148 +          int j = (i + 1) % 3;
149 +          int k = (i + 2) % 3;
150 +          binKE[bin] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
151 +                               angMom[k] * angMom[k] / I(k, k));
152 +          binDof[bin] += 2;
153 +        } else {
154 +          binKE[bin] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
155 +                               angMom[1] * angMom[1] / I(1, 1) +
156 +                               angMom[2] * angMom[2] / I(2, 2));
157 +          binDof[bin] += 3;
158 +        }
159 +      }
160 +    }
161 +    
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 = 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);
174 +        dynamic_cast<Accumulator *>(counts_->accumulator[i])->add(1);
175 +      }
176 +    }
177    }
178 +  
179 +  void RNEMDZ::processStuntDouble(StuntDouble* sd, int bin) {
180 +  }
181  
182    RNEMDR::RNEMDR(SimInfo* info, const std::string& filename,
183                   const std::string& sele, int nrbins)
# Line 138 | 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 157 | Line 217 | namespace OpenMD {
217      data_.push_back(density);
218    }
219  
160 –  void RNEMDR::processStuntDouble(StuntDouble* sd, int bin) {
161 –    RealType mass = sd->getMass();
162 –    Vector3d vel = sd->getVel();
163 –    Vector3d rPos = sd->getPos() - coordinateOrigin_;
164 –    Vector3d aVel = cross(rPos, vel);
220  
221 <    RealType KE = 0.5 * (mass * vel.lengthSquare());
167 <    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;
229 <        int k = (i + 2) % 3;
230 <        KE += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
231 <                     angMom[k] * angMom[k] / I(k, k));
232 <        dof += 2;
233 <      } else {
234 <        KE += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
235 <                     angMom[1] * angMom[1] / I(1, 1) +
236 <                     angMom[2] * angMom[2] / I(2, 2));
237 <        dof += 3;
223 >    Molecule* mol;
224 >    RigidBody* rb;
225 >    StuntDouble* sd;
226 >    SimInfo::MoleculeIterator mi;
227 >    Molecule::RigidBodyIterator rbIter;
228 >    int i;
229 >
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 >    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 <    RealType temp = 2.0 * KE / (dof * PhysicalConstants::kb *
188 <                                PhysicalConstants::energyConvert);
189 <
190 <    RealType rinner = (RealType)bin * binWidth_;
191 <    RealType router = (RealType)(bin+1) * binWidth_;
192 <    RealType den = mass * 3.0 * PhysicalConstants::densityConvert
193 <      / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));  
251 >    // loop over the selected atoms:
252      
253 <    dynamic_cast<Accumulator *>(temperature->accumulator[bin])->add(temp);
254 <    dynamic_cast<VectorAccumulator *>(angularVelocity->accumulator[bin])->add(aVel);
255 <    dynamic_cast<Accumulator *>(density->accumulator[bin])->add(den);
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 (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 +  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 +      std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
331 +      if (amf.size() != 3) {
332 +        sprintf(painCave.errMsg,
333 +                "RNEMDRTheta: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
334 +                "\tthere should be 3 parameters, but %lu were specified.\n",
335 +                amf.size());
336 +        painCave.isFatal = 1;
337 +        simError();      
338 +      }
339 +      fluxVector_.x() = amf[0];
340 +      fluxVector_.y() = amf[1];
341 +      fluxVector_.z() = amf[2];
342 +    } else {
343 +      
344 +      std::string fluxStr = rnemdParams->getFluxType();
345 +      if (fluxStr.find("Lx") != std::string::npos) {
346 +        fluxVector_ = V3X;
347 +      } else if (fluxStr.find("Ly") != std::string::npos) {
348 +        fluxVector_ = V3Y;
349 +      } else {
350 +        fluxVector_ = V3Z;
351 +      }
352 +    }
353 +    
354 +    fluxVector_.normalize();
355 +
356 +    setOutputName(getPrefix(filename) + ".rnemdRTheta");
357 +
358 +    angularVelocity = new OutputData;
359 +    angularVelocity->units = "angstroms^2/fs";
360 +    angularVelocity->title =  "Angular Velocity";  
361 +    angularVelocity->dataType = odtArray2d;
362 +    angularVelocity->dataHandling = odhAverage;
363 +    angularVelocity->accumulatorArray2d.reserve(nBins_);
364 +    for (int i = 0; i < nBins_; i++) {
365 +      angularVelocity->accumulatorArray2d[i].reserve(nAngleBins_);
366 +      for (int j = 0 ; j < nAngleBins_; j++) {      
367 +        angularVelocity->accumulatorArray2d[i][j] = new Accumulator();
368 +      }
369 +    }
370 +    data_.push_back(angularVelocity);
371 +
372 +  }
373 +
374 +
375 +  std::pair<int,int> RNEMDRTheta::getBins(Vector3d pos) {
376 +    std::pair<int,int> result;
377 +
378 +    Vector3d rPos = pos - coordinateOrigin_;
379 +    RealType cosAngle= dot(rPos, fluxVector_) / rPos.length();
380 +
381 +    result.first = int(rPos.length() / binWidth_);
382 +    result.second = int( (nAngleBins_ - 1) * 0.5 * (cosAngle + 1.0) );
383 +    return result;
384 +  }
385 +
386 +  void RNEMDRTheta::processStuntDouble(StuntDouble* sd, int bin) {
387 +  }
388 +
389 +  void RNEMDRTheta::processFrame(int istep) {
390 +
391 +    Molecule* mol;
392 +    RigidBody* rb;
393 +    StuntDouble* sd;
394 +    SimInfo::MoleculeIterator mi;
395 +    Molecule::RigidBodyIterator rbIter;
396 +    int i;
397 +
398 +    vector<vector<Mat3x3d> >  binI;
399 +    vector<vector<Vector3d> > binL;
400 +    vector<vector<int> > binCount;
401 +    
402 +    for (mol = info_->beginMolecule(mi); mol != NULL;
403 +         mol = info_->nextMolecule(mi)) {
404 +      
405 +      // change the positions of atoms which belong to the rigidbodies
406 +      
407 +      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
408 +           rb = mol->nextRigidBody(rbIter)) {
409 +        rb->updateAtomVel();
410 +      }
411 +    }
412 +  
413 +    if (evaluator_.isDynamic()) {
414 +      seleMan_.setSelectionSet(evaluator_.evaluate());
415 +    }
416 +    
417 +    // loop over the selected atoms:
418 +    
419 +    for (sd = seleMan_.beginSelected(i); sd != NULL;
420 +         sd = seleMan_.nextSelected(i)) {
421 +
422 +      // figure out where that object is:
423 +      std::pair<int,int> bins = getBins( sd->getPos() );  
424 +
425 +      if (bins.first >= 0 && bins.first < nBins_)  {
426 +        if (bins.second >= 0 && bins.second < nAngleBins_) {
427 +
428 +          Vector3d rPos = sd->getPos() - coordinateOrigin_;
429 +          Vector3d vel = sd->getVel();      
430 +          RealType m = sd->getMass();
431 +          Vector3d L = m * cross(rPos, vel);
432 +          Mat3x3d I(0.0);
433 +          I = outProduct(rPos, rPos) * m;
434 +          RealType r2 = rPos.lengthSquare();
435 +          I(0, 0) += m * r2;
436 +          I(1, 1) += m * r2;
437 +          I(2, 2) += m * r2;      
438 +
439 +          binI[bins.first][bins.second] += I;
440 +          binL[bins.first][bins.second] += L;
441 +          binCount[bins.first][bins.second]++;
442 +        }
443 +      }
444 +    }
445 +  
446 +    
447 +    for (int i = 0; i < nBins_; i++) {
448 +      for (int j = 0; j < nAngleBins_; j++) {
449 +
450 +        if (binCount[i][j] > 0) {
451 +          Vector3d omega = binI[i][j].inverse() * binL[i][j];
452 +          RealType omegaProj = dot(omega, fluxVector_);
453 +
454 +          dynamic_cast<Accumulator *>(angularVelocity->accumulatorArray2d[i][j])->add(omegaProj);
455 +        }
456 +      }
457 +    }
458 +  }
459 +
460 +  void RNEMDRTheta::writeOutput() {
461 +    
462 +    vector<OutputData*>::iterator i;
463 +    OutputData* outputData;
464 +    
465 +    ofstream outStream(outputFilename_.c_str());
466 +    if (outStream.is_open()) {
467 +      
468 +      //write title
469 +      outStream << "# SPATIAL STATISTICS\n";
470 +      outStream << "#";
471 +      
472 +      for(outputData = beginOutputData(i); outputData;
473 +          outputData = nextOutputData(i)) {
474 +        outStream << "\t" << outputData->title <<
475 +          "(" << outputData->units << ")";
476 +        // add some extra tabs for column alignment
477 +        if (outputData->dataType == odtVector3) outStream << "\t\t";
478 +      }
479 +      
480 +      outStream << std::endl;
481 +      
482 +      outStream.precision(8);
483 +      
484 +      for (int j = 0; j < nBins_; j++) {        
485 +        
486 +        int counts = counts_->accumulator[j]->count();
487 +
488 +        if (counts > 0) {
489 +          for(outputData = beginOutputData(i); outputData;
490 +              outputData = nextOutputData(i)) {
491 +            
492 +            int n = outputData->accumulator[j]->count();
493 +            if (n != 0) {
494 +              writeData( outStream, outputData, j );
495 +            }
496 +          }
497 +          outStream << std::endl;
498 +        }
499 +      }
500 +    }
501 +  }
502   }
503  

Diff Legend

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