ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/brains/Thermo.cpp
(Generate patch)

Comparing:
trunk/OOPSE-3.0/src/brains/Thermo.cpp (file contents), Revision 1490 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/brains/Thermo.cpp (file contents), Revision 1725 by tim, Wed Nov 10 22:01:06 2004 UTC

# Line 1 | Line 1
1   #include <math.h>
2   #include <iostream>
3 +
4   using namespace std;
5  
6   #ifdef IS_MPI
7 +
8   #include <mpi.h>
9   #endif //is_mpi
10  
11 < #include "Thermo.hpp"
12 < #include "SRI.hpp"
13 < #include "Integrator.hpp"
14 < #include "simError.h"
15 < #include "MatVec3.h"
11 > #include "brains/Thermo.hpp"
12 > #include "primitives/SRI.hpp"
13 > #include "integrators/Integrator.hpp"
14 > #include "utils/simError.h"
15 > #include "math/MatVec3.h"
16  
17   #ifdef IS_MPI
18 +
19   #define __C
20 < #include "mpiSimulation.hpp"
20 > #include "brains/mpiSimulation.hpp"
21   #endif // is_mpi
22  
23 < inline double roundMe( double x ){
24 <          return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
22 < }
23 > namespace oopse {
24 > Thermo::Thermo(SimInfo *the_info) { info = the_info; }
25  
26 < Thermo::Thermo( SimInfo* the_info ) {
25 <  info = the_info;
26 <  int baseSeed = the_info->getSeed();
27 <  
28 <  gaussStream = new gaussianSPRNG( baseSeed );
29 < }
26 > Thermo::~Thermo() { }
27  
28 < Thermo::~Thermo(){
29 <  delete gaussStream;
30 < }
28 > double Thermo::getKinetic() {
29 >    const double         e_convert =
30 >                             4.184E - 4; // convert kcal/mol -> (amu A^2)/fs^2
31 >          double         kinetic;
32 >          double         amass;
33 >          double         aVel[3],
34 >                         aJ[3],
35 >                         I[3][3];
36 >          int            i,
37 >                         j,
38 >                         k,
39 >                         kl;
40  
41 < double Thermo::getKinetic(){
41 >          double         kinetic_global;
42 >    vector<StuntDouble *>integrableObjects = info->integrableObjects;
43  
44 <  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
45 <  double kinetic;
39 <  double amass;
40 <  double aVel[3], aJ[3], I[3][3];
41 <  int i, j, k, kl;
44 >    kinetic = 0.0;
45 >    kinetic_global = 0.0;
46  
47 <  double kinetic_global;
48 <  vector<StuntDouble *> integrableObjects = info->integrableObjects;
49 <  
46 <  kinetic = 0.0;
47 <  kinetic_global = 0.0;
47 >    for( kl = 0; kl < integrableObjects.size(); kl++ ) {
48 >        aVel = integrableObjects[kl]->getVel();
49 >        amass = integrableObjects[kl]->getMass();
50  
51 <  for (kl=0; kl<integrableObjects.size(); kl++) {
52 <    integrableObjects[kl]->getVel(aVel);
51 <    amass = integrableObjects[kl]->getMass();
51 >        for( j = 0; j < 3; j++ )
52 >        kinetic += amass * aVel[j] * aVel[j];
53  
54 <   for(j=0; j<3; j++)
55 <      kinetic += amass*aVel[j]*aVel[j];
54 >        if (integrableObjects[kl]->isDirectional()) {
55 >            aJ = integrableObjects[kl]->getJ();
56 >            integrableObjects[kl]->getI(I);
57  
58 <   if (integrableObjects[kl]->isDirectional()){
59 <
60 <      integrableObjects[kl]->getJ( aJ );
61 <      integrableObjects[kl]->getI( I );
58 >            if (integrableObjects[kl]->isLinear()) {
59 >                i = integrableObjects[kl]->linearAxis();
60 >                j = (i + 1) % 3;
61 >                k = (i + 2) % 3;
62 >                kinetic += aJ[j] * aJ[j] / I[j][j] + aJ[k] * aJ[k] / I[k][k];
63 >            } else {
64 >                for( j = 0; j < 3; j++ )
65 >                kinetic += aJ[j] * aJ[j] / I[j][j];
66 >            }
67 >        }
68 >    }
69  
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 –  }
70   #ifdef IS_MPI
71 <  MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
72 <                MPI_SUM, MPI_COMM_WORLD);
73 <  kinetic = kinetic_global;
71 >
72 >    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_DOUBLE, MPI_SUM,
73 >                  MPI_COMM_WORLD);
74 >    kinetic = kinetic_global;
75 >
76   #endif //is_mpi
77 –  
78 –  kinetic = kinetic * 0.5 / e_convert;
77  
78 <  return kinetic;
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;
83 > double Thermo::getPotential() {
84 >    double    potential_local;
85 >    double    potential;
86 >    int       el,
87 >              nSRI;
88 >    Molecule *molecules;
89  
90 <  molecules = info->molecules;
91 <  nSRI = info->n_SRI;
90 >    molecules = info->molecules;
91 >    nSRI = info->n_SRI;
92  
93 <  potential_local = 0.0;
94 <  potential = 0.0;
95 <  potential_local += info->lrPot;
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 <  }
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.
101 >    // Get total potential for entire system from MPI.
102 >
103   #ifdef IS_MPI
104 <  MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
105 <                MPI_SUM, MPI_COMM_WORLD);
104 >
105 >    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
106 >                  MPI_COMM_WORLD);
107 >
108   #else
109 <  potential = potential_local;
109 >
110 >    potential = potential_local;
111 >
112   #endif // is_mpi
113  
114 <  return potential;
114 >    return potential;
115   }
116  
117 < double Thermo::getTotalE(){
117 > double Thermo::getTotalE() {
118 >    double total;
119  
120 <  double total;
121 <
116 <  total = this->getKinetic() + this->getPotential();
117 <  return total;
120 >    total = this->getKinetic() + this->getPotential();
121 >    return total;
122   }
123  
124 < double Thermo::getTemperature(){
124 > double Thermo::getTemperature() {
125 >    const double kb = 1.9872156E - 3; // boltzman's constant in kcal/(mol K)
126 >          double temperature;
127  
128 <  const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
123 <  double temperature;
128 >    temperature
129  
130 <  temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
126 <  return temperature;
127 < }
130 >    = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
131  
132 < double Thermo::getVolume() {
130 <
131 <  return info->boxVol;
132 >    return temperature;
133   }
134  
135 + double Thermo::getVolume() { return info->boxVol; }
136 +
137   double Thermo::getPressure() {
138  
139 <  // 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;
139 >    // Relies on the calculation of the full molecular pressure tensor
140  
141 <  this->getPressureTensor(press);
141 >    const double p_convert = 1.63882576e8;
142 >          double press[3][3];
143 >          double pressure;
144  
145 <  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
145 >    this->getPressureTensor(press);
146  
147 <  return pressure;
147 >    pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
148 >
149 >    return pressure;
150   }
151  
152   double Thermo::getPressureX() {
153  
154 <  // 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;
154 >    // Relies on the calculation of the full molecular pressure tensor
155  
156 <  this->getPressureTensor(press);
156 >    const double p_convert = 1.63882576e8;
157 >          double press[3][3];
158 >          double pressureX;
159  
160 <  pressureX = p_convert * press[0][0];
160 >    this->getPressureTensor(press);
161  
162 <  return pressureX;
162 >    pressureX = p_convert * press[0][0];
163 >
164 >    return pressureX;
165   }
166  
167   double Thermo::getPressureY() {
168  
169 <  // 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;
169 >    // Relies on the calculation of the full molecular pressure tensor
170  
171 <  this->getPressureTensor(press);
171 >    const double p_convert = 1.63882576e8;
172 >          double press[3][3];
173 >          double pressureY;
174  
175 <  pressureY = p_convert * press[1][1];
175 >    this->getPressureTensor(press);
176  
177 <  return pressureY;
177 >    pressureY = p_convert * press[1][1];
178 >
179 >    return pressureY;
180   }
181  
182   double Thermo::getPressureZ() {
183  
184 <  // 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;
184 >    // Relies on the calculation of the full molecular pressure tensor
185  
186 <  this->getPressureTensor(press);
186 >    const double p_convert = 1.63882576e8;
187 >          double press[3][3];
188 >          double pressureZ;
189  
190 <  pressureZ = p_convert * press[2][2];
190 >    this->getPressureTensor(press);
191  
192 <  return pressureZ;
192 >    pressureZ = p_convert * press[2][2];
193 >
194 >    return pressureZ;
195   }
196  
197 + void Thermo::getPressureTensor(double press[3][3]) {
198 +    // returns pressure tensor in units amu*fs^-2*Ang^-1
199 +    // routine derived via viral theorem description in:
200 +    // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
201  
202 < 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
202 >    const double e_convert = 4.184e - 4;
203  
204 <  const double e_convert = 4.184e-4;
204 >          double molmass,
205 >                 volume;
206 >          double vcom[3];
207 >          double p_local[9],
208 >                 p_global[9];
209 >          int    i,
210 >                 j,
211 >                 k;
212  
213 <  double molmass, volume;
214 <  double vcom[3];
215 <  double p_local[9], p_global[9];
216 <  int i, j, k;
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 < 9; i++) {    
208 <    p_local[i] = 0.0;
209 <    p_global[i] = 0.0;
210 <  }
211 <
212 <  // use velocities of integrableObjects and their masses:  
218 >    // use velocities of integrableObjects and their masses:  
219  
220 <  for (i=0; i < info->integrableObjects.size(); i++) {
220 >    for( i = 0; i < info->integrableObjects.size(); i++ ) {
221 >        molmass = info->integrableObjects[i]->getMass();
222  
223 <    molmass = info->integrableObjects[i]->getMass();
217 <    
218 <    info->integrableObjects[i]->getVel(vcom);
219 <    
220 <    p_local[0] += molmass * (vcom[0] * vcom[0]);
221 <    p_local[1] += molmass * (vcom[0] * vcom[1]);
222 <    p_local[2] += molmass * (vcom[0] * vcom[2]);
223 <    p_local[3] += molmass * (vcom[1] * vcom[0]);
224 <    p_local[4] += molmass * (vcom[1] * vcom[1]);
225 <    p_local[5] += molmass * (vcom[1] * vcom[2]);
226 <    p_local[6] += molmass * (vcom[2] * vcom[0]);
227 <    p_local[7] += molmass * (vcom[2] * vcom[1]);
228 <    p_local[8] += molmass * (vcom[2] * vcom[2]);
223 >        vcom = info->integrableObjects[i]->getVel();
224  
225 <  }
226 <
227 <  // Get total for entire system from MPI.
228 <  
229 < #ifdef IS_MPI
230 <  MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
231 < #else
232 <  for (i=0; i<9; i++) {
233 <    p_global[i] = p_local[i];
239 <  }
240 < #endif // is_mpi
241 <
242 <  volume = this->getVolume();
243 <
244 <
245 <
246 <  for(i = 0; i < 3; i++) {
247 <    for (j = 0; j < 3; j++) {
248 <      k = 3*i + j;
249 <      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
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      }
251 –  }
252 – }
235  
236 < void Thermo::velocitize() {
255 <  
256 <  double aVel[3], aJ[3], I[3][3];
257 <  int i, j, l, m, n, vr, vd; // velocity randomizer loop counters
258 <  double vdrift[3];
259 <  double vbar;
260 <  const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
261 <  double av2;
262 <  double kebar;
263 <  double temperature;
264 <  int nobj;
236 >    // Get total for entire system from MPI.
237  
238 <  if (!info->have_target_temp) {
267 <    sprintf( painCave.errMsg,
268 <             "You can't resample the velocities without a targetTemp!\n"
269 <             );
270 <    painCave.isFatal = 1;
271 <    painCave.severity = OOPSE_ERROR;
272 <    simError();
273 <    return;
274 <  }
238 > #ifdef IS_MPI
239  
240 <  nobj = info->integrableObjects.size();
277 <  
278 <  temperature   = info->target_temp;
279 <  
280 <  kebar = kb * temperature * (double)info->ndfRaw /
281 <    ( 2.0 * (double)info->ndf );
282 <  
283 <  for(vr = 0; vr < nobj; vr++){
284 <    
285 <    // uses equipartition theory to solve for vbar in angstrom/fs
240 >    MPI_Allreduce(p_local, p_global, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
241  
287 –    av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
288 –    vbar = sqrt( av2 );
289 –
290 –    // picks random velocities from a gaussian distribution
291 –    // centered on vbar
292 –
293 –    for (j=0; j<3; j++)
294 –      aVel[j] = vbar * gaussStream->getGaussian();
295 –    
296 –    info->integrableObjects[vr]->setVel( aVel );
297 –    
298 –    if(info->integrableObjects[vr]->isDirectional()){
299 –
300 –      info->integrableObjects[vr]->getI( I );
301 –
302 –      if (info->integrableObjects[vr]->isLinear()) {
303 –
304 –        l= info->integrableObjects[vr]->linearAxis();
305 –        m = (l+1)%3;
306 –        n = (l+2)%3;
307 –
308 –        aJ[l] = 0.0;
309 –        vbar = sqrt( 2.0 * kebar * I[m][m] );
310 –        aJ[m] = vbar * gaussStream->getGaussian();
311 –        vbar = sqrt( 2.0 * kebar * I[n][n] );
312 –        aJ[n] = vbar * gaussStream->getGaussian();
313 –        
314 –      } else {
315 –        for (j = 0 ; j < 3; j++) {
316 –          vbar = sqrt( 2.0 * kebar * I[j][j] );
317 –          aJ[j] = vbar * gaussStream->getGaussian();
318 –        }      
319 –      } // else isLinear
320 –
321 –      info->integrableObjects[vr]->setJ( aJ );
322 –      
323 –    }//isDirectional
324 –
325 –  }
326 –
327 –  // Get the Center of Mass drift velocity.
328 –
329 –  getCOMVel(vdrift);
330 –  
331 –  //  Corrects for the center of mass drift.
332 –  // sums all the momentum and divides by total mass.
333 –
334 –  for(vd = 0; vd < nobj; vd++){
335 –    
336 –    info->integrableObjects[vd]->getVel(aVel);
337 –    
338 –    for (j=0; j < 3; j++)
339 –      aVel[j] -= vdrift[j];
340 –        
341 –    info->integrableObjects[vd]->setVel( aVel );
342 –  }
343 –
344 – }
345 –
346 – void Thermo::getCOMVel(double vdrift[3]){
347 –
348 –  double mtot, mtot_local;
349 –  double aVel[3], amass;
350 –  double vdrift_local[3];
351 –  int vd, j;
352 –  int nobj;
353 –
354 –  nobj   = info->integrableObjects.size();
355 –
356 –  mtot_local = 0.0;
357 –  vdrift_local[0] = 0.0;
358 –  vdrift_local[1] = 0.0;
359 –  vdrift_local[2] = 0.0;
360 –  
361 –  for(vd = 0; vd < nobj; vd++){
362 –    
363 –    amass = info->integrableObjects[vd]->getMass();
364 –    info->integrableObjects[vd]->getVel( aVel );
365 –
366 –    for(j = 0; j < 3; j++)
367 –      vdrift_local[j] += aVel[j] * amass;
368 –    
369 –    mtot_local += amass;
370 –  }
371 –
372 – #ifdef IS_MPI
373 –  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
374 –  MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
242   #else
376 –  mtot = mtot_local;
377 –  for(vd = 0; vd < 3; vd++) {
378 –    vdrift[vd] = vdrift_local[vd];
379 –  }
380 – #endif
381 –    
382 –  for (vd = 0; vd < 3; vd++) {
383 –    vdrift[vd] = vdrift[vd] / mtot;
384 –  }
385 –  
386 – }
243  
244 < void Thermo::getCOM(double COM[3]){
244 >    for( i = 0; i < 9; i++ ) {
245 >        p_global[i] = p_local[i];
246 >    }
247  
248 <  double mtot, mtot_local;
391 <  double aPos[3], amass;
392 <  double COM_local[3];
393 <  int i, j;
394 <  int nobj;
248 > #endif // is_mpi
249  
250 <  mtot_local = 0.0;
397 <  COM_local[0] = 0.0;
398 <  COM_local[1] = 0.0;
399 <  COM_local[2] = 0.0;
250 >    volume = this->getVolume();
251  
252 <  nobj = info->integrableObjects.size();
253 <  for(i = 0; i < nobj; i++){
254 <    
404 <    amass = info->integrableObjects[i]->getMass();
405 <    info->integrableObjects[i]->getPos( aPos );
252 >    for( i = 0; i < 3; i++ ) {
253 >        for( j = 0; j < 3; j++ ) {
254 >            k = 3 * i + j;
255  
256 <    for(j = 0; j < 3; j++)
408 <      COM_local[j] += aPos[j] * amass;
409 <    
410 <    mtot_local += amass;
411 <  }
256 >            press
257  
258 < #ifdef IS_MPI
259 <  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
260 <  MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
416 < #else
417 <  mtot = mtot_local;
418 <  for(i = 0; i < 3; i++) {
419 <    COM[i] = COM_local[i];
420 <  }
421 < #endif
422 <    
423 <  for (i = 0; i < 3; i++) {
424 <    COM[i] = COM[i] / mtot;
425 <  }
258 >            [i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
259 >        }
260 >    }
261   }
262 <
428 < void Thermo::removeCOMdrift() {
429 <  double vdrift[3], aVel[3];
430 <  int vd, j, nobj;
431 <
432 <  nobj = info->integrableObjects.size();
433 <
434 <  // Get the Center of Mass drift velocity.
435 <
436 <  getCOMVel(vdrift);
437 <  
438 <  //  Corrects for the center of mass drift.
439 <  // sums all the momentum and divides by total mass.
440 <
441 <  for(vd = 0; vd < nobj; vd++){
442 <    
443 <    info->integrableObjects[vd]->getVel(aVel);
444 <    
445 <    for (j=0; j < 3; j++)
446 <      aVel[j] -= vdrift[j];
447 <        
448 <    info->integrableObjects[vd]->setVel( aVel );
449 <  }
450 < }
262 > } //end namespace oopse

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines