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 574 by gezelter, Tue Jul 8 20:56:10 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(){
41  
42    const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
43 <  double vx2, vy2, vz2;
44 <  double kinetic, v_sqr;
45 <  int kl;
46 <  double jx2, jy2, jz2; // the square of the angular momentums
43 >  double kinetic;
44 >  double amass;
45 >  double aVel[3], aJ[3], I[3][3];
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++ ){
53  
54 <    vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
55 <    vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
56 <    vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();
54 >  for (kl=0; kl<integrableObjects.size(); kl++) {
55 >    integrableObjects[kl]->getVel(aVel);
56 >    amass = integrableObjects[kl]->getMass();
57  
58 <    v_sqr = vx2 + vy2 + vz2;
59 <    kinetic += atoms[kl]->getMass() * v_sqr;
58 >   for(j=0; j<3; j++)
59 >      kinetic += amass*aVel[j]*aVel[j];
60  
61 <    if( atoms[kl]->isDirectional() ){
62 <            
63 <      dAtom = (DirectionalAtom *)atoms[kl];
64 <      
65 <      jx2 = dAtom->getJx() * dAtom->getJx();    
66 <      jy2 = dAtom->getJy() * dAtom->getJy();
67 <      jz2 = dAtom->getJz() * dAtom->getJz();
68 <      
69 <      kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
70 <        + (jz2 / dAtom->getIzz());
71 <    }
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 89 | 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 108 | Line 111 | double Thermo::getPotential(){
111    potential = potential_local;
112   #endif // is_mpi
113  
111 – #ifdef IS_MPI
112 –  /*
113 –  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
114 –  */
115 – #endif
116 –
114    return potential;
115   }
116  
# Line 127 | 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;
141 <  double press[9];
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);
146 –  p = (press[0] + press[4] + press[8]) / 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 <  double volume;
157 <  double Hmat[9];
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 <  entry_plug->getBoxM(Hmat);
162 >  this->getPressureTensor(press);
163  
164 <  // volume = h1 (dot) h2 (cross) h3
164 >  pressureX = p_convert * press[0][0];
165  
166 <  volume = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
163 <         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
164 <         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );
165 <
166 <  return volume;
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[9];
175 <  double pressure;
174 >  double press[3][3];
175 >  double pressureY;
176  
177    this->getPressureTensor(press);
178  
179 <  pressure = p_convert * (press[0] + press[4] + press[8]) / 3.0;
179 >  pressureY = p_convert * press[1][1];
180  
181 <  return pressure;
181 >  return pressureY;
182   }
183  
184 + double Thermo::getPressureZ() {
185  
186 < void Thermo::getPressureTensor(double press[9]){
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 >
200 > void Thermo::getPressureTensor(double press[3][3]){
201    // returns pressure tensor in units amu*fs^-2*Ang^-1
202    // routine derived via viral theorem description in:
203    // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
# Line 192 | Line 207 | void Thermo::getPressureTensor(double press[9]){
207    double molmass, volume;
208    double vcom[3];
209    double p_local[9], p_global[9];
210 <  double theBox[3];
196 <  //double* tau;
197 <  int i, nMols;
198 <  Molecule* molecules;
210 >  int i, j, k;
211  
200 –  nMols = entry_plug->n_mol;
201 –  molecules = entry_plug->molecules;
202 –  //tau = entry_plug->tau;
203 –
204 –  // 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++) {
211 <    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 219 | Line 231 | void Thermo::getPressureTensor(double press[9]){
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 231 | Line 244 | void Thermo::getPressureTensor(double press[9]){
244    }
245   #endif // is_mpi
246  
247 <  volume = entry_plug->boxVol;
247 >  volume = this->getVolume();
248  
249 <  for(i=0; i<9; i++) {
250 <    press[i] = (p_global[i] - entry_plug->tau[i]*e_convert) / volume;
249 >
250 >
251 >  for(i = 0; i < 3; i++) {
252 >    for (j = 0; j < 3; j++) {
253 >      k = 3*i + j;
254 >      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
255 >    }
256    }
257   }
258  
259   void Thermo::velocitize() {
260    
261 <  double x,y;
262 <  double vx, vy, vz;
245 <  double jx, jy, jz;
246 <  int i, vr, vd; // velocity randomizer loop counters
261 >  double aVel[3], aJ[3], I[3][3];
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;
252 –  int n_atoms;
253 –  Atom** atoms;
254 –  DirectionalAtom* dAtom;
268    double temperature;
269 <  int n_oriented;
257 <  int n_constraints;
269 >  int nobj;
270  
271 <  atoms         = entry_plug->atoms;
260 <  n_atoms       = entry_plug->n_atoms;
261 <  temperature   = entry_plug->target_temp;
262 <  n_oriented    = entry_plug->n_oriented;
263 <  n_constraints = entry_plug->n_constraints;
271 >  nobj = info->integrableObjects.size();
272    
273 <  kebar = kb * temperature * (double)entry_plug->ndf /
266 <    ( 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 <
275 < //     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
276 <    
284 >
285      // picks random velocities from a gaussian distribution
286      // centered on vbar
287  
288 <    vx = vbar * gaussStream->getGaussian();
289 <    vy = vbar * gaussStream->getGaussian();
290 <    vz = vbar * gaussStream->getGaussian();
288 >    for (j=0; j<3; j++)
289 >      aVel[j] = vbar * gaussStream->getGaussian();
290 >    
291 >    info->integrableObjects[vr]->setVel( aVel );
292 >    
293 >    if(info->integrableObjects[vr]->isDirectional()){
294  
295 <    atoms[vr]->set_vx( vx );
296 <    atoms[vr]->set_vy( vy );
297 <    atoms[vr]->set_vz( vz );
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 293 | 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 <    vx = atoms[vd]->get_vx();
299 <    vy = atoms[vd]->get_vy();
300 <    vz = atoms[vd]->get_vz();
301 <        
302 <    vx -= vdrift[0];
303 <    vy -= vdrift[1];
304 <    vz -= vdrift[2];
331 >    info->integrableObjects[vd]->getVel(aVel);
332      
333 <    atoms[vd]->set_vx(vx);
334 <    atoms[vd]->set_vy(vy);
308 <    atoms[vd]->set_vz(vz);
309 <  }
310 <  if( n_oriented ){
311 <  
312 <    for( i=0; i<n_atoms; i++ ){
313 <      
314 <      if( atoms[i]->isDirectional() ){
315 <        
316 <        dAtom = (DirectionalAtom *)atoms[i];
317 <
318 <        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
319 <        jx = vbar * gaussStream->getGaussian();
320 <
321 <        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
322 <        jy = vbar * gaussStream->getGaussian();
333 >    for (j=0; j < 3; j++)
334 >      aVel[j] -= vdrift[j];
335          
336 <        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
325 <        jz = vbar * gaussStream->getGaussian();
326 <        
327 <        dAtom->setJx( jx );
328 <        dAtom->setJy( jy );
329 <        dAtom->setJz( jz );
330 <      }
331 <    }  
336 >    info->integrableObjects[vd]->setVel( aVel );
337    }
338 +
339   }
340  
341   void Thermo::getCOMVel(double vdrift[3]){
342  
343    double mtot, mtot_local;
344 +  double aVel[3], amass;
345    double vdrift_local[3];
346 <  int vd, n_atoms;
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
343 <  // We should really fix this before someone pokes an eye out.
349 >  nobj   = info->integrableObjects.size();
350  
345 –  n_atoms = entry_plug->n_atoms;  
346 –  atoms   = entry_plug->atoms;
347 –
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 <    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
359 <    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
360 <    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
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;
363      
364 <    mtot_local += atoms[vd]->getMass();
364 >    mtot_local += amass;
365    }
366  
367   #ifdef IS_MPI
# Line 375 | 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