ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/RNEMDStats.cpp
Revision: 2026
Committed: Wed Oct 22 12:23:59 2014 UTC (10 years, 6 months ago) by gezelter
File size: 16786 byte(s)
Log Message:
Starting to add support for UniformGradient. 
Changed Vector3d input type to a more general std::vector<RealType> input.  This change alters RNEMD and UniformField inputs.

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 */
41
42
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 {
50
51 RNEMDZ::RNEMDZ(SimInfo* info, const std::string& filename,
52 const std::string& sele, int nzbins)
53 : SlabStatistics(info, filename, sele, nzbins) {
54
55 setOutputName(getPrefix(filename) + ".rnemdZ");
56
57 temperature = new OutputData;
58 temperature->units = "K";
59 temperature->title = "Temperature";
60 temperature->dataType = odtReal;
61 temperature->dataHandling = odhAverage;
62 temperature->accumulator.reserve(nBins_);
63 for (int i = 0; i < nBins_; i++)
64 temperature->accumulator.push_back( new Accumulator() );
65 data_.push_back(temperature);
66
67 velocity = new OutputData;
68 velocity->units = "angstroms/fs";
69 velocity->title = "Velocity";
70 velocity->dataType = odtVector3;
71 velocity->dataHandling = odhAverage;
72 velocity->accumulator.reserve(nBins_);
73 for (int i = 0; i < nBins_; i++)
74 velocity->accumulator.push_back( new VectorAccumulator() );
75 data_.push_back(velocity);
76
77 density = new OutputData;
78 density->units = "g cm^-3";
79 density->title = "Density";
80 density->dataType = odtReal;
81 density->dataHandling = odhAverage;
82 density->accumulator.reserve(nBins_);
83 for (int i = 0; i < nBins_; i++)
84 density->accumulator.push_back( new Accumulator() );
85 data_.push_back(density);
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;
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 // loop over the selected atoms:
127
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)
184 : ShellStatistics(info, filename, sele, nrbins) {
185
186
187 setOutputName(getPrefix(filename) + ".rnemdR");
188
189 temperature = new OutputData;
190 temperature->units = "K";
191 temperature->title = "Temperature";
192 temperature->dataType = odtReal;
193 temperature->dataHandling = odhAverage;
194 temperature->accumulator.reserve(nBins_);
195 for (int i = 0; i < nBins_; i++)
196 temperature->accumulator.push_back( new Accumulator() );
197 data_.push_back(temperature);
198
199 angularVelocity = new OutputData;
200 angularVelocity->units = "angstroms^2/fs";
201 angularVelocity->title = "Angular Velocity";
202 angularVelocity->dataType = odtVector3;
203 angularVelocity->dataHandling = odhAverage;
204 angularVelocity->accumulator.reserve(nBins_);
205 for (int i = 0; i < nBins_; i++)
206 angularVelocity->accumulator.push_back( new VectorAccumulator() );
207 data_.push_back(angularVelocity);
208
209 density = new OutputData;
210 density->units = "g cm^-3";
211 density->title = "Density";
212 density->dataType = odtReal;
213 density->dataHandling = odhAverage;
214 density->accumulator.reserve(nBins_);
215 for (int i = 0; i < nBins_; i++)
216 density->accumulator.push_back( new Accumulator() );
217 data_.push_back(density);
218 }
219
220
221 void RNEMDR::processFrame(int istep) {
222
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 // 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 (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

Properties

Name Value
svn:eol-style native