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

Comparing trunk/mdtools/md_code/Thermo.cpp (file contents):
Revision 117 by mmeineke, Tue Sep 24 22:10:55 2002 UTC vs.
Revision 221 by chuckv, Thu Jan 2 20:14:08 2003 UTC

# Line 1 | Line 1
1   #include <cmath>
2 + #include <mpi++.h>
3  
4   #include "Thermo.hpp"
5   #include "SRI.hpp"
# Line 17 | Line 18 | double Thermo::getKinetic(){
18    DirectionalAtom *dAtom;
19  
20    int n_atoms;
21 +  double kinetic_global;
22    Atom** atoms;
23 +
24    
25    n_atoms = entry_plug->n_atoms;
26    atoms = entry_plug->atoms;
27  
28    kinetic = 0.0;
29 +  kinetic_global = 0.0;
30    for( kl=0; kl < n_atoms; kl++ ){
31  
32      vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
# Line 44 | Line 48 | double Thermo::getKinetic(){
48          + (jz2 / dAtom->getIzz());
49      }
50    }
51 <  
51 > #ifdef IS_MPI
52 >  MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
53 >  kinetic = kinetic_global;
54 > #endif
55 >
56    kinetic = kinetic * 0.5 / e_convert;
57  
58    return kinetic;
# Line 53 | Line 61 | double Thermo::getPotential(){
61   double Thermo::getPotential(){
62    
63    double potential;
64 +  double potential_global;
65    int el, nSRI;
66    SRI** sris;
67  
# Line 60 | Line 69 | double Thermo::getPotential(){
69    nSRI = entry_plug->n_SRI;
70  
71    potential = 0.0;
72 <
72 >  potential_global = 0.0;
73    potential += entry_plug->longRange->get_potential();;
74  
75    // std::cerr << "long range potential: " << potential << "\n";
67 –
76    for( el=0; el<nSRI; el++ ){
77      
78      potential += sris[el]->get_potential();
79    }
80  
81 +  // Get total potential for entire system from MPI.
82 + #ifdef IS_MPI
83 +  MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
84 +  potential = potential_global;
85 + #endif
86 +
87    return potential;
88   }
89  
# Line 145 | Line 159 | void Thermo::velocitize() {
159      
160      // picks random velocities from a gaussian distribution
161      // centered on vbar
162 <    
162 > #ifndef USE_SPRNG
163 >    /* If we are using mpi, we need to use the SPRNG random
164 >       generator. The non drand48 generator will just repeat
165 >       the same numbers for every node creating a non-gaussian
166 >       distribution for the simulation. drand48 is fine for the
167 >       single processor version of the code, but SPRNG should
168 >       still be preferred for consistency.
169 >    */
170 > #ifdef IS_MPI
171 > #error "SPRNG random number generator must be used for MPI"
172 > #else
173 > #warning "Using drand48 for random number generation"
174 > #endif  
175      x = drand48();
176      y = drand48();
177      vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
# Line 157 | Line 183 | void Thermo::velocitize() {
183      x = drand48();
184      y = drand48();
185      vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
186 <    
186 > #endif
187 >
188 > #ifdef USE_SPRNG
189 >    vx = vbar * entry_plug->gaussStream->getGaussian();
190 >    vy = vbar * entry_plug->gaussStream->getGaussian();
191 >    vz = vbar * entry_plug->gaussStream->getGaussian();
192 > #endif
193 >
194      atoms[vr]->set_vx( vx );
195      atoms[vr]->set_vy( vy );
196      atoms[vr]->set_vz( vz );
# Line 205 | Line 238 | void Thermo::velocitize() {
238        if( atoms[i]->isDirectional() ){
239          
240          dAtom = (DirectionalAtom *)atoms[i];
241 + #ifdef IS_MPI
242 + #error "SPRNG random number generator must be used for MPI"
243 + #else
244 + #warning "Using drand48 for random number generation"
245 + #endif  
246          
247          vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
248          x = drand48();
# Line 220 | Line 258 | void Thermo::velocitize() {
258          x = drand48();
259          y = drand48();
260          jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
261 + #endif
262 + #ifdef USE_SPRNG
263 +        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
264 +        jx = vbar * entry_plug->gaussStream->getGaussian();
265 +
266 +        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
267 +        jy = vbar * entry_plug->gaussStream->getGaussian();
268 +
269 +        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
270 +        jz = vbar * entry_plug->gaussStream->getGaussian();
271 + #endif
272          
273          dAtom->setJx( jx );
274          dAtom->setJy( jy );

Diff Legend

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