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

Comparing trunk/OOPSE-3.0/src/brains/Thermo.cpp (file contents):
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2235 by tim, Thu May 19 21:31:23 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 53 | Line 53 | double Thermo::getKinetic() {
53  
54   namespace oopse {
55  
56 < double Thermo::getKinetic() {
56 >  double Thermo::getKinetic() {
57      SimInfo::MoleculeIterator miter;
58      std::vector<StuntDouble*>::iterator iiter;
59      Molecule* mol;
# Line 68 | Line 68 | double Thermo::getKinetic() {
68      double kinetic_global = 0.0;
69      
70      for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
71 <        for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
72 <               integrableObject = mol->nextIntegrableObject(iiter)) {
71 >      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
72 >           integrableObject = mol->nextIntegrableObject(iiter)) {
73  
74 <            double mass = integrableObject->getMass();
75 <            Vector3d vel = integrableObject->getVel();
74 >        double mass = integrableObject->getMass();
75 >        Vector3d vel = integrableObject->getVel();
76  
77 <            kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
77 >        kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
78  
79 <            if (integrableObject->isDirectional()) {
80 <                angMom = integrableObject->getJ();
81 <                I = integrableObject->getI();
79 >        if (integrableObject->isDirectional()) {
80 >          angMom = integrableObject->getJ();
81 >          I = integrableObject->getI();
82  
83 <                if (integrableObject->isLinear()) {
84 <                    i = integrableObject->linearAxis();
85 <                    j = (i + 1) % 3;
86 <                    k = (i + 2) % 3;
87 <                    kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
88 <                } else {                        
89 <                    kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
90 <                                    + angMom[2]*angMom[2]/I(2, 2);
91 <                }
92 <            }
83 >          if (integrableObject->isLinear()) {
84 >            i = integrableObject->linearAxis();
85 >            j = (i + 1) % 3;
86 >            k = (i + 2) % 3;
87 >            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
88 >          } else {                        
89 >            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
90 >              + angMom[2]*angMom[2]/I(2, 2);
91 >          }
92 >        }
93              
94 <        }
94 >      }
95      }
96      
97   #ifdef IS_MPI
# Line 105 | Line 105 | double Thermo::getKinetic() {
105      kinetic = kinetic * 0.5 / OOPSEConstant::energyConvert;
106  
107      return kinetic;
108 < }
108 >  }
109  
110 < double Thermo::getPotential() {
110 >  double Thermo::getPotential() {
111      double potential = 0.0;
112      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
113      double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +
114 <                                             curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
114 >      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
115  
116      // Get total potential for entire system from MPI.
117  
# Line 127 | Line 127 | double Thermo::getPotential() {
127   #endif // is_mpi
128  
129      return potential;
130 < }
130 >  }
131  
132 < double Thermo::getTotalE() {
132 >  double Thermo::getTotalE() {
133      double total;
134  
135      total = this->getKinetic() + this->getPotential();
136      return total;
137 < }
137 >  }
138  
139 < double Thermo::getTemperature() {
139 >  double Thermo::getTemperature() {
140      
141      double temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb );
142      return temperature;
143 < }
143 >  }
144  
145 < double Thermo::getVolume() {
145 >  double Thermo::getVolume() {
146      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
147      return curSnapshot->getVolume();
148 < }
148 >  }
149  
150 < double Thermo::getPressure() {
150 >  double Thermo::getPressure() {
151  
152      // Relies on the calculation of the full molecular pressure tensor
153  
# Line 160 | Line 160 | double Thermo::getPressure() {
160      pressure = OOPSEConstant::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
161  
162      return pressure;
163 < }
163 >  }
164  
165 < Mat3x3d Thermo::getPressureTensor() {
165 >  double Thermo::getPressure(int direction) {
166 >
167 >    // Relies on the calculation of the full molecular pressure tensor
168 >
169 >          
170 >    Mat3x3d tensor;
171 >    double pressure;
172 >
173 >    tensor = getPressureTensor();
174 >
175 >    pressure = OOPSEConstant::pressureConvert * tensor(direction, direction);
176 >
177 >    return pressure;
178 >  }
179 >
180 >
181 >
182 >  Mat3x3d Thermo::getPressureTensor() {
183      // returns pressure tensor in units amu*fs^-2*Ang^-1
184      // routine derived via viral theorem description in:
185      // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
# Line 175 | Line 192 | Mat3x3d Thermo::getPressureTensor() {
192      Molecule* mol;
193      StuntDouble* integrableObject;    
194      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
195 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
196 <               integrableObject = mol->nextIntegrableObject(j)) {
195 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
196 >           integrableObject = mol->nextIntegrableObject(j)) {
197  
198 <            double mass = integrableObject->getMass();
199 <            Vector3d vcom = integrableObject->getVel();
200 <            p_local += mass * outProduct(vcom, vcom);        
201 <        }
198 >        double mass = integrableObject->getMass();
199 >        Vector3d vcom = integrableObject->getVel();
200 >        p_local += mass * outProduct(vcom, vcom);        
201 >      }
202      }
203      
204   #ifdef IS_MPI
# Line 197 | Line 214 | Mat3x3d Thermo::getPressureTensor() {
214      pressureTensor =  (p_global + OOPSEConstant::energyConvert* tau)/volume;
215  
216      return pressureTensor;
217 < }
217 >  }
218  
219 < void Thermo::saveStat(){
219 >  void Thermo::saveStat(){
220      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
221      Stats& stat = currSnapshot->statData;
222      
# Line 213 | Line 230 | void Thermo::saveStat(){
230      /**@todo need refactorying*/
231      //Conserved Quantity is set by integrator and time is set by setTime
232      
233 < }
233 >  }
234  
235   } //end namespace oopse

Diff Legend

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