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"

#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.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, 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 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->ndfRaw / 
    ( 2.0 * (double)info->ndf );
  
  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;
  }
  
}

void Thermo::getCOM(double COM[3]){

  double mtot, mtot_local;
  double aPos[3], amass;
  double COM_local[3];
  int i, 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;
  COM_local[0] = 0.0;
  COM_local[1] = 0.0;
  COM_local[2] = 0.0;
  
  for(i = 0; i < n_atoms; i++){
    
    amass = atoms[i]->getMass();
    atoms[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;
  }
}
Revision:	829
Committed:	Tue Oct 28 16:03:37 2003 UTC (20 years, 8 months ago) by gezelter
File size:	9069 byte(s)
Log Message:	replace c++ header stuff with more portable c header stuff Also, mod file fixes and portability changes Some fortran changes will need to be reversed.
#	Content
1	#include <math.h>
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	#include "simError.h"
13
14	#ifdef IS_MPI
15	#define __C
16	#include "mpiSimulation.hpp"
17	#endif // is_mpi
18
19	Thermo::Thermo( SimInfo* the_info ) {
20	info = the_info;
21	int baseSeed = the_info->getSeed();
22
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	double kinetic;
34	double amass;
35	double aVel[3], aJ[3], I[3][3];
36	int j, kl;
37
38	DirectionalAtom *dAtom;
39
40	int n_atoms;
41	double kinetic_global;
42	Atom** atoms;
43
44
45	n_atoms = info->n_atoms;
46	atoms = info->atoms;
47
48	kinetic = 0.0;
49	kinetic_global = 0.0;
50	for( kl=0; kl < n_atoms; kl++ ){
51
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
58	if( atoms[kl]->isDirectional() ){
59
60	dAtom = (DirectionalAtom *)atoms[kl];
61
62	dAtom->getJ( aJ );
63	dAtom->getI( I );
64
65	for (j=0; j<3; j++)
66	kinetic += aJ[j]*aJ[j] / I[j][j];
67
68	}
69	}
70	#ifdef IS_MPI
71	MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
72	MPI_SUM, MPI_COMM_WORLD);
73	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	double potential_local;
84	double potential;
85	int el, nSRI;
86	Molecule* molecules;
87
88	molecules = info->molecules;
89	nSRI = info->n_SRI;
90
91	potential_local = 0.0;
92	potential = 0.0;
93	potential_local += info->lrPot;
94
95	for( el=0; el<info->n_mol; el++ ){
96	potential_local += molecules[el].getPotential();
97	}
98
99	// Get total potential for entire system from MPI.
100	#ifdef IS_MPI
101	MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
102	MPI_SUM, MPI_COMM_WORLD);
103	#else
104	potential = potential_local;
105	#endif // is_mpi
106
107	#ifdef IS_MPI
108	/*
109	std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
110	*/
111	#endif
112
113	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	double Thermo::getTemperature(){
125
126	const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
127	double temperature;
128
129	temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
130	return temperature;
131	}
132
133	double Thermo::getVolume() {
134
135	return info->boxVol;
136	}
137
138	double Thermo::getPressure() {
139
140	// Relies on the calculation of the full molecular pressure tensor
141
142	const double p_convert = 1.63882576e8;
143	double press[3][3];
144	double pressure;
145
146	this->getPressureTensor(press);
147
148	pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
149
150	return pressure;
151	}
152
153	double Thermo::getPressureX() {
154
155	// Relies on the calculation of the full molecular pressure tensor
156
157	const double p_convert = 1.63882576e8;
158	double press[3][3];
159	double pressureX;
160
161	this->getPressureTensor(press);
162
163	pressureX = p_convert * press[0][0];
164
165	return pressureX;
166	}
167
168	double Thermo::getPressureY() {
169
170	// Relies on the calculation of the full molecular pressure tensor
171
172	const double p_convert = 1.63882576e8;
173	double press[3][3];
174	double pressureY;
175
176	this->getPressureTensor(press);
177
178	pressureY = p_convert * press[1][1];
179
180	return pressureY;
181	}
182
183	double Thermo::getPressureZ() {
184
185	// Relies on the calculation of the full molecular pressure tensor
186
187	const double p_convert = 1.63882576e8;
188	double press[3][3];
189	double pressureZ;
190
191	this->getPressureTensor(press);
192
193	pressureZ = p_convert * press[2][2];
194
195	return pressureZ;
196	}
197
198
199	void Thermo::getPressureTensor(double press[3][3]){
200	// returns pressure tensor in units amufs^-2Ang^-1
201	// routine derived via viral theorem description in:
202	// Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
203
204	const double e_convert = 4.184e-4;
205
206	double molmass, volume;
207	double vcom[3];
208	double p_local[9], p_global[9];
209	int i, j, k, nMols;
210	Molecule* molecules;
211
212	nMols = info->n_mol;
213	molecules = info->molecules;
214	//tau = info->tau;
215
216	// use velocities of molecular centers of mass and molecular masses:
217	for (i=0; i < 9; i++) {
218	p_local[i] = 0.0;
219	p_global[i] = 0.0;
220	}
221
222	for (i=0; i < nMols; i++) {
223	molmass = molecules[i].getCOMvel(vcom);
224
225	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	}
235
236	// Get total for entire system from MPI.
237
238	#ifdef IS_MPI
239	MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
240	#else
241	for (i=0; i<9; i++) {
242	p_global[i] = p_local[i];
243	}
244	#endif // is_mpi
245
246	volume = this->getVolume();
247
248	for(i = 0; i < 3; i++) {
249	for (j = 0; j < 3; j++) {
250	k = 3*i + j;
251	press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
252
253	}
254	}
255	}
256
257	void Thermo::velocitize() {
258
259	double aVel[3], aJ[3], I[3][3];
260	int i, j, vr, vd; // velocity randomizer loop counters
261	double vdrift[3];
262	double vbar;
263	const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
264	double av2;
265	double kebar;
266	int n_atoms;
267	Atom** atoms;
268	DirectionalAtom* dAtom;
269	double temperature;
270	int n_oriented;
271	int n_constraints;
272
273	atoms = info->atoms;
274	n_atoms = info->n_atoms;
275	temperature = info->target_temp;
276	n_oriented = info->n_oriented;
277	n_constraints = info->n_constraints;
278
279	kebar = kb * temperature * (double)info->ndfRaw /
280	( 2.0 * (double)info->ndf );
281
282	for(vr = 0; vr < n_atoms; vr++){
283
284	// uses equipartition theory to solve for vbar in angstrom/fs
285
286	av2 = 2.0 * kebar / atoms[vr]->getMass();
287	vbar = sqrt( av2 );
288
289	// vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
290
291	// picks random velocities from a gaussian distribution
292	// centered on vbar
293
294	for (j=0; j<3; j++)
295	aVel[j] = vbar * gaussStream->getGaussian();
296
297	atoms[vr]->setVel( aVel );
298
299	}
300
301	// Get the Center of Mass drift velocity.
302
303	getCOMVel(vdrift);
304
305	// Corrects for the center of mass drift.
306	// sums all the momentum and divides by total mass.
307
308	for(vd = 0; vd < n_atoms; vd++){
309
310	atoms[vd]->getVel(aVel);
311
312	for (j=0; j < 3; j++)
313	aVel[j] -= vdrift[j];
314
315	atoms[vd]->setVel( aVel );
316	}
317	if( n_oriented ){
318
319	for( i=0; i<n_atoms; i++ ){
320
321	if( atoms[i]->isDirectional() ){
322
323	dAtom = (DirectionalAtom *)atoms[i];
324	dAtom->getI( I );
325
326	for (j = 0 ; j < 3; j++) {
327
328	vbar = sqrt( 2.0 * kebar * I[j][j] );
329	aJ[j] = vbar * gaussStream->getGaussian();
330
331	}
332
333	dAtom->setJ( aJ );
334
335	}
336	}
337	}
338	}
339
340	void Thermo::getCOMVel(double vdrift[3]){
341
342	double mtot, mtot_local;
343	double aVel[3], amass;
344	double vdrift_local[3];
345	int vd, n_atoms, j;
346	Atom** atoms;
347
348	// We are very careless here with the distinction between n_atoms and n_local
349	// We should really fix this before someone pokes an eye out.
350
351	n_atoms = info->n_atoms;
352	atoms = info->atoms;
353
354	mtot_local = 0.0;
355	vdrift_local[0] = 0.0;
356	vdrift_local[1] = 0.0;
357	vdrift_local[2] = 0.0;
358
359	for(vd = 0; vd < n_atoms; vd++){
360
361	amass = atoms[vd]->getMass();
362	atoms[vd]->getVel( aVel );
363
364	for(j = 0; j < 3; j++)
365	vdrift_local[j] += aVel[j] * amass;
366
367	mtot_local += amass;
368	}
369
370	#ifdef IS_MPI
371	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
372	MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
373	#else
374	mtot = mtot_local;
375	for(vd = 0; vd < 3; vd++) {
376	vdrift[vd] = vdrift_local[vd];
377	}
378	#endif
379
380	for (vd = 0; vd < 3; vd++) {
381	vdrift[vd] = vdrift[vd] / mtot;
382	}
383
384	}
385
386	void Thermo::getCOM(double COM[3]){
387
388	double mtot, mtot_local;
389	double aPos[3], amass;
390	double COM_local[3];
391	int i, n_atoms, j;
392	Atom** atoms;
393
394	// We are very careless here with the distinction between n_atoms and n_local
395	// We should really fix this before someone pokes an eye out.
396
397	n_atoms = info->n_atoms;
398	atoms = info->atoms;
399
400	mtot_local = 0.0;
401	COM_local[0] = 0.0;
402	COM_local[1] = 0.0;
403	COM_local[2] = 0.0;
404
405	for(i = 0; i < n_atoms; i++){
406
407	amass = atoms[i]->getMass();
408	atoms[i]->getPos( aPos );
409
410	for(j = 0; j < 3; j++)
411	COM_local[j] += aPos[j] * amass;
412
413	mtot_local += amass;
414	}
415
416	#ifdef IS_MPI
417	MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
418	MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
419	#else
420	mtot = mtot_local;
421	for(i = 0; i < 3; i++) {
422	COM[i] = COM_local[i];
423	}
424	#endif
425
426	for (i = 0; i < 3; i++) {
427	COM[i] = COM[i] / mtot;
428	}
429	}