OOPSE/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"

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

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, nMols;
  Molecule* molecules;

  nMols = info->n_mol;
  molecules = info->molecules;
  //tau = info->tau;

  // use velocities of molecular centers of mass and molecular masses:
  for (i=0; i < 9; i++) {    
    p_local[i] = 0.0;
    p_global[i] = 0.0;
  }

  for (i=0; i < nMols; i++) {
    molmass = molecules[i].getCOMvel(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;

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