ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/Snapshot.cpp
Revision: 1896
Committed: Tue Jul 2 20:02:31 2013 UTC (11 years, 10 months ago) by gezelter
File size: 18979 byte(s)
Log Message:
Speedup tests

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 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 */
42
43 /**
44 * @file Snapshot.cpp
45 * @author tlin
46 * @date 11/11/2004
47 * @version 1.0
48 */
49
50 #include "brains/Snapshot.hpp"
51 #include "utils/NumericConstant.hpp"
52 #include "utils/simError.h"
53 #include "utils/Utility.hpp"
54 #include <cstdio>
55
56 namespace OpenMD {
57
58 Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups) :
59 atomData(nAtoms), rigidbodyData(nRigidbodies),
60 cgData(nCutoffGroups, DataStorage::dslPosition),
61 orthoTolerance_(1e-6) {
62
63 frameData.id = -1;
64 frameData.currentTime = 0;
65 frameData.hmat = Mat3x3d(0.0);
66 frameData.invHmat = Mat3x3d(0.0);
67 frameData.orthoRhombic = false;
68 frameData.bondPotential = 0.0;
69 frameData.bendPotential = 0.0;
70 frameData.torsionPotential = 0.0;
71 frameData.inversionPotential = 0.0;
72 frameData.lrPotentials = potVec(0.0);
73 frameData.excludedPotentials = potVec(0.0);
74 frameData.restraintPotential = 0.0;
75 frameData.rawPotential = 0.0;
76 frameData.xyArea = 0.0;
77 frameData.volume = 0.0;
78 frameData.thermostat = make_pair(0.0, 0.0);
79 frameData.electronicThermostat = make_pair(0.0, 0.0);
80 frameData.barostat = Mat3x3d(0.0);
81 frameData.stressTensor = Mat3x3d(0.0);
82 frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
83
84 clearDerivedProperties();
85 }
86
87 Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups,
88 int storageLayout) :
89 atomData(nAtoms, storageLayout),
90 rigidbodyData(nRigidbodies, storageLayout),
91 cgData(nCutoffGroups, DataStorage::dslPosition),
92 orthoTolerance_(1e-6) {
93
94 frameData.id = -1;
95 frameData.currentTime = 0;
96 frameData.hmat = Mat3x3d(0.0);
97 frameData.invHmat = Mat3x3d(0.0);
98 frameData.bBox = Mat3x3d(0.0);
99 frameData.invBbox = Mat3x3d(0.0);
100 frameData.orthoRhombic = false;
101 frameData.bondPotential = 0.0;
102 frameData.bendPotential = 0.0;
103 frameData.torsionPotential = 0.0;
104 frameData.inversionPotential = 0.0;
105 frameData.lrPotentials = potVec(0.0);
106 frameData.excludedPotentials = potVec(0.0);
107 frameData.restraintPotential = 0.0;
108 frameData.rawPotential = 0.0;
109 frameData.xyArea = 0.0;
110 frameData.volume = 0.0;
111 frameData.thermostat = make_pair(0.0, 0.0);
112 frameData.electronicThermostat = make_pair(0.0, 0.0);
113 frameData.barostat = Mat3x3d(0.0);
114 frameData.stressTensor = Mat3x3d(0.0);
115 frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
116
117 clearDerivedProperties();
118 }
119
120 void Snapshot::clearDerivedProperties() {
121 frameData.totalEnergy = 0.0;
122 frameData.translationalKinetic = 0.0;
123 frameData.rotationalKinetic = 0.0;
124 frameData.kineticEnergy = 0.0;
125 frameData.potentialEnergy = 0.0;
126 frameData.shortRangePotential = 0.0;
127 frameData.longRangePotential = 0.0;
128 frameData.pressure = 0.0;
129 frameData.temperature = 0.0;
130 frameData.pressureTensor = Mat3x3d(0.0);
131 frameData.systemDipole = Vector3d(0.0);
132 frameData.convectiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
133 frameData.electronicTemperature = 0.0;
134 frameData.COM = V3Zero;
135 frameData.COMvel = V3Zero;
136 frameData.COMw = V3Zero;
137
138 hasTotalEnergy = false;
139 hasTranslationalKineticEnergy = false;
140 hasRotationalKineticEnergy = false;
141 hasKineticEnergy = false;
142 hasShortRangePotential = false;
143 hasLongRangePotential = false;
144 hasPotentialEnergy = false;
145 hasXYarea = false;
146 hasVolume = false;
147 hasPressure = false;
148 hasTemperature = false;
149 hasElectronicTemperature = false;
150 hasCOM = false;
151 hasCOMvel = false;
152 hasCOMw = false;
153 hasPressureTensor = false;
154 hasSystemDipole = false;
155 hasConvectiveHeatFlux = false;
156 hasInertiaTensor = false;
157 hasGyrationalVolume = false;
158 hasHullVolume = false;
159 hasConservedQuantity = false;
160 hasBoundingBox = false;
161 }
162
163 /** Returns the id of this Snapshot */
164 int Snapshot::getID() {
165 return frameData.id;
166 }
167
168 /** Sets the id of this Snapshot */
169 void Snapshot::setID(int id) {
170 frameData.id = id;
171 }
172
173 int Snapshot::getSize() {
174 return atomData.getSize() + rigidbodyData.getSize();
175 }
176
177 /** Returns the number of atoms */
178 int Snapshot::getNumberOfAtoms() {
179 return atomData.getSize();
180 }
181
182 /** Returns the number of rigid bodies */
183 int Snapshot::getNumberOfRigidBodies() {
184 return rigidbodyData.getSize();
185 }
186
187 /** Returns the number of rigid bodies */
188 int Snapshot::getNumberOfCutoffGroups() {
189 return cgData.getSize();
190 }
191
192 /** Returns the H-Matrix */
193 Mat3x3d Snapshot::getHmat() {
194 return frameData.hmat;
195 }
196
197 /** Sets the H-Matrix */
198 void Snapshot::setHmat(const Mat3x3d& m) {
199 hasVolume = false;
200 frameData.hmat = m;
201 frameData.invHmat = frameData.hmat.inverse();
202
203 //determine whether the box is orthoTolerance or not
204 bool oldOrthoRhombic = frameData.orthoRhombic;
205
206 RealType smallDiag = fabs(frameData.hmat(0, 0));
207 if(smallDiag > fabs(frameData.hmat(1, 1))) smallDiag = fabs(frameData.hmat(1, 1));
208 if(smallDiag > fabs(frameData.hmat(2, 2))) smallDiag = fabs(frameData.hmat(2, 2));
209 RealType tol = smallDiag * orthoTolerance_;
210
211 frameData.orthoRhombic = true;
212
213 for (int i = 0; i < 3; i++ ) {
214 for (int j = 0 ; j < 3; j++) {
215 if (i != j) {
216 if (frameData.orthoRhombic) {
217 if ( fabs(frameData.hmat(i, j)) >= tol)
218 frameData.orthoRhombic = false;
219 }
220 }
221 }
222 }
223
224 if( oldOrthoRhombic != frameData.orthoRhombic){
225
226 // It is finally time to suppress these warnings once and for
227 // all. They were annoying and not very informative.
228
229 // if( frameData.orthoRhombic ) {
230 // sprintf( painCave.errMsg,
231 // "OpenMD is switching from the default Non-Orthorhombic\n"
232 // "\tto the faster Orthorhombic periodic boundary computations.\n"
233 // "\tThis is usually a good thing, but if you want the\n"
234 // "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
235 // "\tvariable ( currently set to %G ) smaller.\n",
236 // orthoTolerance_);
237 // painCave.severity = OPENMD_INFO;
238 // simError();
239 // }
240 // else {
241 // sprintf( painCave.errMsg,
242 // "OpenMD is switching from the faster Orthorhombic to the more\n"
243 // "\tflexible Non-Orthorhombic periodic boundary computations.\n"
244 // "\tThis is usually because the box has deformed under\n"
245 // "\tNPTf integration. If you want to live on the edge with\n"
246 // "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
247 // "\tvariable ( currently set to %G ) larger.\n",
248 // orthoTolerance_);
249 // painCave.severity = OPENMD_WARNING;
250 // simError();
251 // }
252 }
253 }
254
255 /** Returns the inverse H-Matrix */
256 Mat3x3d Snapshot::getInvHmat() {
257 return frameData.invHmat;
258 }
259
260 /** Returns the Bounding Box */
261 Mat3x3d Snapshot::getBoundingBox() {
262 return frameData.bBox;
263 }
264
265 /** Sets the Bounding Box */
266 void Snapshot::setBoundingBox(const Mat3x3d& m) {
267 frameData.bBox = m;
268 frameData.invBbox = frameData.bBox.inverse();
269 hasBoundingBox = true;
270 }
271
272 /** Returns the inverse Bounding Box */
273 Mat3x3d Snapshot::getInvBoundingBox() {
274 return frameData.invBbox;
275 }
276
277 RealType Snapshot::getXYarea() {
278 if (!hasXYarea) {
279 Vector3d x = frameData.hmat.getColumn(0);
280 Vector3d y = frameData.hmat.getColumn(1);
281 frameData.xyArea = cross(x,y).length();
282 hasXYarea = true;
283 }
284 return frameData.xyArea;
285 }
286
287 RealType Snapshot::getVolume() {
288 if (!hasVolume) {
289 frameData.volume = frameData.hmat.determinant();
290 hasVolume = true;
291 }
292 return frameData.volume;
293 }
294
295 void Snapshot::setVolume(RealType vol) {
296 hasVolume = true;
297 frameData.volume = vol;
298 }
299
300
301 /** Wrap a vector according to periodic boundary conditions */
302 void Snapshot::wrapVector(Vector3d& pos) {
303
304 if( !frameData.orthoRhombic ) {
305 Vector3d scaled = frameData.invHmat * pos;
306 for (int i = 0; i < 3; i++) {
307 scaled[i] -= roundMe( scaled[i] );
308 }
309 // calc the wrapped real coordinates from the wrapped scaled coordinates
310 pos = frameData.hmat * scaled;
311 } else {
312 RealType scaled;
313 for (int i=0; i<3; i++) {
314 scaled = pos[i] * frameData.invHmat(i,i);
315 scaled -= roundMe( scaled );
316 pos[i] = scaled * frameData.hmat(i,i);
317 }
318 }
319 }
320
321 /** Scaling a vector to multiples of the periodic box */
322 inline Vector3d Snapshot::scaleVector(Vector3d& pos) {
323
324 Vector3d scaled;
325
326 if( !frameData.orthoRhombic )
327 scaled = frameData.invHmat * pos;
328 else {
329 // calc the scaled coordinates.
330 for (int i=0; i<3; i++)
331 scaled[i] = pos[i] * frameData.invHmat(i, i);
332 }
333
334 return scaled;
335 }
336
337 void Snapshot::setCOM(const Vector3d& com) {
338 frameData.COM = com;
339 hasCOM = true;
340 }
341
342 void Snapshot::setCOMvel(const Vector3d& comVel) {
343 frameData.COMvel = comVel;
344 hasCOMvel = true;
345 }
346
347 void Snapshot::setCOMw(const Vector3d& comw) {
348 frameData.COMw = comw;
349 hasCOMw = true;
350 }
351
352 Vector3d Snapshot::getCOM() {
353 return frameData.COM;
354 }
355
356 Vector3d Snapshot::getCOMvel() {
357 return frameData.COMvel;
358 }
359
360 Vector3d Snapshot::getCOMw() {
361 return frameData.COMw;
362 }
363
364 RealType Snapshot::getTime() {
365 return frameData.currentTime;
366 }
367
368 void Snapshot::increaseTime(RealType dt) {
369 setTime(getTime() + dt);
370 }
371
372 void Snapshot::setTime(RealType time) {
373 frameData.currentTime = time;
374 }
375
376 void Snapshot::setBondPotential(RealType bp) {
377 frameData.bondPotential = bp;
378 }
379
380 void Snapshot::setBendPotential(RealType bp) {
381 frameData.bendPotential = bp;
382 }
383
384 void Snapshot::setTorsionPotential(RealType tp) {
385 frameData.torsionPotential = tp;
386 }
387
388 void Snapshot::setInversionPotential(RealType ip) {
389 frameData.inversionPotential = ip;
390 }
391
392
393 RealType Snapshot::getBondPotential() {
394 return frameData.bondPotential;
395 }
396 RealType Snapshot::getBendPotential() {
397 return frameData.bendPotential;
398 }
399 RealType Snapshot::getTorsionPotential() {
400 return frameData.torsionPotential;
401 }
402 RealType Snapshot::getInversionPotential() {
403 return frameData.inversionPotential;
404 }
405
406 RealType Snapshot::getShortRangePotential() {
407 if (!hasShortRangePotential) {
408 frameData.shortRangePotential = frameData.bondPotential;
409 frameData.shortRangePotential += frameData.bendPotential;
410 frameData.shortRangePotential += frameData.torsionPotential;
411 frameData.shortRangePotential += frameData.inversionPotential;
412 hasShortRangePotential = true;
413 }
414 return frameData.shortRangePotential;
415 }
416
417 void Snapshot::setLongRangePotential(potVec lrPot) {
418 frameData.lrPotentials = lrPot;
419 }
420
421 RealType Snapshot::getLongRangePotential() {
422 if (!hasLongRangePotential) {
423 for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
424 frameData.longRangePotential += frameData.lrPotentials[i];
425 }
426 hasLongRangePotential = true;
427 }
428 return frameData.longRangePotential;
429 }
430
431 potVec Snapshot::getLongRangePotentials() {
432 return frameData.lrPotentials;
433 }
434
435 RealType Snapshot::getPotentialEnergy() {
436 if (!hasPotentialEnergy) {
437 frameData.potentialEnergy = this->getLongRangePotential();
438 frameData.potentialEnergy += this->getShortRangePotential();
439 hasPotentialEnergy = true;
440 }
441 return frameData.potentialEnergy;
442 }
443
444 void Snapshot::setExcludedPotentials(potVec exPot) {
445 frameData.excludedPotentials = exPot;
446 }
447
448 potVec Snapshot::getExcludedPotentials() {
449 return frameData.excludedPotentials;
450 }
451
452 void Snapshot::setRestraintPotential(RealType rp) {
453 frameData.restraintPotential = rp;
454 }
455
456 RealType Snapshot::getRestraintPotential() {
457 return frameData.restraintPotential;
458 }
459
460 void Snapshot::setRawPotential(RealType rp) {
461 frameData.rawPotential = rp;
462 }
463
464 RealType Snapshot::getRawPotential() {
465 return frameData.rawPotential;
466 }
467
468 RealType Snapshot::getTranslationalKineticEnergy() {
469 return frameData.translationalKinetic;
470 }
471
472 RealType Snapshot::getRotationalKineticEnergy() {
473 return frameData.rotationalKinetic;
474 }
475
476 RealType Snapshot::getKineticEnergy() {
477 return frameData.kineticEnergy;
478 }
479
480 void Snapshot::setTranslationalKineticEnergy(RealType tke) {
481 hasTranslationalKineticEnergy = true;
482 frameData.translationalKinetic = tke;
483 }
484
485 void Snapshot::setRotationalKineticEnergy(RealType rke) {
486 hasRotationalKineticEnergy = true;
487 frameData.rotationalKinetic = rke;
488 }
489
490 void Snapshot::setKineticEnergy(RealType ke) {
491 hasKineticEnergy = true;
492 frameData.kineticEnergy = ke;
493 }
494
495 RealType Snapshot::getTotalEnergy() {
496 return frameData.totalEnergy;
497 }
498
499 void Snapshot::setTotalEnergy(RealType te) {
500 hasTotalEnergy = true;
501 frameData.totalEnergy = te;
502 }
503
504 RealType Snapshot::getConservedQuantity() {
505 return frameData.conservedQuantity;
506 }
507
508 void Snapshot::setConservedQuantity(RealType cq) {
509 hasConservedQuantity = true;
510 frameData.conservedQuantity = cq;
511 }
512
513 RealType Snapshot::getTemperature() {
514 return frameData.temperature;
515 }
516
517 void Snapshot::setTemperature(RealType temp) {
518 hasTemperature = true;
519 frameData.temperature = temp;
520 }
521
522 RealType Snapshot::getElectronicTemperature() {
523 return frameData.electronicTemperature;
524 }
525
526 void Snapshot::setElectronicTemperature(RealType eTemp) {
527 hasElectronicTemperature = true;
528 frameData.electronicTemperature = eTemp;
529 }
530
531 RealType Snapshot::getPressure() {
532 return frameData.pressure;
533 }
534
535 void Snapshot::setPressure(RealType pressure) {
536 hasPressure = true;
537 frameData.pressure = pressure;
538 }
539
540 Mat3x3d Snapshot::getPressureTensor() {
541 return frameData.pressureTensor;
542 }
543
544
545 void Snapshot::setPressureTensor(const Mat3x3d& pressureTensor) {
546 hasPressureTensor = true;
547 frameData.pressureTensor = pressureTensor;
548 }
549
550 void Snapshot::setStressTensor(const Mat3x3d& stressTensor) {
551 frameData.stressTensor = stressTensor;
552 }
553
554 Mat3x3d Snapshot::getStressTensor() {
555 return frameData.stressTensor;
556 }
557
558 void Snapshot::setConductiveHeatFlux(const Vector3d& chf) {
559 frameData.conductiveHeatFlux = chf;
560 }
561
562 Vector3d Snapshot::getConductiveHeatFlux() {
563 return frameData.conductiveHeatFlux;
564 }
565
566 Vector3d Snapshot::getConvectiveHeatFlux() {
567 return frameData.convectiveHeatFlux;
568 }
569
570 void Snapshot::setConvectiveHeatFlux(const Vector3d& chf) {
571 hasConvectiveHeatFlux = true;
572 frameData.convectiveHeatFlux = chf;
573 }
574
575 Vector3d Snapshot::getHeatFlux() {
576 // BE CAREFUL WITH UNITS
577 return getConductiveHeatFlux() + getConvectiveHeatFlux();
578 }
579
580 Vector3d Snapshot::getSystemDipole() {
581 return frameData.systemDipole;
582 }
583
584 void Snapshot::setSystemDipole(const Vector3d& bd) {
585 hasSystemDipole = true;
586 frameData.systemDipole = bd;
587 }
588
589 void Snapshot::setThermostat(const pair<RealType, RealType>& thermostat) {
590 frameData.thermostat = thermostat;
591 }
592
593 pair<RealType, RealType> Snapshot::getThermostat() {
594 return frameData.thermostat;
595 }
596
597 void Snapshot::setElectronicThermostat(const pair<RealType, RealType>& eTherm) {
598 frameData.electronicThermostat = eTherm;
599 }
600
601 pair<RealType, RealType> Snapshot::getElectronicThermostat() {
602 return frameData.electronicThermostat;
603 }
604
605 void Snapshot::setBarostat(const Mat3x3d& barostat) {
606 frameData.barostat = barostat;
607 }
608
609 Mat3x3d Snapshot::getBarostat() {
610 return frameData.barostat;
611 }
612
613 void Snapshot::setInertiaTensor(const Mat3x3d& inertiaTensor) {
614 frameData.inertiaTensor = inertiaTensor;
615 hasInertiaTensor = true;
616 }
617
618 Mat3x3d Snapshot::getInertiaTensor() {
619 return frameData.inertiaTensor;
620 }
621
622 void Snapshot::setGyrationalVolume(const RealType gyrationalVolume) {
623 frameData.gyrationalVolume = gyrationalVolume;
624 hasGyrationalVolume = true;
625 }
626
627 RealType Snapshot::getGyrationalVolume() {
628 return frameData.gyrationalVolume;
629 }
630
631 void Snapshot::setHullVolume(const RealType hullVolume) {
632 frameData.hullVolume = hullVolume;
633 hasHullVolume = true;
634 }
635
636 RealType Snapshot::getHullVolume() {
637 return frameData.hullVolume;
638 }
639
640 void Snapshot::setOrthoTolerance(RealType ot) {
641 orthoTolerance_ = ot;
642 }
643 }

Properties

Name Value
svn:keywords Author Id Revision Date