ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Thermo.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Thermo.cpp (file contents):
Revision 378 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
Revision 401 by chuckv, Tue Mar 25 22:54:16 2003 UTC

# Line 10 | Line 10 | using namespace std;
10   #include "Thermo.hpp"
11   #include "SRI.hpp"
12   #include "Integrator.hpp"
13 + #define __C
14 + //#include "mpiSimulation.hpp"
15  
16   #define BASE_SEED 123456789
17  
# Line 77 | Line 79 | double Thermo::getPotential(){
79  
80   double Thermo::getPotential(){
81    
82 +  double potential_local;
83    double potential;
81 –  double potential_global;
84    int el, nSRI;
85    SRI** sris;
86  
87    sris = entry_plug->sr_interactions;
88    nSRI = entry_plug->n_SRI;
89  
90 <  potential = 0.0;
91 <  potential_global = 0.0;
90 <  potential += entry_plug->lrPot;
90 >  potential_local = 0.0;
91 >  potential_local += entry_plug->lrPot;
92  
93 <  for( el=0; el<nSRI; el++ ){
94 <    
94 <    potential += sris[el]->get_potential();
93 >  for( el=0; el<nSRI; el++ ){    
94 >    potential_local += sris[el]->get_potential();
95    }
96  
97    // Get total potential for entire system from MPI.
98   #ifdef IS_MPI
99 <  MPI::COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
100 <  potential = potential_global;
101 <
99 >  MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM);
100 > #else
101 >  potential = potential_local;
102   #endif // is_mpi
103  
104    return potential;
# Line 116 | Line 116 | double Thermo::getTemperature(){
116  
117    const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
118    double temperature;
119 +  int ndf_local, ndf;
120    
121 <  int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
122 <    - entry_plug->n_constraints - 3;
121 >  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
122 >    - entry_plug->n_constraints;
123  
124 + #ifdef IS_MPI
125 +  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
126 + #else
127 +  ndf = ndf_local;
128 + #endif
129 +
130 +  ndf = ndf - 3;
131 +  
132    temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
133    return temperature;
134   }
# Line 139 | Line 148 | void Thermo::velocitize() {
148    double vx, vy, vz;
149    double jx, jy, jz;
150    int i, vr, vd; // velocity randomizer loop counters
151 <  double vdrift[3];
143 <  double mtot = 0.0;
151 >  double *vdrift;
152    double vbar;
153    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
154    double av2;
# Line 185 | Line 193 | void Thermo::velocitize() {
193      atoms[vr]->set_vy( vy );
194      atoms[vr]->set_vz( vz );
195    }
196 +
197 +  // Get the Center of Mass drift velocity.
198 +
199 +  vdrift = getCOMVel();
200    
201    //  Corrects for the center of mass drift.
202    // sums all the momentum and divides by total mass.
191 –  
192 –  mtot = 0.0;
193 –  vdrift[0] = 0.0;
194 –  vdrift[1] = 0.0;
195 –  vdrift[2] = 0.0;
196 –  for(vd = 0; vd < n_atoms; vd++){
197 –    
198 –    vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
199 –    vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
200 –    vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
201 –    
202 –    mtot += atoms[vd]->getMass();
203 –  }
204 –  
205 –  for (vd = 0; vd < 3; vd++) {
206 –    vdrift[vd] = vdrift[vd] / mtot;
207 –  }
208 –  
203  
204    for(vd = 0; vd < n_atoms; vd++){
205      
206      vx = atoms[vd]->get_vx();
207      vy = atoms[vd]->get_vy();
208      vz = atoms[vd]->get_vz();
209 <    
216 <    
209 >        
210      vx -= vdrift[0];
211      vy -= vdrift[1];
212      vz -= vdrift[2];
# Line 246 | Line 239 | void Thermo::velocitize() {
239      }  
240    }
241   }
242 +
243 + double* Thermo::getCOMVel(){
244 +
245 +  double mtot, mtot_local;
246 +  double* vdrift;
247 +  double vdrift_local[3];
248 +  int vd, n_atoms;
249 +  Atom** atoms;
250 +
251 +  vdrift = new double[3];
252 +  // We are very careless here with the distinction between n_atoms and n_local
253 +  // We should really fix this before someone pokes an eye out.
254 +
255 +  n_atoms = entry_plug->n_atoms;  
256 +  atoms   = entry_plug->atoms;
257 +
258 +  mtot_local = 0.0;
259 +  vdrift_local[0] = 0.0;
260 +  vdrift_local[1] = 0.0;
261 +  vdrift_local[2] = 0.0;
262 +  
263 +  for(vd = 0; vd < n_atoms; vd++){
264 +    
265 +    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
266 +    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
267 +    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
268 +    
269 +    mtot_local += atoms[vd]->getMass();
270 +  }
271 +
272 + #ifdef IS_MPI
273 +  MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM);
274 +  MPI::COMM_WORLD.Allreduce(&vdrift_local,&vdrift,3,MPI_DOUBLE,MPI_SUM);
275 + #else
276 +  mtot = mtot_local;
277 +  for(vd = 0; vd < 3; vd++) {
278 +    vdrift[vd] = vdrift_local[vd];
279 +  }
280 + #endif
281 +    
282 +  for (vd = 0; vd < 3; vd++) {
283 +    vdrift[vd] = vdrift[vd] / mtot;
284 +  }
285 +  
286 +  return vdrift;
287 + }
288 +

Diff Legend

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