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

Comparing trunk/src/brains/Thermo.cpp (file contents):
Revision 1665 by chuckv, Fri Sep 23 20:52:24 2011 UTC vs.
Revision 1666 by chuckv, Wed Dec 14 20:21:54 2011 UTC

# Line 32 | Line 32
32   * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33   * research, please cite the appropriate papers when you publish your
34   * work.  Good starting points are:
35 < *                                                                      
36 < * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 < * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
35 > *
36 > * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 > * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39 > * [4]  Vardeman & Gezelter, in progress (2009).
40   */
41 <
41 >
42   #include <math.h>
43   #include <iostream>
44  
# Line 57 | Line 57 | namespace OpenMD {
57      SimInfo::MoleculeIterator miter;
58      std::vector<StuntDouble*>::iterator iiter;
59      Molecule* mol;
60 <    StuntDouble* integrableObject;    
60 >    StuntDouble* integrableObject;
61      Vector3d vel;
62      Vector3d angMom;
63      Mat3x3d I;
# Line 67 | Line 67 | namespace OpenMD {
67      RealType mass;
68      RealType kinetic = 0.0;
69      RealType kinetic_global = 0.0;
70 <    
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)) {
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();
72 >      for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL;
73 >     integrableObject = mol->nextIntegrableObject(iiter)) {
74  
75 <          if (integrableObject->isLinear()) {
76 <            i = integrableObject->linearAxis();
77 <            j = (i + 1) % 3;
78 <            k = (i + 2) % 3;
79 <            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
80 <          } else {                        
81 <            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
82 <              + angMom[2]*angMom[2]/I(2, 2);
83 <          }
84 <        }
85 <            
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 >    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        }
96      }
97 <    
97 >
98   #ifdef IS_MPI
99  
100      MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
# Line 138 | Line 138 | namespace OpenMD {
138    }
139  
140    RealType Thermo::getTemperature() {
141 <    
141 >
142      RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb );
143      return temperature;
144    }
145  
146 <  RealType Thermo::getVolume() {
146 >  RealType Thermo::getVolume() {
147      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
148      return curSnapshot->getVolume();
149    }
# Line 167 | Line 167 | namespace OpenMD {
167  
168      // Relies on the calculation of the full molecular pressure tensor
169  
170 <          
170 >
171      Mat3x3d tensor;
172      RealType pressure;
173  
# Line 189 | Line 189 | namespace OpenMD {
189      SimInfo::MoleculeIterator i;
190      std::vector<StuntDouble*>::iterator j;
191      Molecule* mol;
192 <    StuntDouble* integrableObject;    
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 <        RealType mass = integrableObject->getMass();
198 <        Vector3d vcom = integrableObject->getVel();
199 <        p_local += mass * outProduct(vcom, vcom);        
197 >  RealType mass = integrableObject->getMass();
198 >  Vector3d vcom = integrableObject->getVel();
199 >  p_local += mass * outProduct(vcom, vcom);
200        }
201      }
202 <    
202 >
203   #ifdef IS_MPI
204      MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
205   #else
# Line 211 | Line 211 | namespace OpenMD {
211      Mat3x3d tau = curSnapshot->statData.getTau();
212  
213      pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
214 <    
214 >
215      return pressureTensor;
216    }
217  
# Line 219 | Line 219 | namespace OpenMD {
219    void Thermo::saveStat(){
220      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
221      Stats& stat = currSnapshot->statData;
222 <    
222 >
223      stat[Stats::KINETIC_ENERGY] = getKinetic();
224      stat[Stats::POTENTIAL_ENERGY] = getPotential();
225      stat[Stats::TOTAL_ENERGY] = stat[Stats::KINETIC_ENERGY]  + stat[Stats::POTENTIAL_ENERGY] ;
226      stat[Stats::TEMPERATURE] = getTemperature();
227      stat[Stats::PRESSURE] = getPressure();
228 <    stat[Stats::VOLUME] = getVolume();      
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);      
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      Vector3d GKappa_t = getThermalHelfand();
241      stat[Stats::THERMAL_HELFANDMOMENT_X] = GKappa_t.x();
242      stat[Stats::THERMAL_HELFANDMOMENT_Y] = GKappa_t.y();
# Line 244 | Line 244 | namespace OpenMD {
244  
245      Globals* simParams = info_->getSimParams();
246  
247 <    if (simParams->haveTaggedAtomPair() &&
247 >    if (simParams->haveTaggedAtomPair() &&
248          simParams->havePrintTaggedPairDistance()) {
249        if ( simParams->getPrintTaggedPairDistance()) {
250 <        
250 >
251          std::pair<int, int> tap = simParams->getTaggedAtomPair();
252          Vector3d pos1, pos2, rab;
253  
254 < #ifdef IS_MPI        
254 > #ifdef IS_MPI
255          std::cerr << "tap = " << tap.first << "  " << tap.second << std::endl;
256  
257 <        int mol1 = info_->getGlobalMolMembership(tap.first);
258 <        int mol2 = info_->getGlobalMolMembership(tap.second);
257 >  int mol1 = info_->getGlobalMolMembership(tap.first);
258 >  int mol2 = info_->getGlobalMolMembership(tap.second);
259          std::cerr << "mols = " << mol1 << " " << mol2 << std::endl;
260  
261          int proc1 = info_->getMolToProc(mol1);
# Line 263 | Line 263 | namespace OpenMD {
263  
264          std::cerr << " procs = " << proc1 << " " <<proc2 <<std::endl;
265  
266 <        RealType data[3];
266 >  RealType data[3];
267          if (proc1 == worldRank) {
268            StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first);
269            std::cerr << " on proc " << proc1 << ", sd1 has global index= " << sd1->getGlobalIndex() << std::endl;
270            pos1 = sd1->getPos();
271            data[0] = pos1.x();
272            data[1] = pos1.y();
273 <          data[2] = pos1.z();          
273 >          data[2] = pos1.z();
274            MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
275          } else {
276            MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD);
# Line 284 | Line 284 | namespace OpenMD {
284            pos2 = sd2->getPos();
285            data[0] = pos2.x();
286            data[1] = pos2.y();
287 <          data[2] = pos2.z();          
287 >          data[2] = pos2.z();
288            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
289          } else {
290            MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD);
# Line 295 | Line 295 | namespace OpenMD {
295          StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second);
296          pos1 = at1->getPos();
297          pos2 = at2->getPos();
298 < #endif        
298 > #endif
299          rab = pos2 - pos1;
300          currSnapshot->wrapVector(rab);
301          stat[Stats::TAGGED_PAIR_DISTANCE] =  rab.length();
302        }
303      }
304 <      
304 >
305      /**@todo need refactorying*/
306      //Conserved Quantity is set by integrator and time is set by setTime
307 <    
307 >
308    }
309  
310  
# Line 435 | Line 435 | namespace OpenMD {
435      RealType AvgE_a_ = 0;
436      Vector3d GKappa_t = V3Zero;
437      Vector3d ThermalHelfandMoment;
438 <    
438 >
439      for (mol = info_->beginMolecule(miter); mol != NULL;
440           mol = info_->nextMolecule(miter)) {
441  
# Line 444 | Line 444 | namespace OpenMD {
444  
445          mass = atom->getMass();
446          velocity = atom->getVel();
447 <        kinetic = mass * (velocity[0]*velocity[0] + velocity[1]*velocity[1] +
447 >        kinetic = mass * (velocity[0]*velocity[0] + velocity[1]*velocity[1] +
448                                     velocity[2]*velocity[2]) / PhysicalConstants::energyConvert;
449          potential =  atom->getParticlePot();
450          eatom += (kinetic + potential)/2.0;
# Line 453 | Line 453 | namespace OpenMD {
453  
454     int natoms = info_->getNGlobalAtoms();
455   #ifdef IS_MPI
456 <    
456 >
457      MPI_Allreduce(&eatom, &AvgE_a_, 1, MPI_REALTYPE, MPI_SUM,
458                    MPI_COMM_WORLD);
459 < #else    
459 > #else
460      AvgE_a_ = eatom;
461   #endif
462      AvgE_a_ = AvgE_a_/RealType(natoms);
463 <    
463 >
464      for (mol = info_->beginMolecule(miter); mol != NULL;
465           mol = info_->nextMolecule(miter)) {
466 <      
466 >
467        for (atom = mol->beginAtom(aiter); atom != NULL;
468             atom = mol->nextAtom(aiter)) {
469 <        
469 >
470          /* We think that x_a is relative to the total box and should be a wrapped coordinate */
471          x_a = atom->getPos();
472          currSnapshot->wrapVector(x_a);
473          potential =  atom->getParticlePot();
474 <
475 <        GKappa_t += x_a*(potential-AvgE_a_);
474 >        velocity = atom->getVel();
475 >        kinetic = mass * (velocity[0]*velocity[0] + velocity[1]*velocity[1] +
476 >                           velocity[2]*velocity[2]) / PhysicalConstants::energyConvert;
477 >        eatom += (kinetic + potential)/2.0
478 >        GKappa_t += x_a*(eatom-AvgE_a_);
479          }
480        }
481   #ifdef IS_MPI
482       MPI_Allreduce(GKappa_t.getArrayPointer(), ThermalHelfandMoment.getArrayPointer(), 3,
483                    MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
484 < #else    
484 > #else
485       ThermalHelfandMoment = GKappa_t;
486 < #endif  
486 > #endif
487       return ThermalHelfandMoment;
488 <    
488 >
489   }
490  
491  

Diff Legend

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