ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Snapshot.cpp
(Generate patch)

Comparing branches/development/src/brains/Snapshot.cpp (file contents):
Revision 1715 by gezelter, Tue May 22 21:55:31 2012 UTC vs.
Revision 1838 by gezelter, Tue Jan 22 16:20:11 2013 UTC

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

Diff Legend

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