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

Comparing branches/development/src/brains/Thermo.cpp (file contents):
Revision 1503 by gezelter, Sat Oct 2 19:54:41 2010 UTC vs.
Revision 1723 by gezelter, Thu May 24 20:59:54 2012 UTC

# Line 36 | Line 36
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).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <math.h>
# Line 50 | Line 51
51   #include "primitives/Molecule.hpp"
52   #include "utils/simError.h"
53   #include "utils/PhysicalConstants.hpp"
54 + #include "types/MultipoleAdapter.hpp"
55  
56   namespace OpenMD {
57  
# Line 143 | Line 145 | namespace OpenMD {
145      return temperature;
146    }
147  
148 +  RealType Thermo::getElectronicTemperature() {
149 +    SimInfo::MoleculeIterator miter;
150 +    std::vector<Atom*>::iterator iiter;
151 +    Molecule* mol;
152 +    Atom* atom;    
153 +    RealType cvel;
154 +    RealType cmass;
155 +    RealType kinetic = 0.0;
156 +    RealType kinetic_global = 0.0;
157 +    
158 +    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
159 +      for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
160 +           atom = mol->nextFluctuatingCharge(iiter)) {
161 +        cmass = atom->getChargeMass();
162 +        cvel = atom->getFlucQVel();
163 +        
164 +        kinetic += cmass * cvel * cvel;
165 +        
166 +      }
167 +    }
168 +    
169 + #ifdef IS_MPI
170 +
171 +    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
172 +                  MPI_COMM_WORLD);
173 +    kinetic = kinetic_global;
174 +
175 + #endif //is_mpi
176 +
177 +    kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert;
178 +    return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb );    
179 +  }
180 +
181 +
182 +
183 +
184    RealType Thermo::getVolume() {
185      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
186      return curSnapshot->getVolume();
# Line 208 | Line 246 | namespace OpenMD {
246  
247      RealType volume = this->getVolume();
248      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
249 <    Mat3x3d tau = curSnapshot->statData.getTau();
249 >    Mat3x3d stressTensor = curSnapshot->getStressTensor();
250  
251 <    pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
251 >    pressureTensor =  (p_global +
252 >                       PhysicalConstants::energyConvert * stressTensor)/volume;
253      
254      return pressureTensor;
255    }
# Line 247 | Line 286 | namespace OpenMD {
286      }
287  
288      Globals* simParams = info_->getSimParams();
289 +    // grab the heat flux if desired
290 +    if (simParams->havePrintHeatFlux()) {
291 +      if (simParams->getPrintHeatFlux()){
292 +        Vector3d heatFlux = getHeatFlux();
293 +        stat[Stats::HEATFLUX_X] = heatFlux(0);
294 +        stat[Stats::HEATFLUX_Y] = heatFlux(1);
295 +        stat[Stats::HEATFLUX_Z] = heatFlux(2);
296 +      }
297 +    }
298  
299      if (simParams->haveTaggedAtomPair() &&
300          simParams->havePrintTaggedPairDistance()) {
# Line 363 | Line 411 | namespace OpenMD {
411            }
412          }
413          
414 <        if (atom->isDipole() ) {
414 >        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
415 >        if (ma.isDipole() ) {
416            Vector3d u_i = atom->getElectroFrame().getColumn(2);
417 <          GenericData* data = dynamic_cast<DirectionalAtomType*>(atom->getAtomType())->getPropertyByName("Dipole");
418 <          if (data != NULL) {
419 <            moment = (dynamic_cast<DoubleGenericData*>(data))->getData();
371 <            
372 <            moment *= debyeToCm;
373 <            dipoleVector += u_i * moment;
374 <          }
417 >          moment = ma.getDipoleMoment();
418 >          moment *= debyeToCm;
419 >          dipoleVector += u_i * moment;
420          }
421        }
422      }
# Line 423 | Line 468 | namespace OpenMD {
468  
469      return boxDipole;
470    }
471 +
472 +  // Returns the Heat Flux Vector for the system
473 +  Vector3d Thermo::getHeatFlux(){
474 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
475 +    SimInfo::MoleculeIterator miter;
476 +    std::vector<StuntDouble*>::iterator iiter;
477 +    Molecule* mol;
478 +    StuntDouble* integrableObject;    
479 +    RigidBody::AtomIterator ai;
480 +    Atom* atom;      
481 +    Vector3d vel;
482 +    Vector3d angMom;
483 +    Mat3x3d I;
484 +    int i;
485 +    int j;
486 +    int k;
487 +    RealType mass;
488 +
489 +    Vector3d x_a;
490 +    RealType kinetic;
491 +    RealType potential;
492 +    RealType eatom;
493 +    RealType AvgE_a_ = 0;
494 +    // Convective portion of the heat flux
495 +    Vector3d heatFluxJc = V3Zero;
496 +
497 +    /* Calculate convective portion of the heat flux */
498 +    for (mol = info_->beginMolecule(miter); mol != NULL;
499 +         mol = info_->nextMolecule(miter)) {
500 +      
501 +      for (integrableObject = mol->beginIntegrableObject(iiter);
502 +           integrableObject != NULL;
503 +           integrableObject = mol->nextIntegrableObject(iiter)) {
504 +        
505 +        mass = integrableObject->getMass();
506 +        vel = integrableObject->getVel();
507 +
508 +        kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
509 +        
510 +        if (integrableObject->isDirectional()) {
511 +          angMom = integrableObject->getJ();
512 +          I = integrableObject->getI();
513 +
514 +          if (integrableObject->isLinear()) {
515 +            i = integrableObject->linearAxis();
516 +            j = (i + 1) % 3;
517 +            k = (i + 2) % 3;
518 +            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
519 +          } else {                        
520 +            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
521 +              + angMom[2]*angMom[2]/I(2, 2);
522 +          }
523 +        }
524 +
525 +        potential = 0.0;
526 +
527 +        if (integrableObject->isRigidBody()) {
528 +          RigidBody* rb = dynamic_cast<RigidBody*>(integrableObject);
529 +          for (atom = rb->beginAtom(ai); atom != NULL;
530 +               atom = rb->nextAtom(ai)) {
531 +            potential +=  atom->getParticlePot();
532 +          }          
533 +        } else {
534 +          potential = integrableObject->getParticlePot();
535 +          cerr << "ppot = "  << potential << "\n";
536 +        }
537 +
538 +        potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
539 +        // The potential may not be a 1/2 factor
540 +        eatom = (kinetic + potential)/2.0;  // amu A^2/fs^2
541 +        heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3
542 +        heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3
543 +        heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3
544 +      }
545 +    }
546 +
547 +    std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl;
548 +
549 +    /* The J_v vector is reduced in fortan so everyone has the global
550 +     *  Jv. Jc is computed over the local atoms and must be reduced
551 +     *  among all processors.
552 +     */
553 + #ifdef IS_MPI
554 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
555 +                              MPI::SUM);
556 + #endif
557 +    
558 +    // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
559 +
560 +    Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
561 +      PhysicalConstants::energyConvert;
562 +    
563 +    std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl;
564 +    std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl;
565 +    
566 +    // Correct for the fact the flux is 1/V (Jc + Jv)
567 +    return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
568 +  }
569   } //end namespace OpenMD

Diff Legend

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