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 252 by chuckv, Tue Jan 28 22:16:55 2003 UTC

# Line 1 | Line 1
1   #include <cmath>
2 + #include <iostream>
3  
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 +  gaussStream = new gaussianSPRNG( baseSeed );
20 + }
21 +
22 + Thermo::~Thermo(){
23 +  delete gaussStream;
24 + }
25 +
26   double Thermo::getKinetic(){
27  
28    const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
# Line 17 | Line 34 | double Thermo::getKinetic(){
34    DirectionalAtom *dAtom;
35  
36    int n_atoms;
37 +  double kinetic_global;
38    Atom** atoms;
39 +
40    
41    n_atoms = entry_plug->n_atoms;
42    atoms = entry_plug->atoms;
43  
44    kinetic = 0.0;
45 +  kinetic_global = 0.0;
46    for( kl=0; kl < n_atoms; kl++ ){
47  
48      vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
# Line 44 | Line 64 | double Thermo::getKinetic(){
64          + (jz2 / dAtom->getIzz());
65      }
66    }
67 <  
67 > #ifdef IS_MPI
68 >  MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
69 >  kinetic = kinetic_global;
70 > #endif //is_mpi
71 >
72    kinetic = kinetic * 0.5 / e_convert;
73  
74    return kinetic;
# Line 53 | Line 77 | double Thermo::getPotential(){
77   double Thermo::getPotential(){
78    
79    double potential;
80 +  double potential_global;
81    int el, nSRI;
82    SRI** sris;
83  
# Line 60 | Line 85 | double Thermo::getPotential(){
85    nSRI = entry_plug->n_SRI;
86  
87    potential = 0.0;
88 +  potential_global = 0.0;
89 +  potential += entry_plug->lrPot;
90  
64 –  potential += entry_plug->longRange->get_potential();;
65 –
91    // std::cerr << "long range potential: " << potential << "\n";
67 –
92    for( el=0; el<nSRI; el++ ){
93      
94      potential += sris[el]->get_potential();
95    }
96  
97 +  // Get total potential for entire system from MPI.
98 + #ifdef IS_MPI
99 +  MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
100 +  potential = potential_global;
101 + #endif // is_mpi
102 +
103    return potential;
104   }
105  
# Line 83 | Line 113 | double Thermo::getTemperature(){
113  
114   double Thermo::getTemperature(){
115  
116 <  const double kb = 1.88E-3; // boltzman's constant in kcal/(mol K)
116 >  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
117    double temperature;
118    
119    int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
# Line 95 | Line 125 | double Thermo::getPressure(){
125  
126   double Thermo::getPressure(){
127  
128 <  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
129 <  const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
130 <  const double conv_A_m = 1.0E-10; //convert A -> m
128 > //  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
129 > // const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
130 > //  const double conv_A_m = 1.0E-10; //convert A -> m
131  
132    return 0.0;
133   }
# Line 134 | Line 164 | void Thermo::velocitize() {
164    ndf = ndfRaw - n_constraints - 3;
165    kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
166    
167 +  printf("Entered Velocitize\n");
168    for(vr = 0; vr < n_atoms; vr++){
169      
170      // uses equipartition theory to solve for vbar in angstrom/fs
# Line 145 | Line 176 | void Thermo::velocitize() {
176      
177      // picks random velocities from a gaussian distribution
178      // centered on vbar
179 <    
179 > #ifndef USE_SPRNG
180 >    /* If we are using mpi, we need to use the SPRNG random
181 >       generator. The non drand48 generator will just repeat
182 >       the same numbers for every node creating a non-gaussian
183 >       distribution for the simulation. drand48 is fine for the
184 >       single processor version of the code, but SPRNG should
185 >       still be preferred for consistency.
186 >    */
187 >
188 > #ifdef IS_MPI
189 > #error "SPRNG random number generator must be used for MPI"
190 > #else
191 > #warning "Using drand48 for random number generation"
192 > #endif  // is_mpi
193 >
194      x = drand48();
195      y = drand48();
196      vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
# Line 157 | Line 202 | void Thermo::velocitize() {
202      x = drand48();
203      y = drand48();
204      vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
205 <    
205 >    printf("Setting new velocities vx: %f\n",vx);
206 > #endif // use_spring
207 >
208 > #ifdef USE_SPRNG
209 >    vx = vbar * gaussStream->getGaussian();
210 >    vy = vbar * gaussStream->getGaussian();
211 >    vz = vbar * gaussStream->getGaussian();
212 > #endif // use_spring
213 >
214      atoms[vr]->set_vx( vx );
215      atoms[vr]->set_vy( vy );
216      atoms[vr]->set_vz( vz );
# Line 205 | Line 258 | void Thermo::velocitize() {
258        if( atoms[i]->isDirectional() ){
259          
260          dAtom = (DirectionalAtom *)atoms[i];
261 +
262 + #ifndef USE_SPRNG
263 +
264 + #ifdef IS_MPI
265 + #error "SPRNG random number generator must be used for MPI"
266 + #else  // is_mpi
267 + #warning "Using drand48 for random number generation"
268 + #endif   // is_MPI
269          
270          vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
271          x = drand48();
# Line 220 | Line 281 | void Thermo::velocitize() {
281          x = drand48();
282          y = drand48();
283          jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
284 +
285 + #else //use_sprng
286 +
287 +        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
288 +        jx = vbar * gaussStream->getGaussian();
289 +
290 +        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
291 +        jy = vbar * gaussStream->getGaussian();
292 +
293 +        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
294 +        jz = vbar * gaussStream->getGaussian();
295 + #endif //use_sprng
296          
297          dAtom->setJx( jx );
298          dAtom->setJy( jy );

Diff Legend

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