OOPSE/libmdtools/Thermo.cpp

#include <cmath>
#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


#define BASE_SEED 123456789

Thermo::Thermo( SimInfo* the_entry_plug ) { 
  entry_plug = the_entry_plug;
  int baseSeed = BASE_SEED;
  
  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 vx2, vy2, vz2;
  double kinetic, v_sqr;
  int kl;
  double jx2, jy2, jz2; // the square of the angular momentums

  DirectionalAtom *dAtom;

  int n_atoms;
  double kinetic_global;
  Atom** atoms;

  
  n_atoms = entry_plug->n_atoms;
  atoms = entry_plug->atoms;

  kinetic = 0.0;
  kinetic_global = 0.0;
  for( kl=0; kl < n_atoms; kl++ ){

    vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
    vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
    vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();

    v_sqr = vx2 + vy2 + vz2;
    kinetic += atoms[kl]->getMass() * v_sqr;

    if( atoms[kl]->isDirectional() ){
            
      dAtom = (DirectionalAtom *)atoms[kl];
      
      jx2 = dAtom->getJx() * dAtom->getJx();    
      jy2 = dAtom->getJy() * dAtom->getJy();
      jz2 = dAtom->getJz() * dAtom->getJz();
      
      kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy()) 
        + (jz2 / dAtom->getIzz());
    }
  }
#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 = entry_plug->molecules;
  nSRI = entry_plug->n_SRI;

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

  for( el=0; el<entry_plug->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

#ifdef IS_MPI
  /*
  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
  */
#endif

  return potential;
}

double Thermo::getTotalE(){

  double total;

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

double Thermo::getTemperature(){

  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
  double temperature;
  
  temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb );
  return temperature;
}

double Thermo::getEnthalpy() {

  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
  double u, p, v;
  double press[9];

  u = this->getTotalE();

  this->getPressureTensor(press);
  p = (press[0] + press[4] + press[8]) / 3.0;

  v = this->getVolume();

  return (u + (p*v)/e_convert);
}

double Thermo::getVolume() {

  double volume;
  double Hmat[9];

  entry_plug->getBoxM(Hmat);

  // volume = h1 (dot) h2 (cross) h3

  volume = Hmat[0] * ( (Hmat[4]*Hmat[8]) - (Hmat[7]*Hmat[5]) )
         + Hmat[1] * ( (Hmat[5]*Hmat[6]) - (Hmat[8]*Hmat[3]) )
         + Hmat[2] * ( (Hmat[3]*Hmat[7]) - (Hmat[6]*Hmat[4]) );

  return volume;
}

double Thermo::getPressure() {

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

  this->getPressureTensor(press);

  pressure = p_convert * (press[0] + press[4] + press[8]) / 3.0;

  return pressure;
}


void Thermo::getPressureTensor(double press[9]){
  // 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];
  double theBox[3];
  //double* tau;
  int i, nMols;
  Molecule* molecules;

  nMols = entry_plug->n_mol;
  molecules = entry_plug->molecules;
  //tau = entry_plug->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 = entry_plug->boxVol;

  for(i=0; i<9; i++) {
    press[i] = (p_global[i] - entry_plug->tau[i]*e_convert) / volume;
  }
}

void Thermo::velocitize() {
  
  double x,y;
  double vx, vy, vz;
  double jx, jy, jz;
  int i, 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         = entry_plug->atoms;
  n_atoms       = entry_plug->n_atoms;
  temperature   = entry_plug->target_temp;
  n_oriented    = entry_plug->n_oriented;
  n_constraints = entry_plug->n_constraints;
  
  kebar = kb * temperature * (double)entry_plug->ndf / 
    ( 2.0 * (double)entry_plug->ndfRaw );
  
  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 );
 
