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 11 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
Revision 254 by chuckv, Thu Jan 30 20:03:37 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 + #include <mpi++.h>
8 + #endif //is_mpi
9 +
10   #include "Thermo.hpp"
11   #include "SRI.hpp"
12   #include "LRI.hpp"
13   #include "Integrator.hpp"
14  
15 + #define BASE_SEED 123456789
16  
17 + Thermo::Thermo( SimInfo* the_entry_plug ) {
18 +  entry_plug = the_entry_plug;
19 +  int baseSeed = BASE_SEED;
20 +  
21 +  cerr << "creating thermo stream\n";
22 +  gaussStream = new gaussianSPRNG( baseSeed );
23 +  cerr << "created thermo stream\n";
24 + }
25 +
26 + Thermo::~Thermo(){
27 +  delete gaussStream;
28 + }
29 +
30   double Thermo::getKinetic(){
31  
32    const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
# Line 17 | Line 38 | double Thermo::getKinetic(){
38    DirectionalAtom *dAtom;
39  
40    int n_atoms;
41 +  double kinetic_global;
42    Atom** atoms;
43 +
44    
45    n_atoms = entry_plug->n_atoms;
46    atoms = entry_plug->atoms;
47  
48    kinetic = 0.0;
49 +  kinetic_global = 0.0;
50    for( kl=0; kl < n_atoms; kl++ ){
51  
52      vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
# Line 44 | Line 68 | double Thermo::getKinetic(){
68          + (jz2 / dAtom->getIzz());
69      }
70    }
71 <  
71 > #ifdef IS_MPI
72 >  MPI::COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
73 >  kinetic = kinetic_global;
74 > #endif //is_mpi
75 >
76    kinetic = kinetic * 0.5 / e_convert;
77  
78    return kinetic;
# Line 53 | Line 81 | double Thermo::getPotential(){
81   double Thermo::getPotential(){
82    
83    double potential;
84 +  double potential_global;
85    int el, nSRI;
86    SRI** sris;
87  
# Line 60 | Line 89 | double Thermo::getPotential(){
89    nSRI = entry_plug->n_SRI;
90  
91    potential = 0.0;
92 +  potential_global = 0.0;
93 +  potential += entry_plug->lrPot;
94  
64 –  potential += entry_plug->longRange->get_potential();;
65 –
95    // std::cerr << "long range potential: " << potential << "\n";
67 –
96    for( el=0; el<nSRI; el++ ){
97      
98      potential += 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 + #endif // is_mpi
106 +
107    return potential;
108   }
109  
# Line 83 | Line 117 | double Thermo::getTemperature(){
117  
118   double Thermo::getTemperature(){
119  
120 <  const double kb = 1.88E-3; // boltzman's constant in kcal/(mol K)
120 >  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
121    double temperature;
122    
123    int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
# Line 95 | Line 129 | double Thermo::getPressure(){
129  
130   double Thermo::getPressure(){
131  
132 <  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
133 <  const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
134 <  const double conv_A_m = 1.0E-10; //convert A -> m
132 > //  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
133 > // const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
134 > //  const double conv_A_m = 1.0E-10; //convert A -> m
135  
136    return 0.0;
137   }
# Line 145 | Line 179 | void Thermo::velocitize() {
179      
180      // picks random velocities from a gaussian distribution
181      // centered on vbar
182 <    
182 > #ifndef USE_SPRNG
183 >    /* If we are using mpi, we need to use the SPRNG random
184 >       generator. The non drand48 generator will just repeat
185 >       the same numbers for every node creating a non-gaussian
186 >       distribution for the simulation. drand48 is fine for the
187 >       single processor version of the code, but SPRNG should
188 >       still be preferred for consistency.
189 >    */
190 >
191 > #ifdef IS_MPI
192 > #error "SPRNG random number generator must be used for MPI"
193 > #else
194 >    // warning "Using drand48 for random number generation"
195 > #endif  // is_mpi
196 >
197      x = drand48();
198      y = drand48();
199      vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
# Line 157 | Line 205 | void Thermo::velocitize() {
205      x = drand48();
206      y = drand48();
207      vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
208 <    
208 >
209 > #endif // use_spring
210 >
211 > #ifdef USE_SPRNG
212 >    vx = vbar * gaussStream->getGaussian();
213 >    vy = vbar * gaussStream->getGaussian();
214 >    vz = vbar * gaussStream->getGaussian();
215 > #endif // use_spring
216 >
217      atoms[vr]->set_vx( vx );
218      atoms[vr]->set_vy( vy );
219      atoms[vr]->set_vz( vz );
# Line 205 | Line 261 | void Thermo::velocitize() {
261        if( atoms[i]->isDirectional() ){
262          
263          dAtom = (DirectionalAtom *)atoms[i];
264 +
265 + #ifndef USE_SPRNG
266 +
267 + #ifdef IS_MPI
268 + #error "SPRNG random number generator must be used for MPI"
269 + #else  // is_mpi
270 +        //warning "Using drand48 for random number generation"
271 + #endif   // is_MPI
272          
273          vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
274          x = drand48();
# Line 220 | Line 284 | void Thermo::velocitize() {
284          x = drand48();
285          y = drand48();
286          jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
287 +
288 + #else //use_sprng
289 +
290 +        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
291 +        jx = vbar * gaussStream->getGaussian();
292 +
293 +        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
294 +        jy = vbar * gaussStream->getGaussian();
295 +
296 +        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
297 +        jz = vbar * gaussStream->getGaussian();
298 + #endif //use_sprng
299          
300          dAtom->setJx( jx );
301          dAtom->setJy( jy );

Diff Legend

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