ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/Snapshot.cpp
Revision: 2022
Committed: Fri Sep 26 22:22:28 2014 UTC (10 years, 7 months ago) by gezelter
File size: 19752 byte(s)
Log Message:
Added support for accumulateBoxQuadrupole flag

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

Properties

Name Value
svn:keywords Author Id Revision Date