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 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC

# Line 36 | Line 36
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, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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   /**
# Line 55 | Line 56 | namespace OpenMD {
56  
57   namespace OpenMD {
58  
59 <  void  Snapshot::setHmat(const Mat3x3d& m) {
60 <    hmat_ = m;
61 <    invHmat_ = hmat_.inverse();
59 >  Snapshot::Snapshot(int nAtoms, int nRigidbodies, int nCutoffGroups) :
60 >    atomData(nAtoms), rigidbodyData(nRigidbodies),
61 >    cgData(nCutoffGroups, DataStorage::dslPosition),
62 >    orthoTolerance_(1e-6) {
63      
64 +    frameData.id = -1;                  
65 +    frameData.currentTime = 0;    
66 +    frameData.hmat = Mat3x3d(0.0);            
67 +    frameData.invHmat = Mat3x3d(0.0);          
68 +    frameData.orthoRhombic = false;        
69 +    frameData.bondPotential = 0.0;      
70 +    frameData.bendPotential = 0.0;      
71 +    frameData.torsionPotential = 0.0;  
72 +    frameData.inversionPotential = 0.0;
73 +    frameData.lrPotentials = potVec(0.0);
74 +    frameData.excludedPotentials = potVec(0.0);
75 +    frameData.restraintPotential = 0.0;
76 +    frameData.rawPotential = 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 <    //prepare fortran Hmat
85 <    RealType fortranHmat[9];
86 <    RealType fortranInvHmat[9];
87 <    hmat_.getArray(fortranHmat);
88 <    invHmat_.getArray(fortranInvHmat);
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.volume = 0.0;          
108 >    frameData.thermostat = make_pair(0.0, 0.0);
109 >    frameData.electronicThermostat = make_pair(0.0, 0.0);
110 >    frameData.barostat = Mat3x3d(0.0);              
111 >    frameData.stressTensor = Mat3x3d(0.0);              
112 >    frameData.conductiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
113  
114 +    clearDerivedProperties();
115 +  }
116 +
117 +  void Snapshot::clearDerivedProperties() {
118 +    frameData.totalEnergy = 0.0;    
119 +    frameData.translationalKinetic = 0.0;  
120 +    frameData.rotationalKinetic = 0.0;  
121 +    frameData.kineticEnergy = 0.0;  
122 +    frameData.potentialEnergy = 0.0;
123 +    frameData.shortRangePotential = 0.0;
124 +    frameData.longRangePotential = 0.0;
125 +    frameData.pressure = 0.0;        
126 +    frameData.temperature = 0.0;
127 +    frameData.pressureTensor = Mat3x3d(0.0);  
128 +    frameData.systemDipole = Vector3d(0.0);            
129 +    frameData.convectiveHeatFlux = Vector3d(0.0, 0.0, 0.0);
130 +    frameData.electronicTemperature = 0.0;
131 +    frameData.COM = V3Zero;            
132 +    frameData.COMvel = V3Zero;          
133 +    frameData.COMw = V3Zero;            
134 +
135 +    hasTotalEnergy = false;        
136 +    hasTranslationalKineticEnergy = false;      
137 +    hasRotationalKineticEnergy = false;      
138 +    hasKineticEnergy = false;      
139 +    hasShortRangePotential = false;
140 +    hasLongRangePotential = false;
141 +    hasPotentialEnergy = false;    
142 +    hasVolume = false;        
143 +    hasPressure = false;      
144 +    hasTemperature = false;    
145 +    hasElectronicTemperature = false;
146 +    hasCOM = false;
147 +    hasCOMvel = false;
148 +    hasCOMw = false;
149 +    hasPressureTensor = false;    
150 +    hasSystemDipole = false;      
151 +    hasConvectiveHeatFlux = false;  
152 +    hasInertiaTensor = false;
153 +    hasGyrationalVolume = false;  
154 +    hasHullVolume = false;
155 +    hasConservedQuantity = false;
156 +  }
157 +
158 +  /** Returns the id of this Snapshot */
159 +  int Snapshot::getID() {
160 +    return frameData.id;
161 +  }
162 +  
163 +  /** Sets the id of this Snapshot */
164 +  void Snapshot::setID(int id) {
165 +    frameData.id = id;
166 +  }
167 +  
168 +  int Snapshot::getSize() {
169 +    return atomData.getSize() + rigidbodyData.getSize();
170 +  }
171 +  
172 +  /** Returns the number of atoms */
173 +  int Snapshot::getNumberOfAtoms() {
174 +    return atomData.getSize();
175 +  }
176 +  
177 +  /** Returns the number of rigid bodies */
178 +  int Snapshot::getNumberOfRigidBodies() {
179 +    return rigidbodyData.getSize();
180 +  }
181 +  
182 +  /** Returns the number of rigid bodies */
183 +  int Snapshot::getNumberOfCutoffGroups() {
184 +    return cgData.getSize();
185 +  }
186 +  
187 +  /** Returns the H-Matrix */
188 +  Mat3x3d Snapshot::getHmat() {
189 +    return frameData.hmat;
190 +  }
191 +
192 +  /** Sets the H-Matrix */  
193 +  void Snapshot::setHmat(const Mat3x3d& m) {
194 +    hasVolume = false;
195 +    frameData.hmat = m;
196 +    frameData.invHmat = frameData.hmat.inverse();
197 +    
198      //determine whether the box is orthoTolerance or not
199 <    int oldOrthoRhombic = orthoRhombic_;
199 >    bool oldOrthoRhombic = frameData.orthoRhombic;
200      
201 <    RealType smallDiag = fabs(hmat_(0, 0));
202 <    if(smallDiag > fabs(hmat_(1, 1))) smallDiag = fabs(hmat_(1, 1));
203 <    if(smallDiag > fabs(hmat_(2, 2))) smallDiag = fabs(hmat_(2, 2));    
201 >    RealType smallDiag = fabs(frameData.hmat(0, 0));
202 >    if(smallDiag > fabs(frameData.hmat(1, 1))) smallDiag = fabs(frameData.hmat(1, 1));
203 >    if(smallDiag > fabs(frameData.hmat(2, 2))) smallDiag = fabs(frameData.hmat(2, 2));    
204      RealType tol = smallDiag * orthoTolerance_;
205  
206 <    orthoRhombic_ = 1;
206 >    frameData.orthoRhombic = true;
207  
208      for (int i = 0; i < 3; i++ ) {
209        for (int j = 0 ; j < 3; j++) {
210          if (i != j) {
211 <          if (orthoRhombic_) {
212 <            if ( fabs(hmat_(i, j)) >= tol)
213 <              orthoRhombic_ = 0;
211 >          if (frameData.orthoRhombic) {
212 >            if ( fabs(frameData.hmat(i, j)) >= tol)
213 >              frameData.orthoRhombic = false;
214            }        
215          }
216        }
217      }
218 <
219 <    if( oldOrthoRhombic != orthoRhombic_ ){
220 <
221 <      if( orthoRhombic_ ) {
218 >    
219 >    if( oldOrthoRhombic != frameData.orthoRhombic){
220 >      
221 >      if( frameData.orthoRhombic ) {
222          sprintf( painCave.errMsg,
223                   "OpenMD is switching from the default Non-Orthorhombic\n"
224                   "\tto the faster Orthorhombic periodic boundary computations.\n"
# Line 113 | Line 242 | namespace OpenMD {
242          simError();
243        }
244      }    
245 +  }
246 +  
247 +  /** Returns the inverse H-Matrix */
248 +  Mat3x3d Snapshot::getInvHmat() {
249 +    return frameData.invHmat;
250 +  }
251  
252 <    //notify fortran simulation box has changed
253 <    setFortranBox(fortranHmat, fortranInvHmat, &orthoRhombic_);
252 >  RealType Snapshot::getVolume() {
253 >    if (!hasVolume) {
254 >      frameData.volume = frameData.hmat.determinant();
255 >      hasVolume = true;
256 >    }
257 >    return frameData.volume;
258    }
259  
260 +  void Snapshot::setVolume(RealType vol) {
261 +    hasVolume = true;
262 +    frameData.volume = vol;
263 +  }
264  
265 +  /** Wrap a vector according to periodic boundary conditions */
266    void Snapshot::wrapVector(Vector3d& pos) {
267 <
268 <    int i;
269 <    Vector3d scaled;
270 <
271 <    if( !orthoRhombic_ ){
272 <
273 <      // calc the scaled coordinates.
274 <      scaled = invHmat_* pos;
275 <
276 <      // wrap the scaled coordinates
133 <      for (i = 0; i < 3; ++i) {
134 <        scaled[i] -= roundMe(scaled[i]);
135 <      }
136 <
267 >    
268 >    Vector3d scaled = scaleVector(pos);
269 >    
270 >    for (int i = 0; i < 3; i++)
271 >      scaled[i] -= roundMe(scaled[i]);
272 >    
273 >    if( !frameData.orthoRhombic )
274 >      pos = frameData.hmat * scaled;    
275 >    else {
276 >      
277        // calc the wrapped real coordinates from the wrapped scaled coordinates
278 <      pos = hmat_ * scaled;    
278 >      for (int i=0; i<3; i++) {
279 >        pos[i] = scaled[i] * frameData.hmat(i, i);
280 >      }  
281 >    }
282 >  }
283  
284 <    } else {//if it is orthoRhombic, we could improve efficiency by only caculating the diagonal element
284 >  /** Scaling a vector to multiples of the periodic box */
285 >  inline Vector3d Snapshot::scaleVector(Vector3d& pos) {  
286      
287 <      // calc the scaled coordinates.
143 <      for (i=0; i<3; i++) {
144 <        scaled[i] = pos[i] * invHmat_(i, i);
145 <      }
146 <        
147 <      // wrap the scaled coordinates
148 <      for (i = 0; i < 3; ++i) {
149 <        scaled[i] -= roundMe(scaled[i]);
150 <      }
287 >    Vector3d scaled;
288  
289 <      // calc the wrapped real coordinates from the wrapped scaled coordinates
290 <      for (i=0; i<3; i++) {
291 <        pos[i] = scaled[i] * hmat_(i, i);
292 <      }
293 <        
289 >    if( !frameData.orthoRhombic )
290 >      scaled = frameData.invHmat * pos;
291 >    else {
292 >      // calc the scaled coordinates.
293 >      for (int i=0; i<3; i++)
294 >        scaled[i] = pos[i] * frameData.invHmat(i, i);
295      }
296  
297 +    return scaled;
298    }
299  
300 +  void Snapshot::setCOM(const Vector3d& com) {
301 +    frameData.COM = com;
302 +    hasCOM = true;
303 +  }
304 +  
305 +  void Snapshot::setCOMvel(const Vector3d& comVel) {
306 +    frameData.COMvel = comVel;
307 +    hasCOMvel = true;
308 +  }
309 +  
310 +  void Snapshot::setCOMw(const Vector3d& comw) {
311 +    frameData.COMw = comw;
312 +    hasCOMw = true;
313 +  }
314 +  
315    Vector3d Snapshot::getCOM() {
316 <    if( !hasCOM_ ) {
163 <      sprintf( painCave.errMsg, "COM was requested before COM was computed!\n");
164 <      painCave.severity = OPENMD_ERROR;
165 <      simError();
166 <    }
167 <    return COM_;
316 >    return frameData.COM;
317    }
318    
319    Vector3d Snapshot::getCOMvel() {
320 <    if( !hasCOM_ ) {
172 <      sprintf( painCave.errMsg, "COMvel was requested before COM was computed!\n");
173 <      painCave.severity = OPENMD_ERROR;
174 <      simError();
175 <    }
176 <    return COMvel_;
320 >    return frameData.COMvel;
321    }
322    
323    Vector3d Snapshot::getCOMw() {
324 <    if( !hasCOM_ ) {
181 <      sprintf( painCave.errMsg, "COMw was requested before COM was computed!\n");
182 <      painCave.severity = OPENMD_ERROR;
183 <      simError();
184 <    }
185 <    return COMw_;
324 >    return frameData.COMw;
325    }
326 +  
327 +  RealType Snapshot::getTime() {
328 +    return frameData.currentTime;
329 +  }
330 +  
331 +  void Snapshot::increaseTime(RealType dt) {
332 +    setTime(getTime() + dt);
333 +  }
334 +  
335 +  void Snapshot::setTime(RealType time) {
336 +    frameData.currentTime = time;
337 +  }
338  
339 < }
339 >  void Snapshot::setBondPotential(RealType bp) {
340 >    frameData.bondPotential = bp;
341 >  }
342    
343 +  void Snapshot::setBendPotential(RealType bp) {
344 +    frameData.bendPotential = bp;
345 +  }
346 +  
347 +  void Snapshot::setTorsionPotential(RealType tp) {
348 +    frameData.torsionPotential = tp;
349 +  }
350 +  
351 +  void Snapshot::setInversionPotential(RealType ip) {
352 +    frameData.inversionPotential = ip;
353 +  }
354 +
355 +
356 +  RealType Snapshot::getBondPotential() {
357 +    return frameData.bondPotential;
358 +  }
359 +  RealType Snapshot::getBendPotential() {
360 +    return frameData.bendPotential;
361 +  }
362 +  RealType Snapshot::getTorsionPotential() {
363 +    return frameData.torsionPotential;
364 +  }
365 +  RealType Snapshot::getInversionPotential() {
366 +    return frameData.inversionPotential;
367 +  }
368 +
369 +  RealType Snapshot::getShortRangePotential() {
370 +    if (!hasShortRangePotential) {
371 +      frameData.shortRangePotential = frameData.bondPotential;
372 +      frameData.shortRangePotential += frameData.bendPotential;
373 +      frameData.shortRangePotential += frameData.torsionPotential;
374 +      frameData.shortRangePotential += frameData.inversionPotential;
375 +      hasShortRangePotential = true;
376 +    }
377 +    return frameData.shortRangePotential;
378 +  }
379 +
380 +  void Snapshot::setLongRangePotential(potVec lrPot) {
381 +    frameData.lrPotentials = lrPot;
382 +  }
383 +    
384 +  RealType Snapshot::getLongRangePotential() {
385 +    if (!hasLongRangePotential) {
386 +      for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
387 +        frameData.longRangePotential += frameData.lrPotentials[i];
388 +      }
389 +      hasLongRangePotential = true;
390 +    }  
391 +    return frameData.longRangePotential;
392 +  }
393 +
394 +  potVec Snapshot::getLongRangePotentials() {
395 +    return frameData.lrPotentials;
396 +  }
397 +
398 +  RealType Snapshot::getPotentialEnergy() {
399 +    if (!hasPotentialEnergy) {
400 +      frameData.potentialEnergy = this->getLongRangePotential();
401 +      frameData.potentialEnergy += this->getShortRangePotential();
402 +      hasPotentialEnergy = true;
403 +    }
404 +    return frameData.potentialEnergy;
405 +  }
406 +    
407 +  void Snapshot::setExcludedPotentials(potVec exPot) {
408 +    frameData.excludedPotentials = exPot;
409 +  }
410 +
411 +  potVec Snapshot::getExcludedPotentials() {
412 +    return frameData.excludedPotentials;
413 +  }
414 +      
415 +  void Snapshot::setRestraintPotential(RealType rp) {
416 +    frameData.restraintPotential = rp;
417 +  }
418 +  
419 +  RealType Snapshot::getRestraintPotential() {
420 +    return frameData.restraintPotential;
421 +  }
422 +  
423 +  void Snapshot::setRawPotential(RealType rp) {
424 +    frameData.rawPotential = rp;
425 +  }
426 +  
427 +  RealType Snapshot::getRawPotential() {
428 +    return frameData.rawPotential;
429 +  }
430 +
431 +  RealType Snapshot::getTranslationalKineticEnergy() {
432 +    return frameData.translationalKinetic;
433 +  }
434 +
435 +  RealType Snapshot::getRotationalKineticEnergy() {
436 +    return frameData.rotationalKinetic;
437 +  }
438 +
439 +  RealType Snapshot::getKineticEnergy() {
440 +    return frameData.kineticEnergy;
441 +  }
442 +
443 +  void Snapshot::setTranslationalKineticEnergy(RealType tke) {
444 +    hasTranslationalKineticEnergy = true;
445 +    frameData.translationalKinetic = tke;
446 +  }
447 +
448 +  void Snapshot::setRotationalKineticEnergy(RealType rke) {
449 +    hasRotationalKineticEnergy = true;
450 +    frameData.rotationalKinetic = rke;
451 +  }
452 +
453 +  void Snapshot::setKineticEnergy(RealType ke) {
454 +    hasKineticEnergy = true;
455 +    frameData.kineticEnergy = ke;
456 +  }
457 +
458 +  RealType Snapshot::getTotalEnergy() {
459 +    return frameData.totalEnergy;
460 +  }
461 +
462 +  void Snapshot::setTotalEnergy(RealType te) {
463 +    hasTotalEnergy = true;
464 +    frameData.totalEnergy = te;
465 +  }
466 +
467 +  RealType Snapshot::getConservedQuantity() {
468 +    return frameData.conservedQuantity;
469 +  }
470 +
471 +  void Snapshot::setConservedQuantity(RealType cq) {
472 +    hasConservedQuantity = true;
473 +    frameData.conservedQuantity = cq;
474 +  }
475 +
476 +  RealType Snapshot::getTemperature() {
477 +    return frameData.temperature;
478 +  }
479 +
480 +  void Snapshot::setTemperature(RealType temp) {
481 +    hasTemperature = true;
482 +    frameData.temperature = temp;
483 +  }
484 +
485 +  RealType Snapshot::getElectronicTemperature() {
486 +    return frameData.electronicTemperature;
487 +  }
488 +
489 +  void Snapshot::setElectronicTemperature(RealType eTemp) {
490 +    hasElectronicTemperature = true;
491 +    frameData.electronicTemperature = eTemp;
492 +  }
493 +
494 +  RealType Snapshot::getPressure() {
495 +    return frameData.pressure;
496 +  }
497 +
498 +  void Snapshot::setPressure(RealType pressure) {
499 +    hasPressure = true;
500 +    frameData.pressure = pressure;
501 +  }
502 +
503 +  Mat3x3d Snapshot::getPressureTensor() {
504 +    return frameData.pressureTensor;
505 +  }
506 +
507 +
508 +  void Snapshot::setPressureTensor(const Mat3x3d& pressureTensor) {
509 +    hasPressureTensor = true;
510 +    frameData.pressureTensor = pressureTensor;
511 +  }
512 +
513 +  void Snapshot::setStressTensor(const Mat3x3d& stressTensor) {
514 +    frameData.stressTensor = stressTensor;
515 +  }
516 +
517 +  Mat3x3d  Snapshot::getStressTensor() {
518 +    return frameData.stressTensor;
519 +  }
520 +
521 +  void Snapshot::setConductiveHeatFlux(const Vector3d& chf) {
522 +    frameData.conductiveHeatFlux = chf;
523 +  }
524 +
525 +  Vector3d Snapshot::getConductiveHeatFlux() {
526 +    return frameData.conductiveHeatFlux;
527 +  }
528 +  
529 +  Vector3d Snapshot::getConvectiveHeatFlux() {
530 +    return frameData.convectiveHeatFlux;
531 +  }
532 +
533 +  void Snapshot::setConvectiveHeatFlux(const Vector3d& chf) {    
534 +    hasConvectiveHeatFlux = true;
535 +    frameData.convectiveHeatFlux = chf;
536 +  }
537 +
538 +  Vector3d Snapshot::getHeatFlux() {
539 +    // BE CAREFUL WITH UNITS
540 +    return getConductiveHeatFlux() + getConvectiveHeatFlux();
541 +  }
542 +
543 +  Vector3d Snapshot::getSystemDipole() {
544 +    return frameData.systemDipole;
545 +  }
546 +
547 +  void Snapshot::setSystemDipole(const Vector3d& bd) {    
548 +    hasSystemDipole = true;
549 +    frameData.systemDipole = bd;
550 +  }
551 +
552 +  void Snapshot::setThermostat(const pair<RealType, RealType>& thermostat) {
553 +    frameData.thermostat = thermostat;
554 +  }
555 +
556 +  pair<RealType, RealType> Snapshot::getThermostat() {
557 +    return frameData.thermostat;
558 +  }
559 +
560 +  void Snapshot::setElectronicThermostat(const pair<RealType, RealType>& eTherm) {
561 +    frameData.electronicThermostat = eTherm;
562 +  }
563 +
564 +  pair<RealType, RealType> Snapshot::getElectronicThermostat() {
565 +    return frameData.electronicThermostat;
566 +  }
567 +
568 +  void Snapshot::setBarostat(const Mat3x3d& barostat) {
569 +    frameData.barostat = barostat;
570 +  }
571 +
572 +  Mat3x3d Snapshot::getBarostat() {
573 +    return frameData.barostat;
574 +  }
575 +
576 +  void Snapshot::setInertiaTensor(const Mat3x3d& inertiaTensor) {
577 +    frameData.inertiaTensor = inertiaTensor;
578 +    hasInertiaTensor = true;
579 +  }
580 +
581 +  Mat3x3d Snapshot::getInertiaTensor() {
582 +    return frameData.inertiaTensor;
583 +  }
584 +
585 +  void Snapshot::setGyrationalVolume(const RealType gyrationalVolume) {
586 +    frameData.gyrationalVolume = gyrationalVolume;
587 +    hasGyrationalVolume = true;
588 +  }
589 +
590 +  RealType Snapshot::getGyrationalVolume() {
591 +    return frameData.gyrationalVolume;
592 +  }
593 +
594 +  void Snapshot::setHullVolume(const RealType hullVolume) {
595 +    frameData.hullVolume = hullVolume;
596 +    hasHullVolume = true;
597 +  }
598 +
599 +  RealType Snapshot::getHullVolume() {
600 +    return frameData.hullVolume;
601 +  }
602 +
603 +  void Snapshot::setOrthoTolerance(RealType ot) {
604 +    orthoTolerance_ = ot;
605 +  }
606 + }

Diff Legend

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