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 447 by mmeineke, Thu Apr 3 20:21:54 2003 UTC vs.
Revision 454 by gezelter, Fri Apr 4 01:57:11 2003 UTC

# Line 125 | Line 125 | double Thermo::getTemperature(){
125    return total;
126   }
127  
128 < double Thermo::getTemperature(){
129 <
130 <  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
131 <  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
# Line 141 | Line 138 | double Thermo::getTemperature(){
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 180 | 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;
185 <
186 <  // Degrees of freedom that can contain kinetic energy
187 <  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
188 <    - entry_plug->n_constraints;
204 >  kebar = kb * temperature * (double)this->getNDF() /
205 >    ( 2.0 * (double)this->getNDFraw() );
206    
190 #ifdef IS_MPI
191  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
192  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
193 #else
194  ndfRaw = ndfRaw_local;
195  ndf = ndf_local;
196 #endif
197  ndf = ndf - 3;
198
199  kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
200  
207    for(vr = 0; vr < n_atoms; vr++){
208      
209      // uses equipartition theory to solve for vbar in angstrom/fs
# Line 253 | 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          

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines