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 1133 by gezelter, Mon Apr 26 14:29:18 2004 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# Line 11 | Line 11 | using namespace std;
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
# Line 26 | Line 28 | Thermo::Thermo( SimInfo* 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 200 | Line 205 | void Thermo::getPressureTensor(double press[3][3]){
205    const double e_convert = 4.184e-4;
206  
207    double molmass, volume;
208 <  double vcom[3], pcom[3], fcom[3], scaled[3];
208 >  double vcom[3];
209    double p_local[9], p_global[9];
210 <  int i, j, k, nMols;
206 <  Molecule* molecules;
210 >  int i, j, k;
211  
208 –  nMols = info->n_mol;
209 –  molecules = info->molecules;
210 –  //tau = info->tau;
211 –
212 –  // 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 +  // 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);
223 –    info->integrableObjects[i]->getPos(pcom);
224 –    info->integrableObjects[i]->getFrc(fcom);
225 –
226 –    matVecMul3(info->HmatInv, pcom, scaled);
227 –  
228 –    for(j=0; j<3; j++)
229 –      scaled[j] -= roundMe(scaled[j]);
230 –
231 –    // calc the wrapped real coordinates from the wrapped scaled coordinates
232 –  
233 –    matVecMul3(info->Hmat, scaled, pcom);
224      
225 <    p_local[0] += molmass * (vcom[0] * vcom[0]) + fcom[0]*pcom[0]*eConvert;
226 <    p_local[1] += molmass * (vcom[0] * vcom[1]) + fcom[0]*pcom[1]*eConvert;
227 <    p_local[2] += molmass * (vcom[0] * vcom[2]) + fcom[0]*pcom[2]*eConvert;
228 <    p_local[3] += molmass * (vcom[1] * vcom[0]) + fcom[1]*pcom[0]*eConvert;
229 <    p_local[4] += molmass * (vcom[1] * vcom[1]) + fcom[1]*pcom[1]*eConvert;
230 <    p_local[5] += molmass * (vcom[1] * vcom[2]) + fcom[1]*pcom[2]*eConvert;
231 <    p_local[6] += molmass * (vcom[2] * vcom[0]) + fcom[2]*pcom[0]*eConvert;
232 <    p_local[7] += molmass * (vcom[2] * vcom[1]) + fcom[2]*pcom[1]*eConvert;
233 <    p_local[8] += molmass * (vcom[2] * vcom[2]) + fcom[2]*pcom[2]*eConvert;
234 <    
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 <
238 >  
239   #ifdef IS_MPI
240    MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
241   #else
# Line 255 | Line 245 | void Thermo::getPressureTensor(double press[3][3]){
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] /  volume;
263 <
254 >      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
255      }
256    }
257   }
# Line 452 | Line 443 | void Thermo::removeCOMdrift() {
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