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


#define BASE_SEED 123456789

Thermo::Thermo( SimInfo* the_entry_plug ) { 
  entry_plug = the_entry_plug;
  int baseSeed = BASE_SEED;
  
  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 = entry_plug->n_atoms;
  atoms = entry_plug->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 = entry_plug->molecules;
  nSRI = entry_plug->n_SRI;

  potential_local = 0.0;
  potential = 0.0;
  potential_local += entry_plug->lrPot;

  for( el=0; el<entry_plug->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)entry_plug->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 entry_plug->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;
}


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, l, nMols;
  Molecule* molecules;

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