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 218 by chuckv, Sun Dec 29 19:11:05 2002 UTC vs.
Revision 253 by chuckv, Thu Jan 30 15:20:21 2003 UTC

# Line 1 | Line 1
1   #include <cmath>
2 + #include <iostream>
3 + using namespace std;
4 +
5 + #ifdef IS_MPI
6   #include <mpi++.h>
7 + #endif //is_mpi
8  
9   #include "Thermo.hpp"
10   #include "SRI.hpp"
11   #include "LRI.hpp"
12   #include "Integrator.hpp"
13  
14 + #define BASE_SEED 123456789
15  
16 + Thermo::Thermo( SimInfo* the_entry_plug ) {
17 +  entry_plug = the_entry_plug;
18 +  int baseSeed = BASE_SEED;
19 +  
20 +  cerr << "creating thermo stream\n";
21 +  gaussStream = new gaussianSPRNG( baseSeed );
22 +  cerr << "created thermo stream\n";
23 + }
24 +
25 + Thermo::~Thermo(){
26 +  delete gaussStream;
27 + }
28 +
29   double Thermo::getKinetic(){
30  
31    const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
# Line 51 | Line 70 | double Thermo::getKinetic(){
70   #ifdef IS_MPI
71    MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
72    kinetic = kinetic_global;
73 < #endif
73 > #endif //is_mpi
74  
75    kinetic = kinetic * 0.5 / e_convert;
76  
# Line 70 | Line 89 | double Thermo::getPotential(){
89  
90    potential = 0.0;
91    potential_global = 0.0;
92 <  potential += entry_plug->longRange->get_potential();;
92 >  potential += entry_plug->lrPot;
93  
94    // std::cerr << "long range potential: " << potential << "\n";
95    for( el=0; el<nSRI; el++ ){
# Line 82 | Line 101 | double Thermo::getPotential(){
101   #ifdef IS_MPI
102    MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
103    potential = potential_global;
104 < #endif
104 > #endif // is_mpi
105  
106    return potential;
107   }
# Line 159 | Line 178 | void Thermo::velocitize() {
178      
179      // picks random velocities from a gaussian distribution
180      // centered on vbar
181 <    
181 > #ifndef USE_SPRNG
182 >    /* If we are using mpi, we need to use the SPRNG random
183 >       generator. The non drand48 generator will just repeat
184 >       the same numbers for every node creating a non-gaussian
185 >       distribution for the simulation. drand48 is fine for the
186 >       single processor version of the code, but SPRNG should
187 >       still be preferred for consistency.
188 >    */
189 >
190 > #ifdef IS_MPI
191 > #error "SPRNG random number generator must be used for MPI"
192 > #else
193 > #warning "Using drand48 for random number generation"
194 > #endif  // is_mpi
195 >
196      x = drand48();
197      y = drand48();
198      vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
# Line 171 | Line 204 | void Thermo::velocitize() {
204      x = drand48();
205      y = drand48();
206      vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
207 <    
207 >
208 > #endif // use_spring
209 >
210 > #ifdef USE_SPRNG
211 >    vx = vbar * gaussStream->getGaussian();
212 >    vy = vbar * gaussStream->getGaussian();
213 >    vz = vbar * gaussStream->getGaussian();
214 > #endif // use_spring
215 >
216      atoms[vr]->set_vx( vx );
217      atoms[vr]->set_vy( vy );
218      atoms[vr]->set_vz( vz );
# Line 219 | Line 260 | void Thermo::velocitize() {
260        if( atoms[i]->isDirectional() ){
261          
262          dAtom = (DirectionalAtom *)atoms[i];
263 +
264 + #ifndef USE_SPRNG
265 +
266 + #ifdef IS_MPI
267 + #error "SPRNG random number generator must be used for MPI"
268 + #else  // is_mpi
269 + #warning "Using drand48 for random number generation"
270 + #endif   // is_MPI
271          
272          vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
273          x = drand48();
# Line 234 | Line 283 | void Thermo::velocitize() {
283          x = drand48();
284          y = drand48();
285          jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
286 +
287 + #else //use_sprng
288 +
289 +        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
290 +        jx = vbar * gaussStream->getGaussian();
291 +
292 +        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
293 +        jy = vbar * gaussStream->getGaussian();
294 +
295 +        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
296 +        jz = vbar * gaussStream->getGaussian();
297 + #endif //use_sprng
298          
299          dAtom->setJx( jx );
300          dAtom->setJy( jy );

Diff Legend

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