ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/Thermo.cpp
(Generate patch)

Comparing trunk/src/brains/Thermo.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1292 by chuckv, Fri Sep 12 20:51:22 2008 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 | namespace oopse {
53  
54   namespace oopse {
55  
56 < double Thermo::getKinetic() {
56 >  RealType Thermo::getKinetic() {
57      SimInfo::MoleculeIterator miter;
58      std::vector<StuntDouble*>::iterator iiter;
59      Molecule* mol;
# Line 64 | Line 64 | double Thermo::getKinetic() {
64      int i;
65      int j;
66      int k;
67 <    double kinetic = 0.0;
68 <    double kinetic_global = 0.0;
67 >    RealType mass;
68 >    RealType kinetic = 0.0;
69 >    RealType kinetic_global = 0.0;
70      
71      for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
72 <        for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
73 <               integrableObject = mol->nextIntegrableObject(iiter)) {
72 >      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
73 >           integrableObject = mol->nextIntegrableObject(iiter)) {
74 >        
75 >        mass = integrableObject->getMass();
76 >        vel = integrableObject->getVel();
77 >        
78 >        kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
79 >        
80 >        if (integrableObject->isDirectional()) {
81 >          angMom = integrableObject->getJ();
82 >          I = integrableObject->getI();
83  
84 <            double mass = integrableObject->getMass();
85 <            Vector3d vel = integrableObject->getVel();
86 <
87 <            kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
88 <
89 <            if (integrableObject->isDirectional()) {
90 <                angMom = integrableObject->getJ();
91 <                I = integrableObject->getI();
92 <
93 <                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 <            }
84 >          if (integrableObject->isLinear()) {
85 >            i = integrableObject->linearAxis();
86 >            j = (i + 1) % 3;
87 >            k = (i + 2) % 3;
88 >            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
89 >          } else {                        
90 >            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
91 >              + angMom[2]*angMom[2]/I(2, 2);
92 >          }
93 >        }
94              
95 <        }
95 >      }
96      }
97      
98   #ifdef IS_MPI
99  
100 <    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_DOUBLE, MPI_SUM,
100 >    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
101                    MPI_COMM_WORLD);
102      kinetic = kinetic_global;
103  
# Line 105 | Line 106 | double Thermo::getKinetic() {
106      kinetic = kinetic * 0.5 / OOPSEConstant::energyConvert;
107  
108      return kinetic;
109 < }
109 >  }
110  
111 < double Thermo::getPotential() {
112 <    double potential = 0.0;
111 >  RealType Thermo::getPotential() {
112 >    RealType potential = 0.0;
113      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
114 <    double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] +
114 <                                             curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
114 >    RealType shortRangePot_local =  curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
115  
116      // Get total potential for entire system from MPI.
117  
118   #ifdef IS_MPI
119  
120 <    MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM,
120 >    MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM,
121                    MPI_COMM_WORLD);
122 +    potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
123  
124   #else
125  
126 <    potential = potential_local;
126 >    potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
127  
128   #endif // is_mpi
129  
130      return potential;
131 < }
131 >  }
132  
133 < double Thermo::getTotalE() {
134 <    double total;
133 >  RealType Thermo::getTotalE() {
134 >    RealType total;
135  
136      total = this->getKinetic() + this->getPotential();
137      return total;
138 < }
138 >  }
139  
140 < double Thermo::getTemperature() {
140 >  RealType Thermo::getTemperature() {
141      
142 <    double temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb );
142 >    RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb );
143      return temperature;
144 < }
144 >  }
145  
146 < double Thermo::getVolume() {
146 >  RealType Thermo::getVolume() {
147      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
148      return curSnapshot->getVolume();
149 < }
149 >  }
150  
151 < double Thermo::getPressure() {
151 >  RealType Thermo::getPressure() {
152  
153      // Relies on the calculation of the full molecular pressure tensor
154  
155  
156      Mat3x3d tensor;
157 <    double pressure;
157 >    RealType pressure;
158  
159      tensor = getPressureTensor();
160  
161      pressure = OOPSEConstant::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0;
162  
163      return pressure;
164 < }
164 >  }
165  
166 < Mat3x3d Thermo::getPressureTensor() {
166 >  RealType Thermo::getPressure(int direction) {
167 >
168 >    // Relies on the calculation of the full molecular pressure tensor
169 >
170 >          
171 >    Mat3x3d tensor;
172 >    RealType pressure;
173 >
174 >    tensor = getPressureTensor();
175 >
176 >    pressure = OOPSEConstant::pressureConvert * tensor(direction, direction);
177 >
178 >    return pressure;
179 >  }
180 >
181 >  Mat3x3d Thermo::getPressureTensor() {
182      // returns pressure tensor in units amu*fs^-2*Ang^-1
183      // routine derived via viral theorem description in:
184      // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
# Line 175 | Line 191 | Mat3x3d Thermo::getPressureTensor() {
191      Molecule* mol;
192      StuntDouble* integrableObject;    
193      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
194 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
195 <               integrableObject = mol->nextIntegrableObject(j)) {
194 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
195 >           integrableObject = mol->nextIntegrableObject(j)) {
196  
197 <            double mass = integrableObject->getMass();
198 <            Vector3d vcom = integrableObject->getVel();
199 <            p_local += mass * outProduct(vcom, vcom);        
200 <        }
197 >        RealType mass = integrableObject->getMass();
198 >        Vector3d vcom = integrableObject->getVel();
199 >        p_local += mass * outProduct(vcom, vcom);        
200 >      }
201      }
202      
203   #ifdef IS_MPI
204 <    MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
204 >    MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
205   #else
206      p_global = p_local;
207   #endif // is_mpi
208  
209 <    double volume = this->getVolume();
209 >    RealType volume = this->getVolume();
210      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
211      Mat3x3d tau = curSnapshot->statData.getTau();
212  
213      pressureTensor =  (p_global + OOPSEConstant::energyConvert* tau)/volume;
214 <
214 >    
215      return pressureTensor;
216 < }
216 >  }
217  
218 < void Thermo::saveStat(){
218 >
219 >  void Thermo::saveStat(){
220      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
221      Stats& stat = currSnapshot->statData;
222      
# Line 210 | Line 227 | void Thermo::saveStat(){
227      stat[Stats::PRESSURE] = getPressure();
228      stat[Stats::VOLUME] = getVolume();      
229  
230 +    Mat3x3d tensor =getPressureTensor();
231 +    stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0);      
232 +    stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1);      
233 +    stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2);      
234 +    stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0);      
235 +    stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1);      
236 +    stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2);      
237 +    stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0);      
238 +    stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
239 +    stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
240 +
241 +
242 +    Globals* simParams = info_->getSimParams();
243 +
244 +    if (simParams->haveTaggedAtomPair() &&
245 +        simParams->havePrintTaggedPairDistance()) {
246 +      if ( simParams->getPrintTaggedPairDistance()) {
247 +        
248 +        std::pair<int, int> tap = simParams->getTaggedAtomPair();
249 +        Vector3d pos1, pos2, rab;
250 +
251 + #ifdef IS_MPI        
252 +
253 +        int mol1 = info_->getGlobalMolMembership(tap.first);
254 +        int mol2 = info_->getGlobalMolMembership(tap.second);
255 +        
256 +        int proc1 = info_->getMolToProc(mol1);
257 +        int proc2 = info_->getMolToProc(mol2);
258 +
259 +        RealType data[3];
260 +        if (proc1 == worldRank) {
261 +          StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
262 +          pos1 = sd1->getPos();
263 +          data[0] = pos1.x();
264 +          data[1] = pos1.y();
265 +          data[2] = pos1.z();          
266 +          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
267 +        } else {
268 +          MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
269 +          pos1 = Vector3d(data);
270 +        }
271 +
272 +
273 +        if (proc2 == worldRank) {
274 +          StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second);
275 +          pos2 = sd2->getPos();
276 +          data[0] = pos2.x();
277 +          data[1] = pos2.y();
278 +          data[2] = pos2.z();          
279 +          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
280 +        } else {
281 +          MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
282 +          pos2 = Vector3d(data);
283 +        }
284 + #else
285 +        StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first);
286 +        StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
287 +        pos1 = at1->getPos();
288 +        pos2 = at2->getPos();
289 + #endif        
290 +        rab = pos2 - pos1;
291 +        currSnapshot->wrapVector(rab);
292 +        stat[Stats::TAGGED_PAIR_DISTANCE] =  rab.length();
293 +      }
294 +    }
295 +      
296      /**@todo need refactorying*/
297      //Conserved Quantity is set by integrator and time is set by setTime
298      
299 < }
299 >  }
300  
301   } //end namespace oopse

Diff Legend

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