ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Thermo.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Thermo.cpp (file contents):
Revision 608 by gezelter, Tue Jul 15 14:45:09 2003 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <math.h>
2   #include <iostream>
3   using namespace std;
4  
# Line 10 | Line 10 | using namespace std;
10   #include "SRI.hpp"
11   #include "Integrator.hpp"
12   #include "simError.h"
13 + #include "MatVec3.h"
14 + #include "ConstraintManager.hpp"
15 + #include "Mat3x3d.hpp"
16  
17   #ifdef IS_MPI
18   #define __C
19   #include "mpiSimulation.hpp"
20   #endif // is_mpi
21  
22 + inline double roundMe( double x ){
23 +          return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
24 + }
25  
26 < #define BASE_SEED 123456789
27 <
28 < Thermo::Thermo( SimInfo* the_entry_plug ) {
23 <  entry_plug = the_entry_plug;
24 <  int baseSeed = BASE_SEED;
26 > Thermo::Thermo( SimInfo* the_info ) {
27 >  info = the_info;
28 >  int baseSeed = the_info->getSeed();
29    
30    gaussStream = new gaussianSPRNG( baseSeed );
31 +
32 +  cpIter = info->consMan->createPairIterator();
33   }
34  
35   Thermo::~Thermo(){
36    delete gaussStream;
37 +  delete cpIter;
38   }
39  
40   double Thermo::getKinetic(){
# Line 36 | Line 43 | double Thermo::getKinetic(){
43    double kinetic;
44    double amass;
45    double aVel[3], aJ[3], I[3][3];
46 <  int j, kl;
46 >  int i, j, k, kl;
47  
41 –  DirectionalAtom *dAtom;
42 –
43 –  int n_atoms;
48    double kinetic_global;
49 <  Atom** atoms;
46 <
49 >  vector<StuntDouble *> integrableObjects = info->integrableObjects;
50    
48 –  n_atoms = entry_plug->n_atoms;
49 –  atoms = entry_plug->atoms;
50 –
51    kinetic = 0.0;
52    kinetic_global = 0.0;
53 –  for( kl=0; kl < n_atoms; kl++ ){
54 –    
55 –    atoms[kl]->getVel(aVel);
56 –    amass = atoms[kl]->getMass();
57 –    
58 –    for (j=0; j < 3; j++)
59 –      kinetic += amass * aVel[j] * aVel[j];
53  
54 <    if( atoms[kl]->isDirectional() ){
55 <            
56 <      dAtom = (DirectionalAtom *)atoms[kl];
54 >  for (kl=0; kl<integrableObjects.size(); kl++) {
55 >    integrableObjects[kl]->getVel(aVel);
56 >    amass = integrableObjects[kl]->getMass();
57  
58 <      dAtom->getJ( aJ );
59 <      dAtom->getI( I );
60 <      
61 <      for (j=0; j<3; j++)
62 <        kinetic += aJ[j]*aJ[j] / I[j][j];
63 <      
64 <    }
58 >   for(j=0; j<3; j++)
59 >      kinetic += amass*aVel[j]*aVel[j];
60 >
61 >   if (integrableObjects[kl]->isDirectional()){
62 >
63 >      integrableObjects[kl]->getJ( aJ );
64 >      integrableObjects[kl]->getI( I );
65 >
66 >      if (integrableObjects[kl]->isLinear()) {
67 >        i = integrableObjects[kl]->linearAxis();
68 >        j = (i+1)%3;
69 >        k = (i+2)%3;
70 >        kinetic += aJ[j]*aJ[j]/I[j][j] + aJ[k]*aJ[k]/I[k][k];
71 >      } else {
72 >        for (j=0; j<3; j++)
73 >          kinetic += aJ[j]*aJ[j] / I[j][j];
74 >      }
75 >   }
76    }
77   #ifdef IS_MPI
78    MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
79                  MPI_SUM, MPI_COMM_WORLD);
80    kinetic = kinetic_global;
81   #endif //is_mpi
82 <
82 >  
83    kinetic = kinetic * 0.5 / e_convert;
84  
85    return kinetic;
# Line 88 | Line 92 | double Thermo::getPotential(){
92    int el, nSRI;
93    Molecule* molecules;
94  
95 <  molecules = entry_plug->molecules;
96 <  nSRI = entry_plug->n_SRI;
95 >  molecules = info->molecules;
96 >  nSRI = info->n_SRI;
97  
98    potential_local = 0.0;
99    potential = 0.0;
100 <  potential_local += entry_plug->lrPot;
100 >  potential_local += info->lrPot;
101  
102 <  for( el=0; el<entry_plug->n_mol; el++ ){    
102 >  for( el=0; el<info->n_mol; el++ ){    
103      potential_local += molecules[el].getPotential();
104    }
105  
# Line 107 | Line 111 | double Thermo::getPotential(){
111    potential = potential_local;
112   #endif // is_mpi
113  
110 – #ifdef IS_MPI
111 –  /*
112 –  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
113 –  */
114 – #endif
115 –
114    return potential;
115   }
116  
# Line 126 | Line 124 | double Thermo::getTemperature(){
124  
125   double Thermo::getTemperature(){
126  
127 <  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
127 >  const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
128    double temperature;
129 <  
130 <  temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb );
129 >
130 >  temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
131    return temperature;
132   }
133  
134 < double Thermo::getEnthalpy() {
134 > double Thermo::getVolume() {
135  
136 <  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
137 <  double u, p, v;
140 <  double press[3][3];
136 >  return info->boxVol;
137 > }
138  
139 <  u = this->getTotalE();
139 > double Thermo::getPressure() {
140  
141 +  // Relies on the calculation of the full molecular pressure tensor
142 +  
143 +  const double p_convert = 1.63882576e8;
144 +  double press[3][3];
145 +  double pressure;
146 +
147    this->getPressureTensor(press);
145 –  p = (press[0][0] + press[1][1] + press[2][2]) / 3.0;
148  
149 <  v = this->getVolume();
149 >  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
150  
151 <  return (u + (p*v)/e_convert);
151 >  return pressure;
152   }
153  
154 < double Thermo::getVolume() {
154 > double Thermo::getPressureX() {
155  
156 <  return entry_plug->boxVol;
156 >  // Relies on the calculation of the full molecular pressure tensor
157 >  
158 >  const double p_convert = 1.63882576e8;
159 >  double press[3][3];
160 >  double pressureX;
161 >
162 >  this->getPressureTensor(press);
163 >
164 >  pressureX = p_convert * press[0][0];
165 >
166 >  return pressureX;
167   }
168  
169 < double Thermo::getPressure() {
169 > double Thermo::getPressureY() {
170  
171    // Relies on the calculation of the full molecular pressure tensor
172    
173    const double p_convert = 1.63882576e8;
174    double press[3][3];
175 <  double pressure;
175 >  double pressureY;
176  
177    this->getPressureTensor(press);
178  
179 <  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
179 >  pressureY = p_convert * press[1][1];
180  
181 <  return pressure;
181 >  return pressureY;
182 > }
183 >
184 > double Thermo::getPressureZ() {
185 >
186 >  // Relies on the calculation of the full molecular pressure tensor
187 >  
188 >  const double p_convert = 1.63882576e8;
189 >  double press[3][3];
190 >  double pressureZ;
191 >
192 >  this->getPressureTensor(press);
193 >
194 >  pressureZ = p_convert * press[2][2];
195 >
196 >  return pressureZ;
197   }
198  
199  
# Line 180 | Line 207 | void Thermo::getPressureTensor(double press[3][3]){
207    double molmass, volume;
208    double vcom[3];
209    double p_local[9], p_global[9];
210 <  int i, j, k, l, nMols;
184 <  Molecule* molecules;
210 >  int i, j, k;
211  
186 –  nMols = entry_plug->n_mol;
187 –  molecules = entry_plug->molecules;
188 –  //tau = entry_plug->tau;
189 –
190 –  // use velocities of molecular centers of mass and molecular masses:
212    for (i=0; i < 9; i++) {    
213      p_local[i] = 0.0;
214      p_global[i] = 0.0;
215    }
216  
217 <  for (i=0; i < nMols; i++) {
197 <    molmass = molecules[i].getCOMvel(vcom);
217 >  // use velocities of integrableObjects and their masses:  
218  
219 +  for (i=0; i < info->integrableObjects.size(); i++) {
220 +
221 +    molmass = info->integrableObjects[i]->getMass();
222 +    
223 +    info->integrableObjects[i]->getVel(vcom);
224 +    
225      p_local[0] += molmass * (vcom[0] * vcom[0]);
226      p_local[1] += molmass * (vcom[0] * vcom[1]);
227      p_local[2] += molmass * (vcom[0] * vcom[2]);
# Line 205 | Line 231 | void Thermo::getPressureTensor(double press[3][3]){
231      p_local[6] += molmass * (vcom[2] * vcom[0]);
232      p_local[7] += molmass * (vcom[2] * vcom[1]);
233      p_local[8] += molmass * (vcom[2] * vcom[2]);
234 +
235    }
236  
237    // Get total for entire system from MPI.
238 <
238 >  
239   #ifdef IS_MPI
240    MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
241   #else
# Line 217 | Line 244 | void Thermo::getPressureTensor(double press[3][3]){
244    }
245   #endif // is_mpi
246  
247 <  volume = entry_plug->boxVol;
247 >  volume = this->getVolume();
248  
249 +
250 +
251    for(i = 0; i < 3; i++) {
252      for (j = 0; j < 3; j++) {
253        k = 3*i + j;
254 <      l = 3*j + i;
226 <      press[i][j] = (p_global[k] - entry_plug->tau[l]*e_convert) / volume;
254 >      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
255      }
256    }
257   }
258  
259   void Thermo::velocitize() {
260    
233 –  double x,y;
261    double aVel[3], aJ[3], I[3][3];
262 <  int i, j, vr, vd; // velocity randomizer loop counters
262 >  int i, j, l, m, n, vr, vd; // velocity randomizer loop counters
263    double vdrift[3];
264    double vbar;
265    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
266    double av2;
267    double kebar;
241 –  int n_atoms;
242 –  Atom** atoms;
243 –  DirectionalAtom* dAtom;
268    double temperature;
269 <  int n_oriented;
246 <  int n_constraints;
269 >  int nobj;
270  
271 <  atoms         = entry_plug->atoms;
249 <  n_atoms       = entry_plug->n_atoms;
250 <  temperature   = entry_plug->target_temp;
251 <  n_oriented    = entry_plug->n_oriented;
252 <  n_constraints = entry_plug->n_constraints;
271 >  nobj = info->integrableObjects.size();
272    
273 <  kebar = kb * temperature * (double)entry_plug->ndf /
255 <    ( 2.0 * (double)entry_plug->ndfRaw );
273 >  temperature   = info->target_temp;
274    
275 <  for(vr = 0; vr < n_atoms; vr++){
275 >  kebar = kb * temperature * (double)info->ndfRaw /
276 >    ( 2.0 * (double)info->ndf );
277 >  
278 >  for(vr = 0; vr < nobj; vr++){
279      
280      // uses equipartition theory to solve for vbar in angstrom/fs
281  
282 <    av2 = 2.0 * kebar / atoms[vr]->getMass();
282 >    av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
283      vbar = sqrt( av2 );
284 <
264 < //     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
265 <    
284 >
285      // picks random velocities from a gaussian distribution
286      // centered on vbar
287  
288      for (j=0; j<3; j++)
289        aVel[j] = vbar * gaussStream->getGaussian();
290      
291 <    atoms[vr]->setVel( aVel );
291 >    info->integrableObjects[vr]->setVel( aVel );
292 >    
293 >    if(info->integrableObjects[vr]->isDirectional()){
294  
295 +      info->integrableObjects[vr]->getI( I );
296 +
297 +      if (info->integrableObjects[vr]->isLinear()) {
298 +
299 +        l= info->integrableObjects[vr]->linearAxis();
300 +        m = (l+1)%3;
301 +        n = (l+2)%3;
302 +
303 +        aJ[l] = 0.0;
304 +        vbar = sqrt( 2.0 * kebar * I[m][m] );
305 +        aJ[m] = vbar * gaussStream->getGaussian();
306 +        vbar = sqrt( 2.0 * kebar * I[n][n] );
307 +        aJ[n] = vbar * gaussStream->getGaussian();
308 +        
309 +      } else {
310 +        for (j = 0 ; j < 3; j++) {
311 +          vbar = sqrt( 2.0 * kebar * I[j][j] );
312 +          aJ[j] = vbar * gaussStream->getGaussian();
313 +        }      
314 +      } // else isLinear
315 +
316 +      info->integrableObjects[vr]->setJ( aJ );
317 +      
318 +    }//isDirectional
319 +
320    }
321  
322    // Get the Center of Mass drift velocity.
# Line 280 | Line 326 | void Thermo::velocitize() {
326    //  Corrects for the center of mass drift.
327    // sums all the momentum and divides by total mass.
328  
329 <  for(vd = 0; vd < n_atoms; vd++){
329 >  for(vd = 0; vd < nobj; vd++){
330      
331 <    atoms[vd]->getVel(aVel);
331 >    info->integrableObjects[vd]->getVel(aVel);
332      
333      for (j=0; j < 3; j++)
334        aVel[j] -= vdrift[j];
335          
336 <    atoms[vd]->setVel( aVel );
336 >    info->integrableObjects[vd]->setVel( aVel );
337    }
292 –  if( n_oriented ){
293 –  
294 –    for( i=0; i<n_atoms; i++ ){
295 –      
296 –      if( atoms[i]->isDirectional() ){
297 –        
298 –        dAtom = (DirectionalAtom *)atoms[i];
299 –        dAtom->getI( I );
300 –        
301 –        for (j = 0 ; j < 3; j++) {
338  
303 –          vbar = sqrt( 2.0 * kebar * I[j][j] );
304 –          aJ[j] = vbar * gaussStream->getGaussian();
305 –
306 –        }      
307 –
308 –        dAtom->setJ( aJ );
309 –
310 –      }
311 –    }  
312 –  }
339   }
340  
341   void Thermo::getCOMVel(double vdrift[3]){
# Line 317 | Line 343 | void Thermo::getCOMVel(double vdrift[3]){
343    double mtot, mtot_local;
344    double aVel[3], amass;
345    double vdrift_local[3];
346 <  int vd, n_atoms, j;
347 <  Atom** atoms;
346 >  int vd, j;
347 >  int nobj;
348  
349 <  // We are very careless here with the distinction between n_atoms and n_local
324 <  // We should really fix this before someone pokes an eye out.
349 >  nobj   = info->integrableObjects.size();
350  
326 –  n_atoms = entry_plug->n_atoms;  
327 –  atoms   = entry_plug->atoms;
328 –
351    mtot_local = 0.0;
352    vdrift_local[0] = 0.0;
353    vdrift_local[1] = 0.0;
354    vdrift_local[2] = 0.0;
355    
356 <  for(vd = 0; vd < n_atoms; vd++){
356 >  for(vd = 0; vd < nobj; vd++){
357      
358 <    amass = atoms[vd]->getMass();
359 <    atoms[vd]->getVel( aVel );
358 >    amass = info->integrableObjects[vd]->getMass();
359 >    info->integrableObjects[vd]->getVel( aVel );
360  
361      for(j = 0; j < 3; j++)
362        vdrift_local[j] += aVel[j] * amass;
# Line 358 | Line 380 | void Thermo::getCOMVel(double vdrift[3]){
380    
381   }
382  
383 + void Thermo::getCOM(double COM[3]){
384 +
385 +  double mtot, mtot_local;
386 +  double aPos[3], amass;
387 +  double COM_local[3];
388 +  int i, j;
389 +  int nobj;
390 +
391 +  mtot_local = 0.0;
392 +  COM_local[0] = 0.0;
393 +  COM_local[1] = 0.0;
394 +  COM_local[2] = 0.0;
395 +
396 +  nobj = info->integrableObjects.size();
397 +  for(i = 0; i < nobj; i++){
398 +    
399 +    amass = info->integrableObjects[i]->getMass();
400 +    info->integrableObjects[i]->getPos( aPos );
401 +
402 +    for(j = 0; j < 3; j++)
403 +      COM_local[j] += aPos[j] * amass;
404 +    
405 +    mtot_local += amass;
406 +  }
407 +
408 + #ifdef IS_MPI
409 +  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
410 +  MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
411 + #else
412 +  mtot = mtot_local;
413 +  for(i = 0; i < 3; i++) {
414 +    COM[i] = COM_local[i];
415 +  }
416 + #endif
417 +    
418 +  for (i = 0; i < 3; i++) {
419 +    COM[i] = COM[i] / mtot;
420 +  }
421 + }
422 +
423 + void Thermo::removeCOMdrift() {
424 +  double vdrift[3], aVel[3];
425 +  int vd, j, nobj;
426 +
427 +  nobj = info->integrableObjects.size();
428 +
429 +  // Get the Center of Mass drift velocity.
430 +
431 +  getCOMVel(vdrift);
432 +  
433 +  //  Corrects for the center of mass drift.
434 +  // sums all the momentum and divides by total mass.
435 +
436 +  for(vd = 0; vd < nobj; vd++){
437 +    
438 +    info->integrableObjects[vd]->getVel(aVel);
439 +    
440 +    for (j=0; j < 3; j++)
441 +      aVel[j] -= vdrift[j];
442 +        
443 +    info->integrableObjects[vd]->setVel( aVel );
444 +  }
445 + }
446 +
447 + void Thermo::removeAngularMomentum(){
448 +  Vector3d vcom;
449 +  Vector3d qcom;
450 +  Vector3d pos;
451 +  Vector3d vel;
452 +  double mass;  
453 +  double xx;
454 +  double yy;
455 +  double zz;
456 +  double xy;
457 +  double xz;
458 +  double yz;
459 +  Vector3d localAngMom;
460 +  Vector3d angMom;
461 +  Vector3d omega;
462 +  vector<StuntDouble *> integrableObjects;
463 +  double localInertiaVec[9];
464 +  double inertiaVec[9];
465 +  vector<Vector3d> qMinusQCom;
466 +  vector<Vector3d> vMinusVCom;
467 +  Mat3x3d inertiaMat;
468 +  Mat3x3d inverseInertiaMat;
469 +  
470 +  integrableObjects = info->integrableObjects;
471 +  qMinusQCom.resize(integrableObjects.size());
472 +  vMinusVCom.resize(integrableObjects.size());
473 +  
474 +  getCOM(qcom.vec);
475 +  getCOMVel(vcom.vec);
476 +        
477 +  //initialize components for inertia tensor
478 +  xx = 0.0;
479 +  yy = 0.0;
480 +  zz = 0.0;
481 +  xy = 0.0;
482 +  xz = 0.0;
483 +  yz = 0.0;
484 +  
485 +   //build components of Inertia tensor
486 +  //
487 +  //       [  Ixx -Ixy  -Ixz ]
488 +  //   J = | -Iyx  Iyy  -Iyz |
489 +  //       [ -Izx -Iyz   Izz ]
490 +  //See Fowles and Cassidy Chapter 9 or Goldstein Chapter 5
491 +  for(size_t i = 0; i < integrableObjects.size(); i++){
492 +    integrableObjects[i]->getPos(pos.vec);
493 +    integrableObjects[i]->getVel(vel.vec);
494 +    mass = integrableObjects[i]->getMass();
495 +    
496 +    qMinusQCom[i] = pos - qcom;
497 +    info->wrapVector(qMinusQCom[i].vec);
498 +    
499 +    vMinusVCom[i] = vel - vcom;
500 +
501 +    //compute moment of inertia coefficents
502 +    xx += qMinusQCom[i].x * qMinusQCom[i].x * mass;
503 +    yy += qMinusQCom[i].y * qMinusQCom[i].y * mass;
504 +    zz += qMinusQCom[i].z * qMinusQCom[i].z * mass;
505 +
506 +    // compute products of inertia
507 +    xy += qMinusQCom[i].x * qMinusQCom[i].y * mass;
508 +    xz += qMinusQCom[i].x * qMinusQCom[i].z * mass;
509 +    yz += qMinusQCom[i].y * qMinusQCom[i].z * mass;
510 +
511 +    localAngMom += crossProduct(qMinusQCom[i] , vMinusVCom[i] ) * mass;
512 +    
513 +  }
514 +
515 +  localInertiaVec[0] =yy+zz;
516 +  localInertiaVec[1] = -xy;
517 +  localInertiaVec[2] = -xz;
518 +  localInertiaVec[3] = -xy;
519 +  localInertiaVec[4] = xx+zz;
520 +  localInertiaVec[5] = -yz;
521 +  localInertiaVec[6] = -xz;
522 +  localInertiaVec[7] = -yz;
523 +  localInertiaVec[8] = xx+yy;
524 +
525 +  //Sum and distribute inertia and angmom arrays
526 + #ifdef MPI
527 +
528 +  MPI_Allreduce(localInertiaVec, inertiaVec, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
529 +
530 +  MPI_Allreduce(localAngMom.vec, angMom.vec, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
531 +
532 +  inertiaMat.element[0][0] = inertiaVec[0];
533 +  inertiaMat.element[0][1] = inertiaVec[1];
534 +  inertiaMat.element[0][2] = inertiaVec[2];
535 +
536 +  inertiaMat.element[1][0] = inertiaVec[3];
537 +  inertiaMat.element[1][1] = inertiaVec[4];
538 +  inertiaMat.element[1][2] = inertiaVec[5];
539 +
540 +  inertiaMat.element[2][0] = inertiaVec[6];
541 +  inertiaMat.element[2][1] = inertiaVec[7];
542 +  inertiaMat.element[2][2] = inertiaVec[8];
543 +
544 + #else
545 +
546 +    inertiaMat.element[0][0] = localInertiaVec[0];
547 +    inertiaMat.element[0][1] = localInertiaVec[1];
548 +    inertiaMat.element[0][2] = localInertiaVec[2];
549 +
550 +    inertiaMat.element[1][0] = localInertiaVec[3];
551 +    inertiaMat.element[1][1] = localInertiaVec[4];
552 +    inertiaMat.element[1][2] = localInertiaVec[5];
553 +
554 +    inertiaMat.element[2][0] = localInertiaVec[6];
555 +    inertiaMat.element[2][1] = localInertiaVec[7];
556 +    inertiaMat.element[2][2] = localInertiaVec[8];
557 +  
558 +    angMom     = localAngMom;
559 + #endif
560 +
561 +    //invert the moment of inertia tensor by LU-decomposition / backsolving:
562 +
563 +    inverseInertiaMat = inertiaMat.inverse();
564 +
565 +    //calculate the angular velocities: omega = I^-1 . L
566 +
567 +    omega = inverseInertiaMat * angMom;
568 +
569 +    //subtract out center of mass velocity and angular momentum from
570 +    //particle velocities
571 +
572 +    for(size_t i = 0; i < integrableObjects.size(); i++){
573 +      vel = vMinusVCom[i] - crossProduct(omega, qMinusQCom[i]);
574 +      integrableObjects[i]->setVel(vel.vec);      
575 +    }
576 + }
577 +
578 + double Thermo::getConsEnergy(){
579 +  ConstraintPair* consPair;
580 +  double totConsEnergy;
581 +  double bondLen2;
582 +  double dist;
583 +  double lamda;
584 +  
585 +  totConsEnergy = 0;
586 +  
587 +  for(cpIter->first(); !cpIter->isEnd(); cpIter->next()){
588 +    consPair =  cpIter->currentItem();
589 +    bondLen2 = consPair->getBondLength2();
590 +    lamda = consPair->getLamda();
591 +    //dist = consPair->getDistance();
592 +
593 +    //totConsEnergy += lamda * (dist*dist - bondLen2);
594 +  }
595 +
596 +  return totConsEnergy;
597 + }
598 +
599 +

Diff Legend

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