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 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# 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   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 33 | 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  
38 –  DirectionalAtom *dAtom;
39 –
40 –  int n_atoms;
48    double kinetic_global;
49 <  Atom** atoms;
43 <
49 >  vector<StuntDouble *> integrableObjects = info->integrableObjects;
50    
45 –  n_atoms = info->n_atoms;
46 –  atoms = info->atoms;
47 –
51    kinetic = 0.0;
52    kinetic_global = 0.0;
50 –  for( kl=0; kl < n_atoms; kl++ ){
51 –    
52 –    atoms[kl]->getVel(aVel);
53 –    amass = atoms[kl]->getMass();
54 –    
55 –    for (j=0; j < 3; j++)
56 –      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 104 | Line 111 | double Thermo::getPotential(){
111    potential = potential_local;
112   #endif // is_mpi
113  
107 – #ifdef IS_MPI
108 –  /*
109 –  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
110 –  */
111 – #endif
112 –
114    return potential;
115   }
116  
# Line 125 | Line 126 | double Thermo::getTemperature(){
126  
127    const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
128    double temperature;
129 <  
129 >
130    temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
131    return temperature;
132   }
# Line 206 | 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, nMols;
210 <  Molecule* molecules;
210 >  int i, j, k;
211  
212 –  nMols = info->n_mol;
213 –  molecules = info->molecules;
214 –  //tau = info->tau;
215 –
216 –  // 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++) {
223 <    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 231 | 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 245 | Line 246 | void Thermo::getPressureTensor(double press[3][3]){
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;
252 –
255      }
256    }
257   }
# Line 257 | Line 259 | void Thermo::velocitize() {
259   void Thermo::velocitize() {
260    
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;
266 –  int n_atoms;
267 –  Atom** atoms;
268 –  DirectionalAtom* dAtom;
268    double temperature;
269 <  int n_oriented;
271 <  int n_constraints;
269 >  int nobj;
270  
271 <  atoms         = info->atoms;
272 <  n_atoms       = info->n_atoms;
271 >  nobj = info->integrableObjects.size();
272 >  
273    temperature   = info->target_temp;
276 –  n_oriented    = info->n_oriented;
277 –  n_constraints = info->n_constraints;
274    
275    kebar = kb * temperature * (double)info->ndfRaw /
276      ( 2.0 * (double)info->ndf );
277    
278 <  for(vr = 0; vr < n_atoms; vr++){
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 <
289 < //     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
290 <    
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 305 | 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    }
317 –  if( n_oriented ){
318 –  
319 –    for( i=0; i<n_atoms; i++ ){
320 –      
321 –      if( atoms[i]->isDirectional() ){
322 –        
323 –        dAtom = (DirectionalAtom *)atoms[i];
324 –        dAtom->getI( I );
325 –        
326 –        for (j = 0 ; j < 3; j++) {
338  
328 –          vbar = sqrt( 2.0 * kebar * I[j][j] );
329 –          aJ[j] = vbar * gaussStream->getGaussian();
330 –
331 –        }      
332 –
333 –        dAtom->setJ( aJ );
334 –
335 –      }
336 –    }  
337 –  }
339   }
340  
341   void Thermo::getCOMVel(double vdrift[3]){
# Line 342 | 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
349 <  // We should really fix this before someone pokes an eye out.
349 >  nobj   = info->integrableObjects.size();
350  
351 –  n_atoms = info->n_atoms;  
352 –  atoms   = info->atoms;
353 –
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 388 | Line 385 | void Thermo::getCOM(double COM[3]){
385    double mtot, mtot_local;
386    double aPos[3], amass;
387    double COM_local[3];
388 <  int i, n_atoms, j;
389 <  Atom** atoms;
393 <
394 <  // We are very careless here with the distinction between n_atoms and n_local
395 <  // We should really fix this before someone pokes an eye out.
396 <
397 <  n_atoms = info->n_atoms;  
398 <  atoms   = info->atoms;
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 <  for(i = 0; i < n_atoms; i++){
395 >
396 >  nobj = info->integrableObjects.size();
397 >  for(i = 0; i < nobj; i++){
398      
399 <    amass = atoms[i]->getMass();
400 <    atoms[i]->getPos( aPos );
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;
# Line 427 | Line 419 | void Thermo::getCOM(double COM[3]){
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