OOPSE-1.0/libmdtools/Thermo.cpp

#include <math.h>
#include <iostream>
using namespace std;

#ifdef IS_MPI
#include <mpi.h>
#endif //is_mpi

#include "Thermo.hpp"
#include "SRI.hpp"
#include "Integrator.hpp"
#include "simError.h"
#include "MatVec3.h"

#ifdef IS_MPI
#define __C
#include "mpiSimulation.hpp"
#endif // is_mpi

inline double roundMe( double x ){
          return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
}

Thermo::Thermo( SimInfo* the_info ) { 
  info = the_info;
  int baseSeed = the_info->getSeed();
  
  gaussStream = new gaussianSPRNG( baseSeed );
}

Thermo::~Thermo(){
  delete gaussStream;
}

double Thermo::getKinetic(){

  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
  double kinetic;
  double amass;
  double aVel[3], aJ[3], I[3][3];
  int i, j, k, kl;

  double kinetic_global;
  vector<StuntDouble *> integrableObjects = info->integrableObjects;
  
  kinetic = 0.0;
  kinetic_global = 0.0;

  for (kl=0; kl<integrableObjects.size(); kl++) {
    integrableObjects[kl]->getVel(aVel);
    amass = integrableObjects[kl]->getMass();

   for(j=0; j<3; j++) 
      kinetic += amass*aVel[j]*aVel[j];

   if (integrableObjects[kl]->isDirectional()){
 
      integrableObjects[kl]->getJ( aJ );
      integrableObjects[kl]->getI( I );

      if (integrableObjects[kl]->isLinear()) {
        i = integrableObjects[kl]->linearAxis();
        j = (i+1)%3;
        k = (i+2)%3;
        kinetic += aJ[j]*aJ[j]/I[j][j] + aJ[k]*aJ[k]/I[k][k];
      } else {
        for (j=0; j<3; j++) 
          kinetic += aJ[j]*aJ[j] / I[j][j];
      }
   }
  }
#ifdef IS_MPI
  MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
                MPI_SUM, MPI_COMM_WORLD);
  kinetic = kinetic_global;
#endif //is_mpi
  
  kinetic = kinetic * 0.5 / e_convert;

  return kinetic;
}

double Thermo::getPotential(){
  
  double potential_local;
  double potential;
  int el, nSRI;
  Molecule* molecules;

  molecules = info->molecules;
  nSRI = info->n_SRI;

  potential_local = 0.0;
  potential = 0.0;
  potential_local += info->lrPot;

  for( el=0; el<info->n_mol; el++ ){    
    potential_local += molecules[el].getPotential();
  }

  // Get total potential for entire system from MPI.
#ifdef IS_MPI
  MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
                MPI_SUM, MPI_COMM_WORLD);
#else
  potential = potential_local; 
#endif // is_mpi

  return potential;
}

double Thermo::getTotalE(){

  double total;

  total = this->getKinetic() + this->getPotential();
  return total;
}

double Thermo::getTemperature(){

  const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
  double temperature;

  temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
  return temperature;
}

double Thermo::getVolume() {

  return info->boxVol;
}

double Thermo::getPressure() {

  // Relies on the calculation of the full molecular pressure tensor
  
  const double p_convert = 1.63882576e8;
  double press[3][3];
  double pressure;

  this->getPressureTensor(press);

  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;

  return pressure;
}

double Thermo::getPressureX() {

  // Relies on the calculation of the full molecular pressure tensor
  
  const double p_convert = 1.63882576e8;
  double press[3][3];
  double pressureX;

  this->getPressureTensor(press);

  pressureX = p_convert * press[0][0];

  return pressureX;
}

double Thermo::getPressureY() {

  // Relies on the calculation of the full molecular pressure tensor
  
  const double p_convert = 1.63882576e8;
  double press[3][3];
  double pressureY;

  this->getPressureTensor(press);

  pressureY = p_convert * press[1][1];

  return pressureY;
}

double Thermo::getPressureZ() {

  // Relies on the calculation of the full molecular pressure tensor
  
  const double p_convert = 1.63882576e8;
  double press[3][3];
  double pressureZ;

  this->getPressureTensor(press);

  pressureZ = p_convert * press[2][2];

  return pressureZ;
}


