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"
#include "MatVec3.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], pcom[3], fcom[3], scaled[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 < info->integrableObjects.size(); i++) {

    molmass = info->integrableObjects[i]->getMass();
    
    info->integrableObjects[i]->getVel(vcom);
    info->integrableObjects[i]->getPos(pcom);
    info->integrableObjects[i]->getFrc(fcom);

    matVecMul3(info->HmatInv, pcom, scaled);
  
    for(j=0; j<3; j++)
      scaled[j] -= roundMe(scaled[j]);

    // calc the wrapped real coordinates from the wrapped scaled coordinates
  
    matVecMul3(info->Hmat, scaled, pcom);
    
    p_local[0] += molmass * (vcom[0] * vcom[0]) + fcom[0]*pcom[0]*eConvert; 
    p_local[1] += molmass * (vcom[0] * vcom[1]) + fcom[0]*pcom[1]*eConvert; 
    p_local[2] += molmass * (vcom[0] * vcom[2]) + fcom[0]*pcom[2]*eConvert; 
    p_local[3] += molmass * (vcom[1] * vcom[0]) + fcom[1]*pcom[0]*eConvert; 
    p_local[4] += molmass * (vcom[1] * vcom[1]) + fcom[1]*pcom[1]*eConvert; 
    p_local[5] += molmass * (vcom[1] * vcom[2]) + fcom[1]*pcom[2]*eConvert; 
    p_local[6] += molmass * (vcom[2] * vcom[0]) + fcom[2]*pcom[0]*eConvert; 
    p_local[7] += molmass * (vcom[2] * vcom[1]) + fcom[2]*pcom[1]*eConvert; 
    p_local[8] += molmass * (vcom[2] * vcom[2]) + fcom[2]*pcom[2]*eConvert; 
    
  }

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