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 223 by chuckv, Fri Jan 3 22:04:50 2003 UTC

# Line 1 | Line 1
1   #include <cmath>
2 + #include <mpi++.h>
3  
4   #include "Thermo.hpp"
5   #include "SRI.hpp"
6   #include "LRI.hpp"
7   #include "Integrator.hpp"
8  
9 + #define BASE_SEED 123456789
10  
11 + Thermo::Thermo( SimInfo* the_entry_plug ) {
12 +  entry_plug = the_entry_plug;
13 +  baseSeed = BASE_SEED;
14 +  gaussStream = new gaussianSPRNG( baseSeed );
15 + }
16 +
17 + Thermo::~Thermo(){
18 +  delete gaussStream;
19 + }
20 +
21   double Thermo::getKinetic(){
22  
23    const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
# Line 17 | Line 29 | double Thermo::getKinetic(){
29    DirectionalAtom *dAtom;
30  
31    int n_atoms;
32 +  double kinetic_global;
33    Atom** atoms;
34 +
35    
36    n_atoms = entry_plug->n_atoms;
37    atoms = entry_plug->atoms;
38  
39    kinetic = 0.0;
40 +  kinetic_global = 0.0;
41    for( kl=0; kl < n_atoms; kl++ ){
42  
43      vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
# Line 44 | Line 59 | double Thermo::getKinetic(){
59          + (jz2 / dAtom->getIzz());
60      }
61    }
62 <  
62 > #ifdef IS_MPI
63 >  MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
64 >  kinetic = kinetic_global;
65 > #endif
66 >
67    kinetic = kinetic * 0.5 / e_convert;
68  
69    return kinetic;
# Line 53 | Line 72 | double Thermo::getPotential(){
72   double Thermo::getPotential(){
73    
74    double potential;
75 +  double potential_global;
76    int el, nSRI;
77    SRI** sris;
78  
# Line 60 | Line 80 | double Thermo::getPotential(){
80    nSRI = entry_plug->n_SRI;
81  
82    potential = 0.0;
83 <
83 >  potential_global = 0.0;
84    potential += entry_plug->longRange->get_potential();;
85  
86    // std::cerr << "long range potential: " << potential << "\n";
67 –
87    for( el=0; el<nSRI; el++ ){
88      
89      potential += sris[el]->get_potential();
90    }
91  
92 +  // Get total potential for entire system from MPI.
93 + #ifdef IS_MPI
94 +  MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
95 +  potential = potential_global;
96 + #endif
97 +
98    return potential;
99   }
100  
# Line 145 | Line 170 | void Thermo::velocitize() {
170      
171      // picks random velocities from a gaussian distribution
172      // centered on vbar
173 <    
173 > #ifndef USE_SPRNG
174 >    /* If we are using mpi, we need to use the SPRNG random
175 >       generator. The non drand48 generator will just repeat
176 >       the same numbers for every node creating a non-gaussian
177 >       distribution for the simulation. drand48 is fine for the
178 >       single processor version of the code, but SPRNG should
179 >       still be preferred for consistency.
180 >    */
181 > #ifdef IS_MPI
182 > #error "SPRNG random number generator must be used for MPI"
183 > #else
184 > #warning "Using drand48 for random number generation"
185 > #endif  
186      x = drand48();
187      y = drand48();
188      vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
# Line 157 | Line 194 | void Thermo::velocitize() {
194      x = drand48();
195      y = drand48();
196      vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
197 <    
197 > #endif
198 >
199 > #ifdef USE_SPRNG
200 >    vx = vbar * gaussStream->getGaussian();
201 >    vy = vbar * gaussStream->getGaussian();
202 >    vz = vbar * gaussStream->getGaussian();
203 > #endif
204 >
205      atoms[vr]->set_vx( vx );
206      atoms[vr]->set_vy( vy );
207      atoms[vr]->set_vz( vz );
# Line 205 | Line 249 | void Thermo::velocitize() {
249        if( atoms[i]->isDirectional() ){
250          
251          dAtom = (DirectionalAtom *)atoms[i];
252 + #ifdef IS_MPI
253 + #error "SPRNG random number generator must be used for MPI"
254 + #else
255 + #warning "Using drand48 for random number generation"
256 + #endif  
257          
258          vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
259          x = drand48();
# Line 220 | Line 269 | void Thermo::velocitize() {
269          x = drand48();
270          y = drand48();
271          jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
272 + #endif
273 + #ifdef USE_SPRNG
274 +        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
275 +        jx = vbar * gaussStream->getGaussian();
276 +
277 +        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
278 +        jy = vbar * gaussStream->getGaussian();
279 +
280 +        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
281 +        jz = vbar * gaussStream->getGaussian();
282 + #endif
283          
284          dAtom->setJx( jx );
285          dAtom->setJy( jy );

Diff Legend

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