void Thermo::getPressureTensor(double press[3][3]){
  // returns pressure tensor in units amu*fs^-2*Ang^-1
  // routine derived via viral theorem description in:
  // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322

  const double e_convert = 4.184e-4;

  double molmass, volume;
  double vcom[3];
  double p_local[9], p_global[9];
  int i, j, k;

  for (i=0; i < 9; i++) {    
    p_local[i] = 0.0;
    p_global[i] = 0.0;
  }

  // use velocities of integrableObjects and their masses:  

  for (i=0; i < info->integrableObjects.size(); i++) {

    molmass = info->integrableObjects[i]->getMass();
    
    info->integrableObjects[i]->getVel(vcom);
    
    p_local[0] += molmass * (vcom[0] * vcom[0]); 
    p_local[1] += molmass * (vcom[0] * vcom[1]); 
    p_local[2] += molmass * (vcom[0] * vcom[2]); 
    p_local[3] += molmass * (vcom[1] * vcom[0]); 
    p_local[4] += molmass * (vcom[1] * vcom[1]); 
    p_local[5] += molmass * (vcom[1] * vcom[2]); 
    p_local[6] += molmass * (vcom[2] * vcom[0]); 
    p_local[7] += molmass * (vcom[2] * vcom[1]); 
    p_local[8] += molmass * (vcom[2] * vcom[2]); 

  }

  // Get total for entire system from MPI.
  
#ifdef IS_MPI
  MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#else
  for (i=0; i<9; i++) {
    p_global[i] = p_local[i]; 
  }
#endif // is_mpi

  volume = this->getVolume();


  for(i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      k = 3*i + j;
      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
    }
  }
}

