ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/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 423 by mmeineke, Thu Mar 27 20:12:15 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<entry_plug->n_mol; el++ ){    
98 >    potential_local += entry_plug->molecules[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 140 | Line 153 | void Thermo::velocitize() {
153    double jx, jy, jz;
154    int i, vr, vd; // velocity randomizer loop counters
155    double vdrift[3];
143 –  double mtot = 0.0;
156    double vbar;
157    const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
158    double av2;
159    double kebar;
160 <  int ndf; // number of degrees of freedom
161 <  int ndfRaw; // the raw number of degrees of freedom
160 >  int ndf, ndf_local; // number of degrees of freedom
161 >  int ndfRaw, ndfRaw_local; // the raw number of degrees of freedom
162    int n_atoms;
163    Atom** atoms;
164    DirectionalAtom* dAtom;
# Line 160 | Line 172 | void Thermo::velocitize() {
172    n_oriented    = entry_plug->n_oriented;
173    n_constraints = entry_plug->n_constraints;
174    
175 +  // Raw degrees of freedom that we have to set
176 +  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
177  
178 <  ndfRaw = 3 * n_atoms + 3 * n_oriented;
179 <  ndf = ndfRaw - n_constraints - 3;
178 >  // Degrees of freedom that can contain kinetic energy
179 >  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
180 >    - entry_plug->n_constraints;
181 >  
182 > #ifdef IS_MPI
183 >  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
184 >  MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM);
185 > #else
186 >  ndfRaw = ndfRaw_local;
187 >  ndf = ndf_local;
188 > #endif
189 >  ndf = ndf - 3;
190 >
191    kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
192    
193    for(vr = 0; vr < n_atoms; vr++){
# Line 185 | Line 210 | void Thermo::velocitize() {
210      atoms[vr]->set_vy( vy );
211      atoms[vr]->set_vz( vz );
212    }
213 +
214 +  // Get the Center of Mass drift velocity.
215 +
216 +  getCOMVel(vdrift);
217    
218    //  Corrects for the center of mass drift.
219    // 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 –  
220  
221    for(vd = 0; vd < n_atoms; vd++){
222      
223      vx = atoms[vd]->get_vx();
224      vy = atoms[vd]->get_vy();
225      vz = atoms[vd]->get_vz();
226 <    
216 <    
226 >        
227      vx -= vdrift[0];
228      vy -= vdrift[1];
229      vz -= vdrift[2];
# Line 246 | Line 256 | void Thermo::velocitize() {
256      }  
257    }
258   }
259 +
260 + void Thermo::getCOMVel(double vdrift[3]){
261 +
262 +  double mtot, mtot_local;
263 +  double vdrift_local[3];
264 +  int vd, n_atoms;
265 +  Atom** atoms;
266 +
267 +  // We are very careless here with the distinction between n_atoms and n_local
268 +  // We should really fix this before someone pokes an eye out.
269 +
270 +  n_atoms = entry_plug->n_atoms;  
271 +  atoms   = entry_plug->atoms;
272 +
273 +  mtot_local = 0.0;
274 +  vdrift_local[0] = 0.0;
275 +  vdrift_local[1] = 0.0;
276 +  vdrift_local[2] = 0.0;
277 +  
278 +  for(vd = 0; vd < n_atoms; vd++){
279 +    
280 +    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
281 +    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
282 +    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
283 +    
284 +    mtot_local += atoms[vd]->getMass();
285 +  }
286 +
287 + #ifdef IS_MPI
288 +  MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM);
289 +  MPI::COMM_WORLD.Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM);
290 + #else
291 +  mtot = mtot_local;
292 +  for(vd = 0; vd < 3; vd++) {
293 +    vdrift[vd] = vdrift_local[vd];
294 +  }
295 + #endif
296 +    
297 +  for (vd = 0; vd < 3; vd++) {
298 +    vdrift[vd] = vdrift[vd] / mtot;
299 +  }
300 +  
301 + }
302 +

Diff Legend

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