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

Comparing trunk/OOPSE/libmdtools/Thermo.cpp (file contents):
Revision 445 by gezelter, Thu Apr 3 19:58:24 2003 UTC vs.
Revision 454 by gezelter, Fri Apr 4 01:57:11 2003 UTC

# Line 4 | Line 4 | using namespace std;
4  
5   #ifdef IS_MPI
6   #include <mpi.h>
7 #include <mpi++.h>
7   #endif //is_mpi
8  
9   #include "Thermo.hpp"
# Line 73 | Line 72 | double Thermo::getKinetic(){
72      }
73    }
74   #ifdef IS_MPI
75 <  MPI::COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
75 >  MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
76 >                MPI_SUM, MPI_COMM_WORLD);
77    kinetic = kinetic_global;
78   #endif //is_mpi
79  
# Line 100 | Line 100 | double Thermo::getPotential(){
100      potential_local += molecules[el].getPotential();
101    }
102  
103 #ifdef IS_MPI
104  /*
105  std::cerr << "node " << worldRank << ": before LONG RANGE pot = " << entry_plug->lrPot
106            << "; pot_local = " << potential_local
107            << "; pot = " << potential << "\n";
108  */
109 #endif
110
103    // Get total potential for entire system from MPI.
104   #ifdef IS_MPI
105 <  MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM);
105 >  MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
106 >                MPI_SUM, MPI_COMM_WORLD);
107   #else
108    potential = potential_local;
109   #endif // is_mpi
# Line 132 | Line 125 | double Thermo::getTemperature(){
125    return total;
126   }
127  
128 < double Thermo::getTemperature(){
136 <
137 <  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
138 <  double temperature;
128 > int Thermo::getNDF(){
129    int ndf_local, ndf;
130    
131    ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
132      - entry_plug->n_constraints;
133  
134   #ifdef IS_MPI
135 <  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
135 >  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
136   #else
137    ndf = ndf_local;
138   #endif
139  
140    ndf = ndf - 3;
141 +
142 +  return ndf;
143 + }
144 +
145 + int Thermo::getNDFraw() {
146 +  int ndfRaw_local, ndfRaw;
147 +
148 +  // Raw degrees of freedom that we have to set
149 +  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
150    
151 <  temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
151 > #ifdef IS_MPI
152 >  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
153 > #else
154 >  ndfRaw = ndfRaw_local;
155 > #endif
156 >
157 >  return ndfRaw;
158 > }
159 >
160 >
161 > double Thermo::getTemperature(){
162 >
163 >  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
164 >  double temperature;
165 >  
166 >  temperature = ( 2.0 * this->getKinetic() ) / ( (double)this->getNDF() * kb );
167    return temperature;
168   }
169  
# Line 187 | Line 201 | void Thermo::velocitize() {
201    n_oriented    = entry_plug->n_oriented;
202    n_constraints = entry_plug->n_constraints;
203    
204 <  // Raw degrees of freedom that we have to set
205 <  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
192 <
193 <  // Degrees of freedom that can contain kinetic energy
194 <  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
195 <    - entry_plug->n_constraints;
204 >  kebar = kb * temperature * (double)this->getNDF() /
205 >    ( 2.0 * (double)this->getNDFraw() );
206    
197 #ifdef IS_MPI
198  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
199  MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM);
200 #else
201  ndfRaw = ndfRaw_local;
202  ndf = ndf_local;
203 #endif
204  ndf = ndf - 3;
205
206  kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
207  
207    for(vr = 0; vr < n_atoms; vr++){
208      
209      // uses equipartition theory to solve for vbar in angstrom/fs
# Line 260 | Line 259 | void Thermo::velocitize() {
259  
260          vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
261          jy = vbar * gaussStream->getGaussian();
262 <
262 >        
263          vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
264          jz = vbar * gaussStream->getGaussian();
265          
# Line 300 | Line 299 | void Thermo::getCOMVel(double vdrift[3]){
299    }
300  
301   #ifdef IS_MPI
302 <  MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM);
303 <  MPI::COMM_WORLD.Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM);
302 >  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
303 >  MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
304   #else
305    mtot = mtot_local;
306    for(vd = 0; vd < 3; vd++) {

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines