ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Thermo.cpp
Revision: 1192
Committed: Mon May 24 21:03:30 2004 UTC (20 years, 1 month ago) by gezelter
File size: 10131 byte(s)
Log Message:
Fixes for stress / pressure tensor by cutoff group

File Contents

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