//     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
    
    // picks random velocities from a gaussian distribution
    // centered on vbar

    vx = vbar * gaussStream->getGaussian();
    vy = vbar * gaussStream->getGaussian();
    vz = vbar * gaussStream->getGaussian();

    atoms[vr]->set_vx( vx ); 
    atoms[vr]->set_vy( vy );
    atoms[vr]->set_vz( vz );
  }

  // 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++){
    
    vx = atoms[vd]->get_vx();
    vy = atoms[vd]->get_vy();
    vz = atoms[vd]->get_vz();
        
    vx -= vdrift[0];
    vy -= vdrift[1];
    vz -= vdrift[2];
    
    atoms[vd]->set_vx(vx);
    atoms[vd]->set_vy(vy);
    atoms[vd]->set_vz(vz);
  }
  if( n_oriented ){
  
    for( i=0; i<n_atoms; i++ ){
      
      if( atoms[i]->isDirectional() ){
        
        dAtom = (DirectionalAtom *)atoms[i];

        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
        jx = vbar * gaussStream->getGaussian();

        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
        jy = vbar * gaussStream->getGaussian();
        
        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
        jz = vbar * gaussStream->getGaussian();
        
        dAtom->setJx( jx );
        dAtom->setJy( jy );
        dAtom->setJz( jz );
      }
    }   
  }
}

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

  double mtot, mtot_local;
  double vdrift_local[3];
  int vd, n_atoms;
  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 = entry_plug->n_atoms;  
  atoms   = entry_plug->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++){
    
    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
    
    mtot_local += atoms[vd]->getMass();
  }

#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;
  }
  
}

