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

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 j, kl;

  DirectionalAtom *dAtom;

  int n_atoms;
  double kinetic_global;
  Atom** atoms;

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

  kinetic = 0.0;
  kinetic_global = 0.0;
  for( kl=0; kl < n_atoms; kl++ ){
    
    atoms[kl]->getVel(aVel);
    amass = atoms[kl]->getMass();
    
    for (j=0; j < 3; j++) 
      kinetic += amass * aVel[j] * aVel[j];

    if( atoms[kl]->isDirectional() ){
            
      dAtom = (DirectionalAtom *)atoms[kl];

      dAtom->getJ( aJ );
      dAtom->getI( I );
      
      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

#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)info->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[3][3];

  u = this->getTotalE();

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

  v = this->getVolume();

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

double Thermo::getVolume() {

  return info->boxVol;
}

double Thermo::getPressure() {

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

  this->getPressureTensor(press);

  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;

  return pressure;
}

double Thermo::getPressureX() {

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

  this->getPressureTensor(press);

  pressureX = p_convert * press[0][0];

  return pressureX;
}

double Thermo::getPressureY() {

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

  this->getPressureTensor(press);

  pressureY = p_convert * press[1][1];

  return pressureY;
}

double Thermo::getPressureZ() {

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

  this->getPressureTensor(press);

  pressureZ = p_convert * press[2][2];

  return pressureZ;
}


void Thermo::getPressureTensor(double press[3][3]){
  // returns pressure tensor in units amu*fs^-2*Ang^-1
  // routine derived via viral theorem description in:
  // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322

  const double e_convert = 4.184e-4;

  double molmass, volume;
  double vcom[3];
  double p_local[9], p_global[9];
  int i, j, k, nMols;
  Molecule* molecules;

  nMols = info->n_mol;
  molecules = info->molecules;
  //tau = info->tau;

  // use velocities of molecular centers of mass and molecular masses:
  for (i=0; i < 9; i++) {    
    p_local[i] = 0.0;
    p_global[i] = 0.0;
  }

  for (i=0; i < nMols; i++) {
    molmass = molecules[i].getCOMvel(vcom);

    p_local[0] += molmass * (vcom[0] * vcom[0]); 
    p_local[1] += molmass * (vcom[0] * vcom[1]); 
    p_local[2] += molmass * (vcom[0] * vcom[2]); 
    p_local[3] += molmass * (vcom[1] * vcom[0]); 
    p_local[4] += molmass * (vcom[1] * vcom[1]); 
    p_local[5] += molmass * (vcom[1] * vcom[2]); 
    p_local[6] += molmass * (vcom[2] * vcom[0]); 
    p_local[7] += molmass * (vcom[2] * vcom[1]); 
    p_local[8] += molmass * (vcom[2] * vcom[2]); 
  }

  // Get total for entire system from MPI.

#ifdef IS_MPI
  MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#else
  for (i=0; i<9; i++) {
    p_global[i] = p_local[i]; 
  }
#endif // is_mpi

  volume = this->getVolume();

  for(i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      k = 3*i + j;
      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;

    }
  }
}

void Thermo::velocitize() {
  
  double x,y;
  double aVel[3], aJ[3], I[3][3];
  int i, j, vr, vd; // velocity randomizer loop counters
  double vdrift[3];
  double vbar;
  const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
  double av2;
  double kebar;
  int n_atoms;
  Atom** atoms;
  DirectionalAtom* dAtom;
  double temperature;
  int n_oriented;
  int n_constraints;

  atoms         = info->atoms;
  n_atoms       = info->n_atoms;
  temperature   = info->target_temp;
  n_oriented    = info->n_oriented;
  n_constraints = info->n_constraints;
  
  kebar = kb * temperature * (double)info->ndf / 
    ( 2.0 * (double)info->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

    for (j=0; j<3; j++) 
      aVel[j] = vbar * gaussStream->getGaussian();
    
    atoms[vr]->setVel( aVel );

  }

  // Get the Center of Mass drift velocity.

  getCOMVel(vdrift);
  
  //  Corrects for the center of mass drift.
  // sums all the momentum and divides by total mass.

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

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

        }       

        dAtom->setJ( aJ );

      }
    }   
  }
}

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

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

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

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

  mtot_local = 0.0;
  vdrift_local[0] = 0.0;
  vdrift_local[1] = 0.0;
  vdrift_local[2] = 0.0;
  
  for(vd = 0; vd < n_atoms; vd++){
    
    amass = atoms[vd]->getMass();
    atoms[vd]->getVel( aVel );

    for(j = 0; j < 3; j++) 
      vdrift_local[j] += aVel[j] * amass;
    
    mtot_local += amass;
  }

#ifdef IS_MPI
  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
#else
  mtot = mtot_local;
  for(vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift_local[vd];
  }
#endif
    
  for (vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift[vd] / mtot;
  }
  
}

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