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

Comparing:
branches/mmeineke/OOPSE/libmdtools/Thermo.cpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/Thermo.cpp (file contents), 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  
5   #ifdef IS_MPI
6   #include <mpi.h>
7 – #include <mpi++.h>
7   #endif //is_mpi
8  
9   #include "Thermo.hpp"
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 < #define BASE_SEED 123456789
17 > #ifdef IS_MPI
18 > #define __C
19 > #include "mpiSimulation.hpp"
20 > #endif // is_mpi
21  
22 < Thermo::Thermo( SimInfo* the_entry_plug ) {
23 <  entry_plug = the_entry_plug;
24 <  int baseSeed = BASE_SEED;
22 > inline double roundMe( double x ){
23 >          return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
24 > }
25 >
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  
35 –  DirectionalAtom *dAtom;
36 –
37 –  int n_atoms;
48    double kinetic_global;
49 <  Atom** atoms;
40 <
49 >  vector<StuntDouble *> integrableObjects = info->integrableObjects;
50    
42 –  n_atoms = entry_plug->n_atoms;
43 –  atoms = entry_plug->atoms;
44 –
51    kinetic = 0.0;
52    kinetic_global = 0.0;
47 –  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::COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
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 77 | Line 87 | double Thermo::getPotential(){
87  
88   double Thermo::getPotential(){
89    
90 +  double potential_local;
91    double potential;
81 –  double potential_global;
92    int el, nSRI;
93 <  SRI** sris;
93 >  Molecule* molecules;
94  
95 <  sris = entry_plug->sr_interactions;
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_global = 0.0;
90 <  potential += entry_plug->lrPot;
100 >  potential_local += info->lrPot;
101  
102 <  for( el=0; el<nSRI; el++ ){
103 <    
94 <    potential += sris[el]->get_potential();
102 >  for( el=0; el<info->n_mol; el++ ){    
103 >    potential_local += molecules[el].getPotential();
104    }
105  
106    // Get total potential for entire system from MPI.
107   #ifdef IS_MPI
108 <  MPI::COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
109 <  potential = potential_global;
110 <
108 >  MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
109 >                MPI_SUM, MPI_COMM_WORLD);
110 > #else
111 >  potential = potential_local;
112   #endif // is_mpi
113  
114    return potential;
# Line 114 | 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;
119 –  
120 –  int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
121 –    - entry_plug->n_constraints - 3;
129  
130 <  temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
130 >  temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
131    return temperature;
132   }
133  
134 < double Thermo::getPressure(){
134 > double Thermo::getVolume() {
135  
136 < //  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
137 < // const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
131 < //  const double conv_A_m = 1.0E-10; //convert A -> m
136 >  return info->boxVol;
137 > }
138  
139 <  return 0.0;
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);
148 >
149 >  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
150 >
151 >  return pressure;
152   }
153  
154 + double Thermo::getPressureX() {
155 +
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::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 pressureY;
176 +
177 +  this->getPressureTensor(press);
178 +
179 +  pressureY = p_convert * press[1][1];
180 +
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 +
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
204 +
205 +  const double e_convert = 4.184e-4;
206 +
207 +  double molmass, volume;
208 +  double vcom[3];
209 +  double p_local[9], p_global[9];
210 +  int i, j, k;
211 +
212 +  for (i=0; i < 9; i++) {    
213 +    p_local[i] = 0.0;
214 +    p_global[i] = 0.0;
215 +  }
216 +
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]);
228 +    p_local[3] += molmass * (vcom[1] * vcom[0]);
229 +    p_local[4] += molmass * (vcom[1] * vcom[1]);
230 +    p_local[5] += molmass * (vcom[1] * vcom[2]);
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 +  
239 + #ifdef IS_MPI
240 +  MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
241 + #else
242 +  for (i=0; i<9; i++) {
243 +    p_global[i] = p_local[i];
244 +  }
245 + #endif // is_mpi
246 +
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 +      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;
140 <  double jx, jy, jz;
141 <  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];
143 –  double mtot = 0.0;
264    double vbar;
265    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
266    double av2;
267    double kebar;
148 –  int ndf; // number of degrees of freedom
149 –  int ndfRaw; // the raw number of degrees of freedom
150 –  int n_atoms;
151 –  Atom** atoms;
152 –  DirectionalAtom* dAtom;
268    double temperature;
269 <  int n_oriented;
155 <  int n_constraints;
269 >  int nobj;
270  
271 <  atoms         = entry_plug->atoms;
158 <  n_atoms       = entry_plug->n_atoms;
159 <  temperature   = entry_plug->target_temp;
160 <  n_oriented    = entry_plug->n_oriented;
161 <  n_constraints = entry_plug->n_constraints;
271 >  nobj = info->integrableObjects.size();
272    
273 <
164 <  ndfRaw = 3 * n_atoms + 3 * n_oriented;
165 <  ndf = ndfRaw - n_constraints - 3;
166 <  kebar = kb * temperature * (double)ndf / ( 2.0 * (double)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  
175 – //     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
176 –    
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.
323 +
324 +  getCOMVel(vdrift);
325    
326    //  Corrects for the center of mass drift.
327    // sums all the momentum and divides by total mass.
328 <  
329 <  mtot = 0.0;
193 <  vdrift[0] = 0.0;
194 <  vdrift[1] = 0.0;
195 <  vdrift[2] = 0.0;
196 <  for(vd = 0; vd < n_atoms; vd++){
328 >
329 >  for(vd = 0; vd < nobj; vd++){
330      
331 <    vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
199 <    vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
200 <    vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
331 >    info->integrableObjects[vd]->getVel(aVel);
332      
333 <    mtot += atoms[vd]->getMass();
333 >    for (j=0; j < 3; j++)
334 >      aVel[j] -= vdrift[j];
335 >        
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, j;
347 +  int nobj;
348 +
349 +  nobj   = info->integrableObjects.size();
350 +
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 < nobj; vd++){
357 +    
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 += amass;
365 +  }
366 +
367 + #ifdef IS_MPI
368 +  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
369 +  MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
370 + #else
371 +  mtot = mtot_local;
372 +  for(vd = 0; vd < 3; vd++) {
373 +    vdrift[vd] = vdrift_local[vd];
374 +  }
375 + #endif
376 +    
377    for (vd = 0; vd < 3; vd++) {
378      vdrift[vd] = vdrift[vd] / mtot;
379    }
380    
381 + }
382  
383 <  for(vd = 0; vd < n_atoms; vd++){
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 <    vx = atoms[vd]->get_vx();
400 <    vy = atoms[vd]->get_vy();
401 <    vz = atoms[vd]->get_vz();
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 <    vx -= vdrift[0];
419 <    vy -= vdrift[1];
420 <    vz -= vdrift[2];
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 <    atoms[vd]->set_vx(vx);
439 <    atoms[vd]->set_vy(vy);
440 <    atoms[vd]->set_vz(vz);
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 <  if( n_oriented ){
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 <    for( i=0; i<n_atoms; i++ ){
471 <      
472 <      if( atoms[i]->isDirectional() ){
473 <        
474 <        dAtom = (DirectionalAtom *)atoms[i];
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 <        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
502 <        jx = vbar * gaussStream->getGaussian();
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 <        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
507 <        jy = vbar * gaussStream->getGaussian();
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 <        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
512 <        jz = vbar * gaussStream->getGaussian();
241 <        
242 <        dAtom->setJx( jx );
243 <        dAtom->setJy( jy );
244 <        dAtom->setJz( jz );
245 <      }
246 <    }  
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