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

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], 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:	1133
Committed:	Mon Apr 26 14:29:18 2004 UTC (20 years, 2 months ago) by gezelter
File size:	10349 byte(s)
Log Message:	Fixed a bug in calc_charge_charge.F90
#	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	gezelter	1133	inline double roundMe( double x ){
21			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
22			}
23
24	mmeineke	614	Thermo::Thermo( SimInfo* the_info ) {
25			info = the_info;
26	tim	708	int baseSeed = the_info->getSeed();
27	mmeineke	377
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	gezelter	608	double kinetic;
39			double amass;
40			double aVel[3], aJ[3], I[3][3];
41	tim	1118	int i, j, k, kl;
42	mmeineke	377
43			double kinetic_global;
44	tim	1113	vector<StuntDouble *> integrableObjects = info->integrableObjects;
45	mmeineke	377
46			kinetic = 0.0;
47			kinetic_global = 0.0;
48
49	tim	1113	for (kl=0; kl<integrableObjects.size(); kl++) {
50			integrableObjects[kl]->getVel(aVel);
51			amass = integrableObjects[kl]->getMass();
52	gezelter	608
53	tim	1113	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	tim	1118	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	gezelter	1125	for (j=0; j<3; j++)
68			kinetic += aJ[j]*aJ[j] / I[j][j];
69	tim	1118	}
70	gezelter	1125	}
71	mmeineke	377	}
72			#ifdef IS_MPI
73	mmeineke	447	MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
74			MPI_SUM, MPI_COMM_WORLD);
75	mmeineke	377	kinetic = kinetic_global;
76			#endif //is_mpi
77	gezelter	1125
78	mmeineke	377	kinetic = kinetic * 0.5 / e_convert;
79
80			return kinetic;
81			}
82
83			double Thermo::getPotential(){
84
85	chuckv	401	double potential_local;
86	mmeineke	377	double potential;
87			int el, nSRI;
88	mmeineke	428	Molecule* molecules;
89	mmeineke	377
90	mmeineke	614	molecules = info->molecules;
91			nSRI = info->n_SRI;
92	mmeineke	377
93	chuckv	401	potential_local = 0.0;
94	chuckv	438	potential = 0.0;
95	mmeineke	614	potential_local += info->lrPot;
96	mmeineke	377
97	mmeineke	614	for( el=0; el<info->n_mol; el++ ){
98	mmeineke	428	potential_local += molecules[el].getPotential();
99	mmeineke	377	}
100
101			// Get total potential for entire system from MPI.
102			#ifdef IS_MPI
103	mmeineke	447	MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
104			MPI_SUM, MPI_COMM_WORLD);
105	chuckv	401	#else
106			potential = potential_local;
107	mmeineke	377	#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	gezelter	454	double Thermo::getTemperature(){
121
122	tim	763	const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
123	gezelter	454	double temperature;
124	tim	1113
125	mmeineke	614	temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
126	mmeineke	377	return temperature;
127			}
128
129	gezelter	484	double Thermo::getVolume() {
130	gezelter	574
131	mmeineke	614	return info->boxVol;
132	gezelter	484	}
133
134	gezelter	483	double Thermo::getPressure() {
135	gezelter	574
136	gezelter	483	// Relies on the calculation of the full molecular pressure tensor
137
138			const double p_convert = 1.63882576e8;
139	gezelter	588	double press[3][3];
140	gezelter	483	double pressure;
141
142			this->getPressureTensor(press);
143
144	gezelter	588	pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
145	gezelter	483
146			return pressure;
147			}
148
149	mmeineke	755	double Thermo::getPressureX() {
150	gezelter	483
151	mmeineke	755	// 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	gezelter	588	void Thermo::getPressureTensor(double press[3][3]){
196	gezelter	483	// returns pressure tensor in units amufs^-2Ang^-1
197	gezelter	445	// routine derived via viral theorem description in:
198			// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
199	mmeineke	377
200	gezelter	477	const double e_convert = 4.184e-4;
201	gezelter	483
202			double molmass, volume;
203	tim	1131	double vcom[3], pcom[3], fcom[3], scaled[3];
204	gezelter	483	double p_local[9], p_global[9];
205	mmeineke	614	int i, j, k, nMols;
206	gezelter	468	Molecule* molecules;
207
208	mmeineke	614	nMols = info->n_mol;
209			molecules = info->molecules;
210			//tau = info->tau;
211	gezelter	468
212			// use velocities of molecular centers of mass and molecular masses:
213	gezelter	483	for (i=0; i < 9; i++) {
214			p_local[i] = 0.0;
215			p_global[i] = 0.0;
216			}
217	gezelter	475
218	tim	1131	for (i=0; i < info->integrableObjects.size(); i++) {
219	gezelter	483
220	tim	1131	molmass = info->integrableObjects[i]->getMass();
221
222			info->integrableObjects[i]->getVel(vcom);
223			info->integrableObjects[i]->getPos(pcom);
224			info->integrableObjects[i]->getFrc(fcom);
225
226			matVecMul3(info->HmatInv, pcom, scaled);
227
228			for(j=0; j<3; j++)
229			scaled[j] -= roundMe(scaled[j]);
230
231			// calc the wrapped real coordinates from the wrapped scaled coordinates
232
233			matVecMul3(info->Hmat, scaled, pcom);
234
235			p_local[0] += molmass * (vcom[0] * vcom[0]) + fcom[0]pcom[0]eConvert;
236			p_local[1] += molmass * (vcom[0] * vcom[1]) + fcom[0]pcom[1]eConvert;
237			p_local[2] += molmass * (vcom[0] * vcom[2]) + fcom[0]pcom[2]eConvert;
238			p_local[3] += molmass * (vcom[1] * vcom[0]) + fcom[1]pcom[0]eConvert;
239			p_local[4] += molmass * (vcom[1] * vcom[1]) + fcom[1]pcom[1]eConvert;
240			p_local[5] += molmass * (vcom[1] * vcom[2]) + fcom[1]pcom[2]eConvert;
241			p_local[6] += molmass * (vcom[2] * vcom[0]) + fcom[2]pcom[0]eConvert;
242			p_local[7] += molmass * (vcom[2] * vcom[1]) + fcom[2]pcom[1]eConvert;
243			p_local[8] += molmass * (vcom[2] * vcom[2]) + fcom[2]pcom[2]eConvert;
244
245	gezelter	468	}
246
247			// Get total for entire system from MPI.
248	chuckv	479
249	gezelter	468	#ifdef IS_MPI
250	gezelter	483	MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
251	gezelter	468	#else
252	gezelter	483	for (i=0; i<9; i++) {
253			p_global[i] = p_local[i];
254			}
255	gezelter	468	#endif // is_mpi
256
257	gezelter	611	volume = this->getVolume();
258	gezelter	468
259	gezelter	588	for(i = 0; i < 3; i++) {
260			for (j = 0; j < 3; j++) {
261			k = 3*i + j;
262	tim	1131	press[i][j] = p_global[k] / volume;
263	mmeineke	614
264	gezelter	588	}
265	gezelter	483	}
266	mmeineke	377	}
267
268			void Thermo::velocitize() {
269
270	gezelter	608	double aVel[3], aJ[3], I[3][3];
271	tim	1127	int i, j, l, m, n, vr, vd; // velocity randomizer loop counters
272	chuckv	403	double vdrift[3];
273	mmeineke	377	double vbar;
274			const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
275			double av2;
276			double kebar;
277			double temperature;
278	tim	1127	int nobj;
279	mmeineke	377
280	tim	1127	nobj = info->integrableObjects.size();
281
282	mmeineke	614	temperature = info->target_temp;
283	mmeineke	377
284	tim	763	kebar = kb * temperature * (double)info->ndfRaw /
285			( 2.0 * (double)info->ndf );
286	chuckv	403
287	tim	1127	for(vr = 0; vr < nobj; vr++){
288	mmeineke	377
289			// uses equipartition theory to solve for vbar in angstrom/fs
290
291	tim	1127	av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
292	mmeineke	377	vbar = sqrt( av2 );
293	mmeineke	853
294	mmeineke	377	// picks random velocities from a gaussian distribution
295			// centered on vbar
296
297	gezelter	608	for (j=0; j<3; j++)
298			aVel[j] = vbar * gaussStream->getGaussian();
299
300	tim	1127	info->integrableObjects[vr]->setVel( aVel );
301
302			if(info->integrableObjects[vr]->isDirectional()){
303	mmeineke	377
304	tim	1127	info->integrableObjects[vr]->getI( I );
305
306			if (info->integrableObjects[vr]->isLinear()) {
307
308			l= info->integrableObjects[vr]->linearAxis();
309			m = (l+1)%3;
310			n = (l+2)%3;
311
312			aJ[l] = 0.0;
313			vbar = sqrt( 2.0 * kebar * I[m][m] );
314			aJ[m] = vbar * gaussStream->getGaussian();
315			vbar = sqrt( 2.0 * kebar * I[n][n] );
316			aJ[n] = vbar * gaussStream->getGaussian();
317
318			} else {
319			for (j = 0 ; j < 3; j++) {
320			vbar = sqrt( 2.0 * kebar * I[j][j] );
321			aJ[j] = vbar * gaussStream->getGaussian();
322			}
323			} // else isLinear
324
325			info->integrableObjects[vr]->setJ( aJ );
326
327			}//isDirectional
328
329	mmeineke	377	}
330	chuckv	401
331			// Get the Center of Mass drift velocity.
332
333	chuckv	403	getCOMVel(vdrift);
334	mmeineke	377
335			// Corrects for the center of mass drift.
336			// sums all the momentum and divides by total mass.
337
338	tim	1127	for(vd = 0; vd < nobj; vd++){
339	mmeineke	377
340	tim	1127	info->integrableObjects[vd]->getVel(aVel);
341	gezelter	608
342			for (j=0; j < 3; j++)
343			aVel[j] -= vdrift[j];
344	chuckv	401
345	tim	1127	info->integrableObjects[vd]->setVel( aVel );
346	mmeineke	377	}
347
348			}
349	chuckv	401
350	chuckv	403	void Thermo::getCOMVel(double vdrift[3]){
351	chuckv	401
352			double mtot, mtot_local;
353	gezelter	608	double aVel[3], amass;
354	chuckv	401	double vdrift_local[3];
355	tim	1127	int vd, j;
356			int nobj;
357	chuckv	401
358	tim	1127	nobj = info->integrableObjects.size();
359	chuckv	401
360			mtot_local = 0.0;
361			vdrift_local[0] = 0.0;
362			vdrift_local[1] = 0.0;
363			vdrift_local[2] = 0.0;
364
365	tim	1127	for(vd = 0; vd < nobj; vd++){
366	chuckv	401
367	tim	1127	amass = info->integrableObjects[vd]->getMass();
368			info->integrableObjects[vd]->getVel( aVel );
369	gezelter	608
370			for(j = 0; j < 3; j++)
371			vdrift_local[j] += aVel[j] * amass;
372	chuckv	401
373	gezelter	608	mtot_local += amass;
374	chuckv	401	}
375
376			#ifdef IS_MPI
377	mmeineke	447	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
378			MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
379	chuckv	401	#else
380			mtot = mtot_local;
381			for(vd = 0; vd < 3; vd++) {
382			vdrift[vd] = vdrift_local[vd];
383			}
384			#endif
385
386			for (vd = 0; vd < 3; vd++) {
387			vdrift[vd] = vdrift[vd] / mtot;
388			}
389
390			}
391
392	tim	763	void Thermo::getCOM(double COM[3]){
393
394			double mtot, mtot_local;
395			double aPos[3], amass;
396			double COM_local[3];
397	tim	1127	int i, j;
398			int nobj;
399	tim	763
400			mtot_local = 0.0;
401			COM_local[0] = 0.0;
402			COM_local[1] = 0.0;
403			COM_local[2] = 0.0;
404	tim	1127
405			nobj = info->integrableObjects.size();
406			for(i = 0; i < nobj; i++){
407	tim	763
408	tim	1127	amass = info->integrableObjects[i]->getMass();
409			info->integrableObjects[i]->getPos( aPos );
410	tim	763
411			for(j = 0; j < 3; j++)
412			COM_local[j] += aPos[j] * amass;
413
414			mtot_local += amass;
415			}
416
417			#ifdef IS_MPI
418			MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
419			MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
420			#else
421			mtot = mtot_local;
422			for(i = 0; i < 3; i++) {
423			COM[i] = COM_local[i];
424			}
425			#endif
426
427			for (i = 0; i < 3; i++) {
428			COM[i] = COM[i] / mtot;
429			}
430			}
431	tim	1127
432			void Thermo::removeCOMdrift() {
433			double vdrift[3], aVel[3];
434			int vd, j, nobj;
435
436			nobj = info->integrableObjects.size();
437
438			// Get the Center of Mass drift velocity.
439
440			getCOMVel(vdrift);
441
442			// Corrects for the center of mass drift.
443			// sums all the momentum and divides by total mass.
444
445			for(vd = 0; vd < nobj; vd++){
446
447			info->integrableObjects[vd]->getVel(aVel);
448
449			for (j=0; j < 3; j++)
450			aVel[j] -= vdrift[j];
451
452			info->integrableObjects[vd]->setVel( aVel );
453			}
454	gezelter	1133	}