Revision:	574
Committed:	Tue Jul 8 20:56:10 2003 UTC (21 years ago) by gezelter
File size:	8632 byte(s)
Log Message:	NPTi
#	User	Rev	Content
1	mmeineke	377	#include <cmath>
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	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	402
20	mmeineke	377	#define BASE_SEED 123456789
21
22			Thermo::Thermo( SimInfo* the_entry_plug ) {
23			entry_plug = the_entry_plug;
24			int baseSeed = BASE_SEED;
25
26			gaussStream = new gaussianSPRNG( baseSeed );
27			}
28
29			Thermo::~Thermo(){
30			delete gaussStream;
31			}
32
33			double Thermo::getKinetic(){
34
35			const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
36			double vx2, vy2, vz2;
37			double kinetic, v_sqr;
38			int kl;
39			double jx2, jy2, jz2; // the square of the angular momentums
40
41			DirectionalAtom *dAtom;
42
43			int n_atoms;
44			double kinetic_global;
45			Atom** atoms;
46
47
48			n_atoms = entry_plug->n_atoms;
49			atoms = entry_plug->atoms;
50
51			kinetic = 0.0;
52			kinetic_global = 0.0;
53			for( kl=0; kl < n_atoms; kl++ ){
54
55			vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
56			vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
57			vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();
58
59			v_sqr = vx2 + vy2 + vz2;
60			kinetic += atoms[kl]->getMass() * v_sqr;
61
62			if( atoms[kl]->isDirectional() ){
63
64			dAtom = (DirectionalAtom *)atoms[kl];
65
66			jx2 = dAtom->getJx() * dAtom->getJx();
67			jy2 = dAtom->getJy() * dAtom->getJy();
68			jz2 = dAtom->getJz() * dAtom->getJz();
69
70			kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
71			+ (jz2 / dAtom->getIzz());
72			}
73			}
74			#ifdef IS_MPI
75	mmeineke	447	MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
76			MPI_SUM, MPI_COMM_WORLD);
77	mmeineke	377	kinetic = kinetic_global;
78			#endif //is_mpi
79
80			kinetic = kinetic * 0.5 / e_convert;
81
82			return kinetic;
83			}
84
85			double Thermo::getPotential(){
86
87	chuckv	401	double potential_local;
88	mmeineke	377	double potential;
89			int el, nSRI;
90	mmeineke	428	Molecule* molecules;
91	mmeineke	377
92	mmeineke	428	molecules = entry_plug->molecules;
93	mmeineke	377	nSRI = entry_plug->n_SRI;
94
95	chuckv	401	potential_local = 0.0;
96	chuckv	438	potential = 0.0;
97	chuckv	401	potential_local += entry_plug->lrPot;
98	mmeineke	377
99	mmeineke	423	for( el=0; el<entry_plug->n_mol; el++ ){
100	mmeineke	428	potential_local += molecules[el].getPotential();
101	mmeineke	377	}
102
103			// Get total potential for entire system from MPI.
104			#ifdef IS_MPI
105	mmeineke	447	MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
106			MPI_SUM, MPI_COMM_WORLD);
107	chuckv	401	#else
108			potential = potential_local;
109	mmeineke	377	#endif // is_mpi
110
111	chuckv	438	#ifdef IS_MPI
112			/*
113			std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
114			*/
115			#endif
116
117	mmeineke	377	return potential;
118			}
119
120			double Thermo::getTotalE(){
121
122			double total;
123
124			total = this->getKinetic() + this->getPotential();
125			return total;
126			}
127
128	gezelter	454	double Thermo::getTemperature(){
129
130			const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
131			double temperature;
132
133	gezelter	458	temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb );
134	mmeineke	377	return temperature;
135			}
136
137	gezelter	484	double Thermo::getEnthalpy() {
138
139			const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
140			double u, p, v;
141			double press[9];
142
143			u = this->getTotalE();
144
145			this->getPressureTensor(press);
146			p = (press[0] + press[4] + press[8]) / 3.0;
147
148			v = this->getVolume();
149
150			return (u + (p*v)/e_convert);
151			}
152
153			double Thermo::getVolume() {
154	gezelter	574
155			double volume;
156			double Hmat[9];
157
158			entry_plug->getBoxM(Hmat);
159
160			// volume = h1 (dot) h2 (cross) h3
161
162			volume = Hmat[0] * ( (Hmat[4]Hmat[8]) - (Hmat[7]Hmat[5]) )
163			+ Hmat[1] * ( (Hmat[5]Hmat[6]) - (Hmat[8]Hmat[3]) )
164			+ Hmat[2] * ( (Hmat[3]Hmat[7]) - (Hmat[6]Hmat[4]) );
165
166			return volume;
167	gezelter	484	}
168
169	gezelter	483	double Thermo::getPressure() {
170	gezelter	574
171	gezelter	483	// Relies on the calculation of the full molecular pressure tensor
172
173			const double p_convert = 1.63882576e8;
174			double press[9];
175			double pressure;
176
177			this->getPressureTensor(press);
178
179			pressure = p_convert * (press[0] + press[4] + press[8]) / 3.0;
180
181			return pressure;
182			}
183
184
185			void Thermo::getPressureTensor(double press[9]){
186			// returns pressure tensor in units amufs^-2Ang^-1
187	gezelter	445	// routine derived via viral theorem description in:
188			// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
189	mmeineke	377
190	gezelter	477	const double e_convert = 4.184e-4;
191	gezelter	483
192			double molmass, volume;
193	gezelter	468	double vcom[3];
194	gezelter	483	double p_local[9], p_global[9];
195	gezelter	468	double theBox[3];
196	mmeineke	486	//double* tau;
197	gezelter	468	int i, nMols;
198			Molecule* molecules;
199
200			nMols = entry_plug->n_mol;
201			molecules = entry_plug->molecules;
202	mmeineke	486	//tau = entry_plug->tau;
203	gezelter	468
204			// use velocities of molecular centers of mass and molecular masses:
205	gezelter	483	for (i=0; i < 9; i++) {
206			p_local[i] = 0.0;
207			p_global[i] = 0.0;
208			}
209	gezelter	475
210	gezelter	468	for (i=0; i < nMols; i++) {
211	gezelter	475	molmass = molecules[i].getCOMvel(vcom);
212	gezelter	483
213			p_local[0] += molmass * (vcom[0] * vcom[0]);
214			p_local[1] += molmass * (vcom[0] * vcom[1]);
215			p_local[2] += molmass * (vcom[0] * vcom[2]);
216			p_local[3] += molmass * (vcom[1] * vcom[0]);
217			p_local[4] += molmass * (vcom[1] * vcom[1]);
218			p_local[5] += molmass * (vcom[1] * vcom[2]);
219			p_local[6] += molmass * (vcom[2] * vcom[0]);
220			p_local[7] += molmass * (vcom[2] * vcom[1]);
221			p_local[8] += molmass * (vcom[2] * vcom[2]);
222	gezelter	468	}
223
224			// Get total for entire system from MPI.
225	chuckv	479
226	gezelter	468	#ifdef IS_MPI
227	gezelter	483	MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
228	gezelter	468	#else
229	gezelter	483	for (i=0; i<9; i++) {
230			p_global[i] = p_local[i];
231			}
232	gezelter	468	#endif // is_mpi
233
234	mmeineke	572	volume = entry_plug->boxVol;
235	gezelter	468
236	gezelter	483	for(i=0; i<9; i++) {
237	mmeineke	486	press[i] = (p_global[i] - entry_plug->tau[i]*e_convert) / volume;
238	gezelter	483	}
239	mmeineke	377	}
240
241			void Thermo::velocitize() {
242
243			double x,y;
244			double vx, vy, vz;
245			double jx, jy, jz;
246			int i, vr, vd; // velocity randomizer loop counters
247	chuckv	403	double vdrift[3];
248	mmeineke	377	double vbar;
249			const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
250			double av2;
251			double kebar;
252			int n_atoms;
253			Atom** atoms;
254			DirectionalAtom* dAtom;
255			double temperature;
256			int n_oriented;
257			int n_constraints;
258
259			atoms = entry_plug->atoms;
260			n_atoms = entry_plug->n_atoms;
261			temperature = entry_plug->target_temp;
262			n_oriented = entry_plug->n_oriented;
263			n_constraints = entry_plug->n_constraints;
264
265	gezelter	458	kebar = kb * temperature * (double)entry_plug->ndf /
266			( 2.0 * (double)entry_plug->ndfRaw );
267	chuckv	403
268	mmeineke	377	for(vr = 0; vr < n_atoms; vr++){
269
270			// uses equipartition theory to solve for vbar in angstrom/fs
271
272			av2 = 2.0 * kebar / atoms[vr]->getMass();
273			vbar = sqrt( av2 );
274	gezelter	444
275	mmeineke	377	// vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
276
277			// picks random velocities from a gaussian distribution
278			// centered on vbar
279
280			vx = vbar * gaussStream->getGaussian();
281			vy = vbar * gaussStream->getGaussian();
282			vz = vbar * gaussStream->getGaussian();
283
284			atoms[vr]->set_vx( vx );
285			atoms[vr]->set_vy( vy );
286			atoms[vr]->set_vz( vz );
287			}
288	chuckv	401
289			// Get the Center of Mass drift velocity.
290
291	chuckv	403	getCOMVel(vdrift);
292	mmeineke	377
293			// Corrects for the center of mass drift.
294			// sums all the momentum and divides by total mass.
295
296			for(vd = 0; vd < n_atoms; vd++){
297
298			vx = atoms[vd]->get_vx();
299			vy = atoms[vd]->get_vy();
300			vz = atoms[vd]->get_vz();
301	chuckv	401
302	mmeineke	377	vx -= vdrift[0];
303			vy -= vdrift[1];
304			vz -= vdrift[2];
305
306			atoms[vd]->set_vx(vx);
307			atoms[vd]->set_vy(vy);
308			atoms[vd]->set_vz(vz);
309			}
310			if( n_oriented ){
311
312			for( i=0; i<n_atoms; i++ ){
313
314			if( atoms[i]->isDirectional() ){
315
316			dAtom = (DirectionalAtom *)atoms[i];
317
318			vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
319			jx = vbar * gaussStream->getGaussian();
320
321			vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
322			jy = vbar * gaussStream->getGaussian();
323	gezelter	454
324	mmeineke	377	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
325			jz = vbar * gaussStream->getGaussian();
326
327			dAtom->setJx( jx );
328			dAtom->setJy( jy );
329			dAtom->setJz( jz );
330			}
331			}
332			}
333			}
334	chuckv	401
335	chuckv	403	void Thermo::getCOMVel(double vdrift[3]){
336	chuckv	401
337			double mtot, mtot_local;
338			double vdrift_local[3];
339			int vd, n_atoms;
340			Atom** atoms;
341
342			// We are very careless here with the distinction between n_atoms and n_local
343			// We should really fix this before someone pokes an eye out.
344
345			n_atoms = entry_plug->n_atoms;
346			atoms = entry_plug->atoms;
347
348			mtot_local = 0.0;
349			vdrift_local[0] = 0.0;
350			vdrift_local[1] = 0.0;
351			vdrift_local[2] = 0.0;
352
353			for(vd = 0; vd < n_atoms; vd++){
354
355			vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
356			vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
357			vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
358
359			mtot_local += atoms[vd]->getMass();
360			}
361
362			#ifdef IS_MPI
363	mmeineke	447	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
364			MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
365	chuckv	401	#else
366			mtot = mtot_local;
367			for(vd = 0; vd < 3; vd++) {
368			vdrift[vd] = vdrift_local[vd];
369			}
370			#endif
371
372			for (vd = 0; vd < 3; vd++) {
373			vdrift[vd] = vdrift[vd] / mtot;
374			}
375
376			}
377