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 1127 by tim, Tue Apr 20 16:56:40 2004 UTC vs.
Revision 1253 by gezelter, Tue Jun 8 16:49:46 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  
15   #ifdef IS_MPI
16   #define __C
17   #include "mpiSimulation.hpp"
18   #endif // is_mpi
19  
20 + inline double roundMe( double x ){
21 +          return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
22 + }
23 +
24   Thermo::Thermo( SimInfo* the_info ) {
25    info = the_info;
26    int baseSeed = the_info->getSeed();
# Line 197 | Line 202 | void Thermo::getPressureTensor(double press[3][3]){
202    double molmass, volume;
203    double vcom[3];
204    double p_local[9], p_global[9];
205 <  int i, j, k, nMols;
201 <  Molecule* molecules;
205 >  int i, j, k;
206  
203  nMols = info->n_mol;
204  molecules = info->molecules;
205  //tau = info->tau;
206
207  // use velocities of molecular centers of mass and molecular masses:
207    for (i=0; i < 9; i++) {    
208      p_local[i] = 0.0;
209      p_global[i] = 0.0;
210    }
211  
212 <  for (i=0; i < nMols; i++) {
214 <    molmass = molecules[i].getCOMvel(vcom);
212 >  // use velocities of integrableObjects and their masses:  
213  
214 +  for (i=0; i < info->integrableObjects.size(); i++) {
215 +
216 +    molmass = info->integrableObjects[i]->getMass();
217 +    
218 +    info->integrableObjects[i]->getVel(vcom);
219 +    
220      p_local[0] += molmass * (vcom[0] * vcom[0]);
221      p_local[1] += molmass * (vcom[0] * vcom[1]);
222      p_local[2] += molmass * (vcom[0] * vcom[2]);
# Line 222 | Line 226 | void Thermo::getPressureTensor(double press[3][3]){
226      p_local[6] += molmass * (vcom[2] * vcom[0]);
227      p_local[7] += molmass * (vcom[2] * vcom[1]);
228      p_local[8] += molmass * (vcom[2] * vcom[2]);
229 +
230    }
231  
232    // Get total for entire system from MPI.
233 <
233 >  
234   #ifdef IS_MPI
235    MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
236   #else
# Line 236 | Line 241 | void Thermo::getPressureTensor(double press[3][3]){
241  
242    volume = this->getVolume();
243  
244 +
245 +
246    for(i = 0; i < 3; i++) {
247      for (j = 0; j < 3; j++) {
248        k = 3*i + j;
249        press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
243
250      }
251    }
252   }
# Line 431 | Line 437 | void Thermo::removeCOMdrift() {
437          
438      info->integrableObjects[vd]->setVel( aVel );
439    }
440 < }
440 > }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines