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

Comparing:
branches/mmeineke/OOPSE/libmdtools/Thermo.cpp (file contents), Revision 377 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
trunk/OOPSE/libmdtools/Thermo.cpp (file contents), Revision 402 by mmeineke, Wed Mar 26 14:55:50 2003 UTC

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

Diff Legend

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