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, 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;
  int n_atoms;
  Atom** atoms;
  DirectionalAtom* dAtom;
  double temperature;
  int n_oriented;
  int n_constraints;

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

    av2 = 2.0 * kebar / atoms[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();
    
    atoms[vr]->setVel( aVel );

  }

  // 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 < n_atoms; vd++){
    
    atoms[vd]->getVel(aVel);
    
    for (j=0; j < 3; j++) 
      aVel[j] -= vdrift[j];
        
    atoms[vd]->setVel( aVel );
  }
  if( n_oriented ){
  
    for( i=0; i<n_atoms; i++ ){
      
      if( atoms[i]->isDirectional() ){
        
        dAtom = (DirectionalAtom *)atoms[i];
        dAtom->getI( I );
        
        for (j = 0 ; j < 3; j++) {

          vbar = sqrt( 2.0 * kebar * I[j][j] );
          aJ[j] = vbar * gaussStream->getGaussian();

        }       

        dAtom->setJ( aJ );

      }
    }   
  }
}

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

  double mtot, mtot_local;
  double aVel[3], amass;
  double vdrift_local[3];
  int vd, n_atoms, j;
  Atom** atoms;

  // We are very careless here with the distinction between n_atoms and n_local
  // We should really fix this before someone pokes an eye out.

  n_atoms = info->n_atoms;  
  atoms   = info->atoms;

  mtot_local = 0.0;
  vdrift_local[0] = 0.0;
  vdrift_local[1] = 0.0;
  vdrift_local[2] = 0.0;
  
  for(vd = 0; vd < n_atoms; vd++){
    
    amass = atoms[vd]->getMass();
    atoms[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, n_atoms, j;
  Atom** atoms;

  // We are very careless here with the distinction between n_atoms and n_local
  // We should really fix this before someone pokes an eye out.

  n_atoms = info->n_atoms;  
  atoms   = info->atoms;

  mtot_local = 0.0;
  COM_local[0] = 0.0;
  COM_local[1] = 0.0;
  COM_local[2] = 0.0;
  
  for(i = 0; i < n_atoms; i++){
    
    amass = atoms[i]->getMass();
    atoms[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;
  }
}
Revision:	1125
Committed:	Mon Apr 19 22:13:01 2004 UTC (20 years, 2 months ago) by gezelter
File size:	9069 byte(s)
Log Message:	Fixed a charge bug
#	User	Rev	Content
1	gezelter	829	#include <math.h>
2	mmeineke	377	#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	chuckv	438	#include "simError.h"
13	mmeineke	402
14			#ifdef IS_MPI
15	chuckv	401	#define __C
16	mmeineke	402	#include "mpiSimulation.hpp"
17			#endif // is_mpi
18	mmeineke	377
19	mmeineke	614	Thermo::Thermo( SimInfo* the_info ) {
20			info = the_info;
21	tim	708	int baseSeed = the_info->getSeed();
22	mmeineke	377
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	gezelter	608	double kinetic;
34			double amass;
35			double aVel[3], aJ[3], I[3][3];
36	tim	1118	int i, j, k, kl;
37	mmeineke	377
38			double kinetic_global;
39	tim	1113	vector<StuntDouble *> integrableObjects = info->integrableObjects;
40	mmeineke	377
41			kinetic = 0.0;
42			kinetic_global = 0.0;
43
44	tim	1113	for (kl=0; kl<integrableObjects.size(); kl++) {
45			integrableObjects[kl]->getVel(aVel);
46			amass = integrableObjects[kl]->getMass();
47	gezelter	608
48	tim	1113	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	tim	1118	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	gezelter	1125	for (j=0; j<3; j++)
63			kinetic += aJ[j]*aJ[j] / I[j][j];
64	tim	1118	}
65	gezelter	1125	}
66	mmeineke	377	}
67			#ifdef IS_MPI
68	mmeineke	447	MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
69			MPI_SUM, MPI_COMM_WORLD);
70	mmeineke	377	kinetic = kinetic_global;
71			#endif //is_mpi
72	gezelter	1125
73	mmeineke	377	kinetic = kinetic * 0.5 / e_convert;
74
75			return kinetic;
76			}
77
78			double Thermo::getPotential(){
79
80	chuckv	401	double potential_local;
81	mmeineke	377	double potential;
82			int el, nSRI;
83	mmeineke	428	Molecule* molecules;
84	mmeineke	377
85	mmeineke	614	molecules = info->molecules;
86			nSRI = info->n_SRI;
87	mmeineke	377
88	chuckv	401	potential_local = 0.0;
89	chuckv	438	potential = 0.0;
90	mmeineke	614	potential_local += info->lrPot;
91	mmeineke	377
92	mmeineke	614	for( el=0; el<info->n_mol; el++ ){
93	mmeineke	428	potential_local += molecules[el].getPotential();
94	mmeineke	377	}
95
96			// Get total potential for entire system from MPI.
97			#ifdef IS_MPI
98	mmeineke	447	MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
99			MPI_SUM, MPI_COMM_WORLD);
100	chuckv	401	#else
101			potential = potential_local;
102	mmeineke	377	#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	gezelter	454	double Thermo::getTemperature(){
116
117	tim	763	const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
118	gezelter	454	double temperature;
119	tim	1113
120	mmeineke	614	temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
121	mmeineke	377	return temperature;
122			}
123
124	gezelter	484	double Thermo::getVolume() {
125	gezelter	574
126	mmeineke	614	return info->boxVol;
127	gezelter	484	}
128
129	gezelter	483	double Thermo::getPressure() {
130	gezelter	574
131	gezelter	483	// Relies on the calculation of the full molecular pressure tensor
132
133			const double p_convert = 1.63882576e8;
134	gezelter	588	double press[3][3];
135	gezelter	483	double pressure;
136
137			this->getPressureTensor(press);
138
139	gezelter	588	pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
140	gezelter	483
141			return pressure;
142			}
143
144	mmeineke	755	double Thermo::getPressureX() {
145	gezelter	483
146	mmeineke	755	// 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	gezelter	588	void Thermo::getPressureTensor(double press[3][3]){
191	gezelter	483	// returns pressure tensor in units amufs^-2Ang^-1
192	gezelter	445	// routine derived via viral theorem description in:
193			// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
194	mmeineke	377
195	gezelter	477	const double e_convert = 4.184e-4;
196	gezelter	483
197			double molmass, volume;
198	gezelter	468	double vcom[3];
199	gezelter	483	double p_local[9], p_global[9];
200	mmeineke	614	int i, j, k, nMols;
201	gezelter	468	Molecule* molecules;
202
203	mmeineke	614	nMols = info->n_mol;
204			molecules = info->molecules;
205			//tau = info->tau;
206	gezelter	468
207			// use velocities of molecular centers of mass and molecular masses:
208	gezelter	483	for (i=0; i < 9; i++) {
209			p_local[i] = 0.0;
210			p_global[i] = 0.0;
211			}
212	gezelter	475
213	gezelter	468	for (i=0; i < nMols; i++) {
214	gezelter	475	molmass = molecules[i].getCOMvel(vcom);
215	gezelter	483
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	gezelter	468	}
226
227			// Get total for entire system from MPI.
228	chuckv	479
229	gezelter	468	#ifdef IS_MPI
230	gezelter	483	MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
231	gezelter	468	#else
232	gezelter	483	for (i=0; i<9; i++) {
233			p_global[i] = p_local[i];
234			}
235	gezelter	468	#endif // is_mpi
236
237	gezelter	611	volume = this->getVolume();
238	gezelter	468
239	gezelter	588	for(i = 0; i < 3; i++) {
240			for (j = 0; j < 3; j++) {
241			k = 3*i + j;
242	mmeineke	614	press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
243
244	gezelter	588	}
245	gezelter	483	}
246	mmeineke	377	}
247
248			void Thermo::velocitize() {
249
250	gezelter	608	double aVel[3], aJ[3], I[3][3];
251			int i, j, vr, vd; // velocity randomizer loop counters
252	chuckv	403	double vdrift[3];
253	mmeineke	377	double vbar;
254			const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
255			double av2;
256			double kebar;
257			int n_atoms;
258			Atom** atoms;
259			DirectionalAtom* dAtom;
260			double temperature;
261			int n_oriented;
262			int n_constraints;
263
264	mmeineke	614	atoms = info->atoms;
265			n_atoms = info->n_atoms;
266			temperature = info->target_temp;
267			n_oriented = info->n_oriented;
268			n_constraints = info->n_constraints;
269	mmeineke	377
270	tim	763	kebar = kb * temperature * (double)info->ndfRaw /
271			( 2.0 * (double)info->ndf );
272	chuckv	403
273	mmeineke	377	for(vr = 0; vr < n_atoms; vr++){
274
275			// uses equipartition theory to solve for vbar in angstrom/fs
276
277			av2 = 2.0 * kebar / atoms[vr]->getMass();
278			vbar = sqrt( av2 );
279	mmeineke	853
280	mmeineke	377	// picks random velocities from a gaussian distribution
281			// centered on vbar
282
283	gezelter	608	for (j=0; j<3; j++)
284			aVel[j] = vbar * gaussStream->getGaussian();
285
286			atoms[vr]->setVel( aVel );
287	mmeineke	377
288			}
289	chuckv	401
290			// Get the Center of Mass drift velocity.
291
292	chuckv	403	getCOMVel(vdrift);
293	mmeineke	377
294			// Corrects for the center of mass drift.
295			// sums all the momentum and divides by total mass.
296
297			for(vd = 0; vd < n_atoms; vd++){
298
299	gezelter	608	atoms[vd]->getVel(aVel);
300
301			for (j=0; j < 3; j++)
302			aVel[j] -= vdrift[j];
303	chuckv	401
304	gezelter	608	atoms[vd]->setVel( aVel );
305	mmeineke	377	}
306			if( n_oriented ){
307
308			for( i=0; i<n_atoms; i++ ){
309
310			if( atoms[i]->isDirectional() ){
311
312			dAtom = (DirectionalAtom *)atoms[i];
313	gezelter	608	dAtom->getI( I );
314
315			for (j = 0 ; j < 3; j++) {
316	mmeineke	377
317	gezelter	608	vbar = sqrt( 2.0 * kebar * I[j][j] );
318			aJ[j] = vbar * gaussStream->getGaussian();
319	mmeineke	377
320	gezelter	608	}
321
322			dAtom->setJ( aJ );
323
324	mmeineke	377	}
325			}
326			}
327			}
328	chuckv	401
329	chuckv	403	void Thermo::getCOMVel(double vdrift[3]){
330	chuckv	401
331			double mtot, mtot_local;
332	gezelter	608	double aVel[3], amass;
333	chuckv	401	double vdrift_local[3];
334	gezelter	608	int vd, n_atoms, j;
335	chuckv	401	Atom** atoms;
336
337			// We are very careless here with the distinction between n_atoms and n_local
338			// We should really fix this before someone pokes an eye out.
339
340	mmeineke	614	n_atoms = info->n_atoms;
341			atoms = info->atoms;
342	chuckv	401
343			mtot_local = 0.0;
344			vdrift_local[0] = 0.0;
345			vdrift_local[1] = 0.0;
346			vdrift_local[2] = 0.0;
347
348			for(vd = 0; vd < n_atoms; vd++){
349
350	gezelter	608	amass = atoms[vd]->getMass();
351			atoms[vd]->getVel( aVel );
352
353			for(j = 0; j < 3; j++)
354			vdrift_local[j] += aVel[j] * amass;
355	chuckv	401
356	gezelter	608	mtot_local += amass;
357	chuckv	401	}
358
359			#ifdef IS_MPI
360	mmeineke	447	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
361			MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
362	chuckv	401	#else
363			mtot = mtot_local;
364			for(vd = 0; vd < 3; vd++) {
365			vdrift[vd] = vdrift_local[vd];
366			}
367			#endif
368
369			for (vd = 0; vd < 3; vd++) {
370			vdrift[vd] = vdrift[vd] / mtot;
371			}
372
373			}
374
375	tim	763	void Thermo::getCOM(double COM[3]){
376
377			double mtot, mtot_local;
378			double aPos[3], amass;
379			double COM_local[3];
380			int i, n_atoms, j;
381			Atom** atoms;
382
383			// We are very careless here with the distinction between n_atoms and n_local
384			// We should really fix this before someone pokes an eye out.
385
386			n_atoms = info->n_atoms;
387			atoms = info->atoms;
388
389			mtot_local = 0.0;
390			COM_local[0] = 0.0;
391			COM_local[1] = 0.0;
392			COM_local[2] = 0.0;
393
394			for(i = 0; i < n_atoms; i++){
395
396			amass = atoms[i]->getMass();
397			atoms[i]->getPos( aPos );
398
399			for(j = 0; j < 3; j++)
400			COM_local[j] += aPos[j] * amass;
401
402			mtot_local += amass;
403			}
404
405			#ifdef IS_MPI
406			MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
407			MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
408			#else
409			mtot = mtot_local;
410			for(i = 0; i < 3; i++) {
411			COM[i] = COM_local[i];
412			}
413			#endif
414
415			for (i = 0; i < 3; i++) {
416			COM[i] = COM[i] / mtot;
417			}
418			}