void Thermo::velocitize() {
  
  double aVel[3], aJ[3], I[3][3];
  int i, j, l, m, n, vr, vd; // velocity randomizer loop counters
  double vdrift[3];
  double vbar;
  const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
  double av2;
  double kebar;
  double temperature;
  int nobj;

  if (!info->have_target_temp) {
    sprintf( painCave.errMsg,
             "You can't resample the velocities without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    painCave.severity = OOPSE_ERROR;
    simError();
    return;
  }

  nobj = info->integrableObjects.size();
  
  temperature   = info->target_temp;
  
  kebar = kb * temperature * (double)info->ndfRaw / 
    ( 2.0 * (double)info->ndf );
  
  for(vr = 0; vr < nobj; vr++){
    
    // uses equipartition theory to solve for vbar in angstrom/fs

    av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
    vbar = sqrt( av2 );

    // picks random velocities from a gaussian distribution
    // centered on vbar

    for (j=0; j<3; j++) 
      aVel[j] = vbar * gaussStream->getGaussian();
    
    info->integrableObjects[vr]->setVel( aVel );
    
    if(info->integrableObjects[vr]->isDirectional()){

      info->integrableObjects[vr]->getI( I );

      if (info->integrableObjects[vr]->isLinear()) {

        l= info->integrableObjects[vr]->linearAxis();
        m = (l+1)%3;
        n = (l+2)%3;

        aJ[l] = 0.0;
        vbar = sqrt( 2.0 * kebar * I[m][m] );
        aJ[m] = vbar * gaussStream->getGaussian();
        vbar = sqrt( 2.0 * kebar * I[n][n] );
        aJ[n] = vbar * gaussStream->getGaussian();
        
      } else {
        for (j = 0 ; j < 3; j++) {
          vbar = sqrt( 2.0 * kebar * I[j][j] );
          aJ[j] = vbar * gaussStream->getGaussian();
        }       
      } // else isLinear

      info->integrableObjects[vr]->setJ( aJ ); 
      
    }//isDirectional 

  }

  // Get the Center of Mass drift velocity.

  getCOMVel(vdrift);
  
  //  Corrects for the center of mass drift.
  // sums all the momentum and divides by total mass.

  for(vd = 0; vd < nobj; vd++){
    
    info->integrableObjects[vd]->getVel(aVel);
    
    for (j=0; j < 3; j++) 
      aVel[j] -= vdrift[j];
        
    info->integrableObjects[vd]->setVel( aVel );
  }

}

void Thermo::getCOMVel(double vdrift[3]){

  double mtot, mtot_local;
  double aVel[3], amass;
  double vdrift_local[3];
  int vd, j;
  int nobj;

  nobj   = info->integrableObjects.size();

  mtot_local = 0.0;
  vdrift_local[0] = 0.0;
  vdrift_local[1] = 0.0;
  vdrift_local[2] = 0.0;
  
  for(vd = 0; vd < nobj; vd++){
    
    amass = info->integrableObjects[vd]->getMass();
    info->integrableObjects[vd]->getVel( aVel );

    for(j = 0; j < 3; j++) 
      vdrift_local[j] += aVel[j] * amass;
    
    mtot_local += amass;
  }

#ifdef IS_MPI
  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  mtot = mtot_local;
  for(vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift_local[vd];
  }
#endif
    
  for (vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift[vd] / mtot;
  }
  
}

void Thermo::getCOM(double COM[3]){

  double mtot, mtot_local;
  double aPos[3], amass;
  double COM_local[3];
  int i, j;
  int nobj;

  mtot_local = 0.0;
  COM_local[0] = 0.0;
  COM_local[1] = 0.0;
  COM_local[2] = 0.0;

  nobj = info->integrableObjects.size();
  for(i = 0; i < nobj; i++){
    
    amass = info->integrableObjects[i]->getMass();
    info->integrableObjects[i]->getPos( aPos );

    for(j = 0; j < 3; j++) 
      COM_local[j] += aPos[j] * amass;
    
    mtot_local += amass;
  }

#ifdef IS_MPI
  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  mtot = mtot_local;
  for(i = 0; i < 3; i++) {
    COM[i] = COM_local[i];
  }
#endif
    
  for (i = 0; i < 3; i++) {
    COM[i] = COM[i] / mtot;
  }
}

void Thermo::removeCOMdrift() {
  double vdrift[3], aVel[3];
  int vd, j, nobj;

  nobj = info->integrableObjects.size();

  // Get the Center of Mass drift velocity.

  getCOMVel(vdrift);
  
  //  Corrects for the center of mass drift.
  // sums all the momentum and divides by total mass.

  for(vd = 0; vd < nobj; vd++){
    
    info->integrableObjects[vd]->getVel(aVel);
    
    for (j=0; j < 3; j++) 
      aVel[j] -= vdrift[j];
        
    info->integrableObjects[vd]->setVel( aVel );
  }
}
Revision:	1415
Committed:	Mon Jul 26 17:50:57 2004 UTC (19 years, 11 months ago) by gezelter
File size:	9909 byte(s)
Log Message:	Fixing required keyword annoyances
#	Content
1	#include <math.h>
2	#include <iostream>
3	using namespace std;
4
5	#ifdef IS_MPI
6	#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
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();
27
28	gaussStream = new gaussianSPRNG( baseSeed );
29	}
30
31	Thermo::~Thermo(){
32	delete gaussStream;
33	}
34
35	double Thermo::getKinetic(){
36
37	const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
38	double kinetic;
39	double amass;
40	double aVel[3], aJ[3], I[3][3];
41	int i, j, k, kl;
42
43	double kinetic_global;
44	vector<StuntDouble *> integrableObjects = info->integrableObjects;
45
46	kinetic = 0.0;
47	kinetic_global = 0.0;
48
49	for (kl=0; kl<integrableObjects.size(); kl++) {
50	integrableObjects[kl]->getVel(aVel);
51	amass = integrableObjects[kl]->getMass();
52
53	for(j=0; j<3; j++)
54	kinetic += amassaVel[j]aVel[j];
55
56	if (integrableObjects[kl]->isDirectional()){
57
58	integrableObjects[kl]->getJ( aJ );
59	integrableObjects[kl]->getI( I );
60
61	if (integrableObjects[kl]->isLinear()) {
62	i = integrableObjects[kl]->linearAxis();
63	j = (i+1)%3;
64	k = (i+2)%3;
65	kinetic += aJ[j]aJ[j]/I[j][j] + aJ[k]aJ[k]/I[k][k];
66	} else {
67	for (j=0; j<3; j++)
68	kinetic += aJ[j]*aJ[j] / I[j][j];
69	}
70	}
71	}
72	#ifdef IS_MPI
73	MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
74	MPI_SUM, MPI_COMM_WORLD);
75	kinetic = kinetic_global;
76	#endif //is_mpi
77
78	kinetic = kinetic * 0.5 / e_convert;
79
80	return kinetic;
81	}
82
83	double Thermo::getPotential(){
84
85	double potential_local;
86	double potential;
87	int el, nSRI;
88	Molecule* molecules;
89
90	molecules = info->molecules;
91	nSRI = info->n_SRI;
92
93	potential_local = 0.0;
94	potential = 0.0;
95	potential_local += info->lrPot;
96
97	for( el=0; el<info->n_mol; el++ ){
98	potential_local += molecules[el].getPotential();
99	}
100
101	// Get total potential for entire system from MPI.
102	#ifdef IS_MPI
103	MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
104	MPI_SUM, MPI_COMM_WORLD);
105	#else
106	potential = potential_local;
107	#endif // is_mpi
108
109	return potential;
110	}
111
112	double Thermo::getTotalE(){
113
114	double total;
115
116	total = this->getKinetic() + this->getPotential();
117	return total;
118	}
119
120	double Thermo::getTemperature(){
121
122	const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
123	double temperature;
124
125	temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
126	return temperature;
127	}
128
129	double Thermo::getVolume() {
130
131	return info->boxVol;
132	}
133
134	double Thermo::getPressure() {
135
136	// Relies on the calculation of the full molecular pressure tensor
137
138	const double p_convert = 1.63882576e8;
139	double press[3][3];
140	double pressure;
141
142	this->getPressureTensor(press);
143
144	pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
145
146	return pressure;
147	}
148
149	double Thermo::getPressureX() {
150
151	// Relies on the calculation of the full molecular pressure tensor
152
153	const double p_convert = 1.63882576e8;
154	double press[3][3];
155	double pressureX;
156
157	this->getPressureTensor(press);
158
159	pressureX = p_convert * press[0][0];
160
161	return pressureX;
162	}
163
164	double Thermo::getPressureY() {
165
166	// Relies on the calculation of the full molecular pressure tensor
167
168	const double p_convert = 1.63882576e8;
169	double press[3][3];
170	double pressureY;
171
172	this->getPressureTensor(press);
173
174	pressureY = p_convert * press[1][1];
175
176	return pressureY;
177	}
178
179	double Thermo::getPressureZ() {
180
181	// Relies on the calculation of the full molecular pressure tensor
182
183	const double p_convert = 1.63882576e8;
184	double press[3][3];
185	double pressureZ;
186
187	this->getPressureTensor(press);
188
189	pressureZ = p_convert * press[2][2];
190
191	return pressureZ;
192	}
193
194
195	void Thermo::getPressureTensor(double press[3][3]){
196	// returns pressure tensor in units amufs^-2Ang^-1
197	// routine derived via viral theorem description in:
198	// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
199
200	const double e_convert = 4.184e-4;
201
202	double molmass, volume;
203	double vcom[3];
204	double p_local[9], p_global[9];
205	int i, j, k;
206
207	for (i=0; i < 9; i++) {
208	p_local[i] = 0.0;
209	p_global[i] = 0.0;
210	}
211
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]);
223	p_local[3] += molmass * (vcom[1] * vcom[0]);
224	p_local[4] += molmass * (vcom[1] * vcom[1]);
225	p_local[5] += molmass * (vcom[1] * vcom[2]);
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
234	#ifdef IS_MPI
235	MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
236	#else
237	for (i=0; i<9; i++) {
238	p_global[i] = p_local[i];
239	}
240	#endif // is_mpi
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;
250	}
251	}
252	}
253
254	void Thermo::velocitize() {
255
256	double aVel[3], aJ[3], I[3][3];
257	int i, j, l, m, n, vr, vd; // velocity randomizer loop counters
258	double vdrift[3];
259	double vbar;
260	const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
261	double av2;
262	double kebar;
263	double temperature;
264	int nobj;
265
266	if (!info->have_target_temp) {
267	sprintf( painCave.errMsg,
268	"You can't resample the velocities without a targetTemp!\n"
269	);
270	painCave.isFatal = 1;
271	painCave.severity = OOPSE_ERROR;
272	simError();
273	return;
274	}
275
276	nobj = info->integrableObjects.size();
277
278	temperature = info->target_temp;
279
280	kebar = kb * temperature * (double)info->ndfRaw /
281	( 2.0 * (double)info->ndf );
282
283	for(vr = 0; vr < nobj; vr++){
284
285	// uses equipartition theory to solve for vbar in angstrom/fs
286
287	av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
288	vbar = sqrt( av2 );
289
290	// picks random velocities from a gaussian distribution
291	// centered on vbar
292
293	for (j=0; j<3; j++)
294	aVel[j] = vbar * gaussStream->getGaussian();
295
296	info->integrableObjects[vr]->setVel( aVel );
297
298	if(info->integrableObjects[vr]->isDirectional()){
299
300	info->integrableObjects[vr]->getI( I );
301
302	if (info->integrableObjects[vr]->isLinear()) {
303
304	l= info->integrableObjects[vr]->linearAxis();
305	m = (l+1)%3;
306	n = (l+2)%3;
307
308	aJ[l] = 0.0;
309	vbar = sqrt( 2.0 * kebar * I[m][m] );
310	aJ[m] = vbar * gaussStream->getGaussian();
311	vbar = sqrt( 2.0 * kebar * I[n][n] );
312	aJ[n] = vbar * gaussStream->getGaussian();
313
314	} else {
315	for (j = 0 ; j < 3; j++) {
316	vbar = sqrt( 2.0 * kebar * I[j][j] );
317	aJ[j] = vbar * gaussStream->getGaussian();
318	}
319	} // else isLinear
320
321	info->integrableObjects[vr]->setJ( aJ );
322
323	}//isDirectional
324
325	}
326
327	// Get the Center of Mass drift velocity.
328
329	getCOMVel(vdrift);
330
331	// Corrects for the center of mass drift.
332	// sums all the momentum and divides by total mass.
333
334	for(vd = 0; vd < nobj; vd++){
335
336	info->integrableObjects[vd]->getVel(aVel);
337
338	for (j=0; j < 3; j++)
339	aVel[j] -= vdrift[j];
340
341	info->integrableObjects[vd]->setVel( aVel );
342	}
343
344	}
345
346	void Thermo::getCOMVel(double vdrift[3]){
347
348	double mtot, mtot_local;
349	double aVel[3], amass;
350	double vdrift_local[3];
351	int vd, j;
352	int nobj;
353
354	nobj = info->integrableObjects.size();
355
356	mtot_local = 0.0;
357	vdrift_local[0] = 0.0;
358	vdrift_local[1] = 0.0;
359	vdrift_local[2] = 0.0;
360
361	for(vd = 0; vd < nobj; vd++){
362
363	amass = info->integrableObjects[vd]->getMass();
364	info->integrableObjects[vd]->getVel( aVel );
365
366	for(j = 0; j < 3; j++)
367	vdrift_local[j] += aVel[j] * amass;
368
369	mtot_local += amass;
370	}
371
372	#ifdef IS_MPI
373	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
374	MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
375	#else
376	mtot = mtot_local;
377	for(vd = 0; vd < 3; vd++) {
378	vdrift[vd] = vdrift_local[vd];
379	}
380	#endif
381
382	for (vd = 0; vd < 3; vd++) {
383	vdrift[vd] = vdrift[vd] / mtot;
384	}
385
386	}
387
388	void Thermo::getCOM(double COM[3]){
389
390	double mtot, mtot_local;
391	double aPos[3], amass;
392	double COM_local[3];
393	int i, j;
394	int nobj;
395
396	mtot_local = 0.0;
397	COM_local[0] = 0.0;
398	COM_local[1] = 0.0;
399	COM_local[2] = 0.0;
400
401	nobj = info->integrableObjects.size();
402	for(i = 0; i < nobj; i++){
403
404	amass = info->integrableObjects[i]->getMass();
405	info->integrableObjects[i]->getPos( aPos );
406
407	for(j = 0; j < 3; j++)
408	COM_local[j] += aPos[j] * amass;
409
410	mtot_local += amass;
411	}
412
413	#ifdef IS_MPI
414	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
415	MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
416	#else
417	mtot = mtot_local;
418	for(i = 0; i < 3; i++) {
419	COM[i] = COM_local[i];
420	}
421	#endif
422
423	for (i = 0; i < 3; i++) {
424	COM[i] = COM[i] / mtot;
425	}
426	}
427
428	void Thermo::removeCOMdrift() {
429	double vdrift[3], aVel[3];
430	int vd, j, nobj;
431
432	nobj = info->integrableObjects.size();
433
434	// Get the Center of Mass drift velocity.
435
436	getCOMVel(vdrift);
437
438	// Corrects for the center of mass drift.
439	// sums all the momentum and divides by total mass.
440
441	for(vd = 0; vd < nobj; vd++){
442
443	info->integrableObjects[vd]->getVel(aVel);
444
445	for (j=0; j < 3; j++)
446	aVel[j] -= vdrift[j];
447
448	info->integrableObjects[vd]->setVel( aVel );
449	}
450	}