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"
#include "ConstraintManager.hpp"
#include "Mat3x3d.hpp"

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

  cpIter = info->consMan->createPairIterator();
}

Thermo::~Thermo(){
  delete gaussStream;
  delete cpIter;
}

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;

  for (i=0; i < 9; i++) {    
    p_local[i] = 0.0;
    p_global[i] = 0.0;
  }

  // use velocities of integrableObjects and their masses:  

  for (i=0; i < info->integrableObjects.size(); i++) {

    molmass = info->integrableObjects[i]->getMass();
    
    info->integrableObjects[i]->getVel(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 );
  }
}

void Thermo::removeAngularMomentum(){
  Vector3d vcom;
  Vector3d qcom;
  Vector3d pos;
  Vector3d vel;
  double mass;  
  double xx;
  double yy;
  double zz;
  double xy;
  double xz;
  double yz;
  Vector3d localAngMom;
  Vector3d angMom;
  Vector3d omega;
  vector<StuntDouble *> integrableObjects;
  double localInertiaVec[9];
  double inertiaVec[9];
  vector<Vector3d> qMinusQCom;
  vector<Vector3d> vMinusVCom;
  Mat3x3d inertiaMat;
  Mat3x3d inverseInertiaMat;
  
  integrableObjects = info->integrableObjects;
  qMinusQCom.resize(integrableObjects.size());
  vMinusVCom.resize(integrableObjects.size());
  
  getCOM(qcom.vec);
  getCOMVel(vcom.vec);
        
  //initialize components for inertia tensor
  xx = 0.0;
  yy = 0.0;
  zz = 0.0;
  xy = 0.0;
  xz = 0.0;
  yz = 0.0;
  
   //build components of Inertia tensor
  //
  //       [  Ixx -Ixy  -Ixz ]
  //   J = | -Iyx  Iyy  -Iyz |
  //       [ -Izx -Iyz   Izz ]
  //See Fowles and Cassidy Chapter 9 or Goldstein Chapter 5
  for(size_t i = 0; i < integrableObjects.size(); i++){
    integrableObjects[i]->getPos(pos.vec);
    integrableObjects[i]->getVel(vel.vec);
    mass = integrableObjects[i]->getMass();
    
    qMinusQCom[i] = pos - qcom; 
    info->wrapVector(qMinusQCom[i].vec);
    
    vMinusVCom[i] = vel - vcom;

    //compute moment of inertia coefficents
    xx += qMinusQCom[i].x * qMinusQCom[i].x * mass;
    yy += qMinusQCom[i].y * qMinusQCom[i].y * mass;
    zz += qMinusQCom[i].z * qMinusQCom[i].z * mass;

    // compute products of inertia
    xy += qMinusQCom[i].x * qMinusQCom[i].y * mass;
    xz += qMinusQCom[i].x * qMinusQCom[i].z * mass;
    yz += qMinusQCom[i].y * qMinusQCom[i].z * mass;

    localAngMom += crossProduct(qMinusQCom[i] , vMinusVCom[i] ) * mass;
    
  }

  localInertiaVec[0] =yy+zz;
  localInertiaVec[1] = -xy;
  localInertiaVec[2] = -xz;
  localInertiaVec[3] = -xy;
  localInertiaVec[4] = xx+zz;
  localInertiaVec[5] = -yz;
  localInertiaVec[6] = -xz;
  localInertiaVec[7] = -yz;
  localInertiaVec[8] = xx+yy;

  //Sum and distribute inertia and angmom arrays
#ifdef MPI

  MPI_Allreduce(localInertiaVec, inertiaVec, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  MPI_Allreduce(localAngMom.vec, angMom.vec, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  inertiaMat.element[0][0] = inertiaVec[0];
  inertiaMat.element[0][1] = inertiaVec[1];
  inertiaMat.element[0][2] = inertiaVec[2];

  inertiaMat.element[1][0] = inertiaVec[3];
  inertiaMat.element[1][1] = inertiaVec[4];
  inertiaMat.element[1][2] = inertiaVec[5];

  inertiaMat.element[2][0] = inertiaVec[6];
  inertiaMat.element[2][1] = inertiaVec[7];
  inertiaMat.element[2][2] = inertiaVec[8];

#else

    inertiaMat.element[0][0] = localInertiaVec[0];
    inertiaMat.element[0][1] = localInertiaVec[1];
    inertiaMat.element[0][2] = localInertiaVec[2];

    inertiaMat.element[1][0] = localInertiaVec[3];
    inertiaMat.element[1][1] = localInertiaVec[4];
    inertiaMat.element[1][2] = localInertiaVec[5];

    inertiaMat.element[2][0] = localInertiaVec[6];
    inertiaMat.element[2][1] = localInertiaVec[7];
    inertiaMat.element[2][2] = localInertiaVec[8];
   
    angMom     = localAngMom;
#endif

    //invert the moment of inertia tensor by LU-decomposition / backsolving:

    inverseInertiaMat = inertiaMat.inverse();

    //calculate the angular velocities: omega = I^-1 . L

    omega = inverseInertiaMat * angMom;

    //subtract out center of mass velocity and angular momentum from
    //particle velocities

    for(size_t i = 0; i < integrableObjects.size(); i++){
      vel = vMinusVCom[i] - crossProduct(omega, qMinusQCom[i]);
      integrableObjects[i]->setVel(vel.vec);      
    }
}

double Thermo::getConsEnergy(){
  ConstraintPair* consPair;
  double totConsEnergy;
  double bondLen2;
  double dist;
  double lamda;
  
  totConsEnergy = 0;
  
  for(cpIter->first(); !cpIter->isEnd(); cpIter->next()){
    consPair =  cpIter->currentItem();
    bondLen2 = consPair->getBondLength2();
    lamda = consPair->getLamda();
    //dist = consPair->getDistance();

    //totConsEnergy += lamda * (dist*dist - bondLen2);
  }

  return totConsEnergy;
}


Revision:	1452
Committed:	Mon Aug 23 15:11:36 2004 UTC (19 years, 10 months ago) by tim
File size:	13914 byte(s)
Log Message:	* empty log message *
#	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	tim	1452	#include "ConstraintManager.hpp"
15			#include "Mat3x3d.hpp"
16	mmeineke	402
17			#ifdef IS_MPI
18	chuckv	401	#define __C
19	mmeineke	402	#include "mpiSimulation.hpp"
20			#endif // is_mpi
21	mmeineke	377
22	gezelter	1133	inline double roundMe( double x ){
23			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
24			}
25
26	mmeineke	614	Thermo::Thermo( SimInfo* the_info ) {
27			info = the_info;
28	tim	708	int baseSeed = the_info->getSeed();
29	mmeineke	377
30			gaussStream = new gaussianSPRNG( baseSeed );
31	tim	1452
32			cpIter = info->consMan->createPairIterator();
33	mmeineke	377	}
34
35			Thermo::~Thermo(){
36			delete gaussStream;
37	tim	1452	delete cpIter;
38	mmeineke	377	}
39
40			double Thermo::getKinetic(){
41
42			const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
43	gezelter	608	double kinetic;
44			double amass;
45			double aVel[3], aJ[3], I[3][3];
46	tim	1118	int i, j, k, kl;
47	mmeineke	377
48			double kinetic_global;
49	tim	1113	vector<StuntDouble *> integrableObjects = info->integrableObjects;
50	mmeineke	377
51			kinetic = 0.0;
52			kinetic_global = 0.0;
53
54	tim	1113	for (kl=0; kl<integrableObjects.size(); kl++) {
55			integrableObjects[kl]->getVel(aVel);
56			amass = integrableObjects[kl]->getMass();
57	gezelter	608
58	tim	1113	for(j=0; j<3; j++)
59			kinetic += amassaVel[j]aVel[j];
60
61			if (integrableObjects[kl]->isDirectional()){
62
63			integrableObjects[kl]->getJ( aJ );
64			integrableObjects[kl]->getI( I );
65
66	tim	1118	if (integrableObjects[kl]->isLinear()) {
67			i = integrableObjects[kl]->linearAxis();
68			j = (i+1)%3;
69			k = (i+2)%3;
70			kinetic += aJ[j]aJ[j]/I[j][j] + aJ[k]aJ[k]/I[k][k];
71			} else {
72	gezelter	1125	for (j=0; j<3; j++)
73			kinetic += aJ[j]*aJ[j] / I[j][j];
74	tim	1118	}
75	gezelter	1125	}
76	mmeineke	377	}
77			#ifdef IS_MPI
78	mmeineke	447	MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
79			MPI_SUM, MPI_COMM_WORLD);
80	mmeineke	377	kinetic = kinetic_global;
81			#endif //is_mpi
82	gezelter	1125
83	mmeineke	377	kinetic = kinetic * 0.5 / e_convert;
84
85			return kinetic;
86			}
87
88			double Thermo::getPotential(){
89
90	chuckv	401	double potential_local;
91	mmeineke	377	double potential;
92			int el, nSRI;
93	mmeineke	428	Molecule* molecules;
94	mmeineke	377
95	mmeineke	614	molecules = info->molecules;
96			nSRI = info->n_SRI;
97	mmeineke	377
98	chuckv	401	potential_local = 0.0;
99	chuckv	438	potential = 0.0;
100	mmeineke	614	potential_local += info->lrPot;
101	mmeineke	377
102	mmeineke	614	for( el=0; el<info->n_mol; el++ ){
103	mmeineke	428	potential_local += molecules[el].getPotential();
104	mmeineke	377	}
105
106			// Get total potential for entire system from MPI.
107			#ifdef IS_MPI
108	mmeineke	447	MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
109			MPI_SUM, MPI_COMM_WORLD);
110	chuckv	401	#else
111			potential = potential_local;
112	mmeineke	377	#endif // is_mpi
113
114			return potential;
115			}
116
117			double Thermo::getTotalE(){
118
119			double total;
120
121			total = this->getKinetic() + this->getPotential();
122			return total;
123			}
124
125	gezelter	454	double Thermo::getTemperature(){
126
127	tim	763	const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
128	gezelter	454	double temperature;
129	tim	1113
130	mmeineke	614	temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
131	mmeineke	377	return temperature;
132			}
133
134	gezelter	484	double Thermo::getVolume() {
135	gezelter	574
136	mmeineke	614	return info->boxVol;
137	gezelter	484	}
138
139	gezelter	483	double Thermo::getPressure() {
140	gezelter	574
141	gezelter	483	// Relies on the calculation of the full molecular pressure tensor
142
143			const double p_convert = 1.63882576e8;
144	gezelter	588	double press[3][3];
145	gezelter	483	double pressure;
146
147			this->getPressureTensor(press);
148
149	gezelter	588	pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
150	gezelter	483
151			return pressure;
152			}
153
154	mmeineke	755	double Thermo::getPressureX() {
155	gezelter	483
156	mmeineke	755	// Relies on the calculation of the full molecular pressure tensor
157
158			const double p_convert = 1.63882576e8;
159			double press[3][3];
160			double pressureX;
161
162			this->getPressureTensor(press);
163
164			pressureX = p_convert * press[0][0];
165
166			return pressureX;
167			}
168
169			double Thermo::getPressureY() {
170
171			// Relies on the calculation of the full molecular pressure tensor
172
173			const double p_convert = 1.63882576e8;
174			double press[3][3];
175			double pressureY;
176
177			this->getPressureTensor(press);
178
179			pressureY = p_convert * press[1][1];
180
181			return pressureY;
182			}
183
184			double Thermo::getPressureZ() {
185
186			// Relies on the calculation of the full molecular pressure tensor
187
188			const double p_convert = 1.63882576e8;
189			double press[3][3];
190			double pressureZ;
191
192			this->getPressureTensor(press);
193
194			pressureZ = p_convert * press[2][2];
195
196			return pressureZ;
197			}
198
199
200	gezelter	588	void Thermo::getPressureTensor(double press[3][3]){
201	gezelter	483	// returns pressure tensor in units amufs^-2Ang^-1
202	gezelter	445	// routine derived via viral theorem description in:
203			// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
204	mmeineke	377
205	gezelter	477	const double e_convert = 4.184e-4;
206	gezelter	483
207			double molmass, volume;
208	chrisfen	1251	double vcom[3];
209	gezelter	483	double p_local[9], p_global[9];
210	chrisfen	1251	int i, j, k;
211	gezelter	468
212	gezelter	483	for (i=0; i < 9; i++) {
213			p_local[i] = 0.0;
214			p_global[i] = 0.0;
215			}
216	gezelter	475
217	chrisfen	1251	// use velocities of integrableObjects and their masses:
218
219	tim	1131	for (i=0; i < info->integrableObjects.size(); i++) {
220	gezelter	483
221	tim	1131	molmass = info->integrableObjects[i]->getMass();
222
223			info->integrableObjects[i]->getVel(vcom);
224
225	gezelter	1192	p_local[0] += molmass * (vcom[0] * vcom[0]);
226			p_local[1] += molmass * (vcom[0] * vcom[1]);
227			p_local[2] += molmass * (vcom[0] * vcom[2]);
228			p_local[3] += molmass * (vcom[1] * vcom[0]);
229			p_local[4] += molmass * (vcom[1] * vcom[1]);
230			p_local[5] += molmass * (vcom[1] * vcom[2]);
231			p_local[6] += molmass * (vcom[2] * vcom[0]);
232			p_local[7] += molmass * (vcom[2] * vcom[1]);
233			p_local[8] += molmass * (vcom[2] * vcom[2]);
234	gezelter	1253
235	gezelter	468	}
236
237			// Get total for entire system from MPI.
238	chrisfen	1251
239	gezelter	468	#ifdef IS_MPI
240	gezelter	483	MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
241	gezelter	468	#else
242	gezelter	483	for (i=0; i<9; i++) {
243			p_global[i] = p_local[i];
244			}
245	gezelter	468	#endif // is_mpi
246
247	gezelter	611	volume = this->getVolume();
248	gezelter	468
249	gezelter	1253
250
251	gezelter	588	for(i = 0; i < 3; i++) {
252			for (j = 0; j < 3; j++) {
253			k = 3*i + j;
254	gezelter	1192	press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
255	gezelter	588	}
256	gezelter	483	}
257	mmeineke	377	}
258
259			void Thermo::velocitize() {
260
261	gezelter	608	double aVel[3], aJ[3], I[3][3];
262	tim	1127	int i, j, l, m, n, vr, vd; // velocity randomizer loop counters
263	chuckv	403	double vdrift[3];
264	mmeineke	377	double vbar;
265			const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
266			double av2;
267			double kebar;
268			double temperature;
269	tim	1127	int nobj;
270	mmeineke	377
271	tim	1127	nobj = info->integrableObjects.size();
272
273	mmeineke	614	temperature = info->target_temp;
274	mmeineke	377
275	tim	763	kebar = kb * temperature * (double)info->ndfRaw /
276			( 2.0 * (double)info->ndf );
277	chuckv	403
278	tim	1127	for(vr = 0; vr < nobj; vr++){
279	mmeineke	377
280			// uses equipartition theory to solve for vbar in angstrom/fs
281
282	tim	1127	av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
283	mmeineke	377	vbar = sqrt( av2 );
284	mmeineke	853
285	mmeineke	377	// picks random velocities from a gaussian distribution
286			// centered on vbar
287
288	gezelter	608	for (j=0; j<3; j++)
289			aVel[j] = vbar * gaussStream->getGaussian();
290
291	tim	1127	info->integrableObjects[vr]->setVel( aVel );
292
293			if(info->integrableObjects[vr]->isDirectional()){
294	mmeineke	377
295	tim	1127	info->integrableObjects[vr]->getI( I );
296
297			if (info->integrableObjects[vr]->isLinear()) {
298
299			l= info->integrableObjects[vr]->linearAxis();
300			m = (l+1)%3;
301			n = (l+2)%3;
302
303			aJ[l] = 0.0;
304			vbar = sqrt( 2.0 * kebar * I[m][m] );
305			aJ[m] = vbar * gaussStream->getGaussian();
306			vbar = sqrt( 2.0 * kebar * I[n][n] );
307			aJ[n] = vbar * gaussStream->getGaussian();
308
309			} else {
310			for (j = 0 ; j < 3; j++) {
311			vbar = sqrt( 2.0 * kebar * I[j][j] );
312			aJ[j] = vbar * gaussStream->getGaussian();
313			}
314			} // else isLinear
315
316			info->integrableObjects[vr]->setJ( aJ );
317
318			}//isDirectional
319
320	mmeineke	377	}
321	chuckv	401
322			// Get the Center of Mass drift velocity.
323
324	chuckv	403	getCOMVel(vdrift);
325	mmeineke	377
326			// Corrects for the center of mass drift.
327			// sums all the momentum and divides by total mass.
328
329	tim	1127	for(vd = 0; vd < nobj; vd++){
330	mmeineke	377
331	tim	1127	info->integrableObjects[vd]->getVel(aVel);
332	gezelter	608
333			for (j=0; j < 3; j++)
334			aVel[j] -= vdrift[j];
335	chuckv	401
336	tim	1127	info->integrableObjects[vd]->setVel( aVel );
337	mmeineke	377	}
338
339			}
340	chuckv	401
341	chuckv	403	void Thermo::getCOMVel(double vdrift[3]){
342	chuckv	401
343			double mtot, mtot_local;
344	gezelter	608	double aVel[3], amass;
345	chuckv	401	double vdrift_local[3];
346	tim	1127	int vd, j;
347			int nobj;
348	chuckv	401
349	tim	1127	nobj = info->integrableObjects.size();
350	chuckv	401
351			mtot_local = 0.0;
352			vdrift_local[0] = 0.0;
353			vdrift_local[1] = 0.0;
354			vdrift_local[2] = 0.0;
355
356	tim	1127	for(vd = 0; vd < nobj; vd++){
357	chuckv	401
358	tim	1127	amass = info->integrableObjects[vd]->getMass();
359			info->integrableObjects[vd]->getVel( aVel );
360	gezelter	608
361			for(j = 0; j < 3; j++)
362			vdrift_local[j] += aVel[j] * amass;
363	chuckv	401
364	gezelter	608	mtot_local += amass;
365	chuckv	401	}
366
367			#ifdef IS_MPI
368	mmeineke	447	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
369			MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
370	chuckv	401	#else
371			mtot = mtot_local;
372			for(vd = 0; vd < 3; vd++) {
373			vdrift[vd] = vdrift_local[vd];
374			}
375			#endif
376
377			for (vd = 0; vd < 3; vd++) {
378			vdrift[vd] = vdrift[vd] / mtot;
379			}
380
381			}
382
383	tim	763	void Thermo::getCOM(double COM[3]){
384
385			double mtot, mtot_local;
386			double aPos[3], amass;
387			double COM_local[3];
388	tim	1127	int i, j;
389			int nobj;
390	tim	763
391			mtot_local = 0.0;
392			COM_local[0] = 0.0;
393			COM_local[1] = 0.0;
394			COM_local[2] = 0.0;
395	tim	1127
396			nobj = info->integrableObjects.size();
397			for(i = 0; i < nobj; i++){
398	tim	763
399	tim	1127	amass = info->integrableObjects[i]->getMass();
400			info->integrableObjects[i]->getPos( aPos );
401	tim	763
402			for(j = 0; j < 3; j++)
403			COM_local[j] += aPos[j] * amass;
404
405			mtot_local += amass;
406			}
407
408			#ifdef IS_MPI
409			MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
410			MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
411			#else
412			mtot = mtot_local;
413			for(i = 0; i < 3; i++) {
414			COM[i] = COM_local[i];
415			}
416			#endif
417
418			for (i = 0; i < 3; i++) {
419			COM[i] = COM[i] / mtot;
420			}
421			}
422	tim	1127
423			void Thermo::removeCOMdrift() {
424			double vdrift[3], aVel[3];
425			int vd, j, nobj;
426
427			nobj = info->integrableObjects.size();
428
429			// Get the Center of Mass drift velocity.
430
431			getCOMVel(vdrift);
432
433			// Corrects for the center of mass drift.
434			// sums all the momentum and divides by total mass.
435
436			for(vd = 0; vd < nobj; vd++){
437
438			info->integrableObjects[vd]->getVel(aVel);
439
440			for (j=0; j < 3; j++)
441			aVel[j] -= vdrift[j];
442
443			info->integrableObjects[vd]->setVel( aVel );
444			}
445	gezelter	1133	}
446	tim	1452
447			void Thermo::removeAngularMomentum(){
448			Vector3d vcom;
449			Vector3d qcom;
450			Vector3d pos;
451			Vector3d vel;
452			double mass;
453			double xx;
454			double yy;
455			double zz;
456			double xy;
457			double xz;
458			double yz;
459			Vector3d localAngMom;
460			Vector3d angMom;
461			Vector3d omega;
462			vector<StuntDouble *> integrableObjects;
463			double localInertiaVec[9];
464			double inertiaVec[9];
465			vector<Vector3d> qMinusQCom;
466			vector<Vector3d> vMinusVCom;
467			Mat3x3d inertiaMat;
468			Mat3x3d inverseInertiaMat;
469
470			integrableObjects = info->integrableObjects;
471			qMinusQCom.resize(integrableObjects.size());
472			vMinusVCom.resize(integrableObjects.size());
473
474			getCOM(qcom.vec);
475			getCOMVel(vcom.vec);
476
477			//initialize components for inertia tensor
478			xx = 0.0;
479			yy = 0.0;
480			zz = 0.0;
481			xy = 0.0;
482			xz = 0.0;
483			yz = 0.0;
484
485			//build components of Inertia tensor
486			//
487			// [ Ixx -Ixy -Ixz ]
488			// J = \| -Iyx Iyy -Iyz \|
489			// [ -Izx -Iyz Izz ]
490			//See Fowles and Cassidy Chapter 9 or Goldstein Chapter 5
491			for(size_t i = 0; i < integrableObjects.size(); i++){
492			integrableObjects[i]->getPos(pos.vec);
493			integrableObjects[i]->getVel(vel.vec);
494			mass = integrableObjects[i]->getMass();
495
496			qMinusQCom[i] = pos - qcom;
497			info->wrapVector(qMinusQCom[i].vec);
498
499			vMinusVCom[i] = vel - vcom;
500
501			//compute moment of inertia coefficents
502			xx += qMinusQCom[i].x * qMinusQCom[i].x * mass;
503			yy += qMinusQCom[i].y * qMinusQCom[i].y * mass;
504			zz += qMinusQCom[i].z * qMinusQCom[i].z * mass;
505
506			// compute products of inertia
507			xy += qMinusQCom[i].x * qMinusQCom[i].y * mass;
508			xz += qMinusQCom[i].x * qMinusQCom[i].z * mass;
509			yz += qMinusQCom[i].y * qMinusQCom[i].z * mass;
510
511			localAngMom += crossProduct(qMinusQCom[i] , vMinusVCom[i] ) * mass;
512
513			}
514
515			localInertiaVec[0] =yy+zz;
516			localInertiaVec[1] = -xy;
517			localInertiaVec[2] = -xz;
518			localInertiaVec[3] = -xy;
519			localInertiaVec[4] = xx+zz;
520			localInertiaVec[5] = -yz;
521			localInertiaVec[6] = -xz;
522			localInertiaVec[7] = -yz;
523			localInertiaVec[8] = xx+yy;
524
525			//Sum and distribute inertia and angmom arrays
526			#ifdef MPI
527
528			MPI_Allreduce(localInertiaVec, inertiaVec, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
529
530			MPI_Allreduce(localAngMom.vec, angMom.vec, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
531
532			inertiaMat.element[0][0] = inertiaVec[0];
533			inertiaMat.element[0][1] = inertiaVec[1];
534			inertiaMat.element[0][2] = inertiaVec[2];
535
536			inertiaMat.element[1][0] = inertiaVec[3];
537			inertiaMat.element[1][1] = inertiaVec[4];
538			inertiaMat.element[1][2] = inertiaVec[5];
539
540			inertiaMat.element[2][0] = inertiaVec[6];
541			inertiaMat.element[2][1] = inertiaVec[7];
542			inertiaMat.element[2][2] = inertiaVec[8];
543
544			#else
545
546			inertiaMat.element[0][0] = localInertiaVec[0];
547			inertiaMat.element[0][1] = localInertiaVec[1];
548			inertiaMat.element[0][2] = localInertiaVec[2];
549
550			inertiaMat.element[1][0] = localInertiaVec[3];
551			inertiaMat.element[1][1] = localInertiaVec[4];
552			inertiaMat.element[1][2] = localInertiaVec[5];
553
554			inertiaMat.element[2][0] = localInertiaVec[6];
555			inertiaMat.element[2][1] = localInertiaVec[7];
556			inertiaMat.element[2][2] = localInertiaVec[8];
557
558			angMom = localAngMom;
559			#endif
560
561			//invert the moment of inertia tensor by LU-decomposition / backsolving:
562
563			inverseInertiaMat = inertiaMat.inverse();
564
565			//calculate the angular velocities: omega = I^-1 . L
566
567			omega = inverseInertiaMat * angMom;
568
569			//subtract out center of mass velocity and angular momentum from
570			//particle velocities
571
572			for(size_t i = 0; i < integrableObjects.size(); i++){
573			vel = vMinusVCom[i] - crossProduct(omega, qMinusQCom[i]);
574			integrableObjects[i]->setVel(vel.vec);
575			}
576			}
577
578			double Thermo::getConsEnergy(){
579			ConstraintPair* consPair;
580			double totConsEnergy;
581			double bondLen2;
582			double dist;
583			double lamda;
584
585			totConsEnergy = 0;
586
587			for(cpIter->first(); !cpIter->isEnd(); cpIter->next()){
588			consPair = cpIter->currentItem();
589			bondLen2 = consPair->getBondLength2();
590			lamda = consPair->getLamda();
591			//dist = consPair->getDistance();
592
593			//totConsEnergy += lamda * (dist*dist - bondLen2);
594			}
595
596			return totConsEnergy;
597			}
598
599