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 1709 by gezelter, Tue May 15 13:04:08 2012 UTC vs.
Revision 1760 by gezelter, Thu Jun 21 19:26:46 2012 UTC

# Line 51 | 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 111 | Line 112 | namespace OpenMD {
112  
113    RealType Thermo::getPotential() {
114      RealType potential = 0.0;
114 –    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115 –    RealType shortRangePot_local =  curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ;
115  
116 <    // Get total potential for entire system from MPI.
117 <
119 < #ifdef IS_MPI
120 <
121 <    MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM,
122 <                  MPI_COMM_WORLD);
123 <    potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
124 <
125 < #else
126 <
127 <    potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
128 <
129 < #endif // is_mpi
130 <
116 >    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
117 >    potential = curSnapshot->getShortRangePotential() + curSnapshot->getLongRangePotential();
118      return potential;
119    }
120  
# Line 142 | Line 129 | namespace OpenMD {
129      
130      RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb );
131      return temperature;
132 +  }
133 +
134 +  RealType Thermo::getElectronicTemperature() {
135 +    SimInfo::MoleculeIterator miter;
136 +    std::vector<Atom*>::iterator iiter;
137 +    Molecule* mol;
138 +    Atom* atom;    
139 +    RealType cvel;
140 +    RealType cmass;
141 +    RealType kinetic = 0.0;
142 +    RealType kinetic_global = 0.0;
143 +    
144 +    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
145 +      for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
146 +           atom = mol->nextFluctuatingCharge(iiter)) {
147 +        cmass = atom->getChargeMass();
148 +        cvel = atom->getFlucQVel();
149 +        
150 +        kinetic += cmass * cvel * cvel;
151 +        
152 +      }
153 +    }
154 +    
155 + #ifdef IS_MPI
156 +
157 +    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
158 +                  MPI_COMM_WORLD);
159 +    kinetic = kinetic_global;
160 +
161 + #endif //is_mpi
162 +
163 +    kinetic = kinetic * 0.5;
164 +    return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb );    
165    }
166  
167 +
168 +
169 +
170    RealType Thermo::getVolume() {
171      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
172      return curSnapshot->getVolume();
# Line 209 | Line 232 | namespace OpenMD {
232  
233      RealType volume = this->getVolume();
234      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
235 <    Mat3x3d tau = curSnapshot->getTau();
235 >    Mat3x3d stressTensor = curSnapshot->getStressTensor();
236  
237 <    pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
237 >    pressureTensor =  (p_global +
238 >                       PhysicalConstants::energyConvert * stressTensor)/volume;
239      
240      return pressureTensor;
241    }
# Line 248 | Line 272 | namespace OpenMD {
272      }
273  
274      Globals* simParams = info_->getSimParams();
275 +    // grab the heat flux if desired
276 +    if (simParams->havePrintHeatFlux()) {
277 +      if (simParams->getPrintHeatFlux()){
278 +        Vector3d heatFlux = getHeatFlux();
279 +        stat[Stats::HEATFLUX_X] = heatFlux(0);
280 +        stat[Stats::HEATFLUX_Y] = heatFlux(1);
281 +        stat[Stats::HEATFLUX_Z] = heatFlux(2);
282 +      }
283 +    }
284  
285      if (simParams->haveTaggedAtomPair() &&
286          simParams->havePrintTaggedPairDistance()) {
# Line 364 | Line 397 | namespace OpenMD {
397            }
398          }
399          
400 <        if (atom->isDipole() ) {
400 >        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
401 >        if (ma.isDipole() ) {
402            Vector3d u_i = atom->getElectroFrame().getColumn(2);
403 <          GenericData* data = dynamic_cast<DirectionalAtomType*>(atom->getAtomType())->getPropertyByName("Dipole");
404 <          if (data != NULL) {
405 <            moment = (dynamic_cast<DoubleGenericData*>(data))->getData();
372 <            
373 <            moment *= debyeToCm;
374 <            dipoleVector += u_i * moment;
375 <          }
403 >          moment = ma.getDipoleMoment();
404 >          moment *= debyeToCm;
405 >          dipoleVector += u_i * moment;
406          }
407        }
408      }
# Line 423 | Line 453 | namespace OpenMD {
453      boxDipole += (pPos - nPos) * chg_value;
454  
455      return boxDipole;
456 +  }
457 +
458 +  // Returns the Heat Flux Vector for the system
459 +  Vector3d Thermo::getHeatFlux(){
460 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
461 +    SimInfo::MoleculeIterator miter;
462 +    std::vector<StuntDouble*>::iterator iiter;
463 +    Molecule* mol;
464 +    StuntDouble* integrableObject;    
465 +    RigidBody::AtomIterator ai;
466 +    Atom* atom;      
467 +    Vector3d vel;
468 +    Vector3d angMom;
469 +    Mat3x3d I;
470 +    int i;
471 +    int j;
472 +    int k;
473 +    RealType mass;
474 +
475 +    Vector3d x_a;
476 +    RealType kinetic;
477 +    RealType potential;
478 +    RealType eatom;
479 +    RealType AvgE_a_ = 0;
480 +    // Convective portion of the heat flux
481 +    Vector3d heatFluxJc = V3Zero;
482 +
483 +    /* Calculate convective portion of the heat flux */
484 +    for (mol = info_->beginMolecule(miter); mol != NULL;
485 +         mol = info_->nextMolecule(miter)) {
486 +      
487 +      for (integrableObject = mol->beginIntegrableObject(iiter);
488 +           integrableObject != NULL;
489 +           integrableObject = mol->nextIntegrableObject(iiter)) {
490 +        
491 +        mass = integrableObject->getMass();
492 +        vel = integrableObject->getVel();
493 +
494 +        kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
495 +        
496 +        if (integrableObject->isDirectional()) {
497 +          angMom = integrableObject->getJ();
498 +          I = integrableObject->getI();
499 +
500 +          if (integrableObject->isLinear()) {
501 +            i = integrableObject->linearAxis();
502 +            j = (i + 1) % 3;
503 +            k = (i + 2) % 3;
504 +            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
505 +          } else {                        
506 +            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
507 +              + angMom[2]*angMom[2]/I(2, 2);
508 +          }
509 +        }
510 +
511 +        potential = 0.0;
512 +
513 +        if (integrableObject->isRigidBody()) {
514 +          RigidBody* rb = dynamic_cast<RigidBody*>(integrableObject);
515 +          for (atom = rb->beginAtom(ai); atom != NULL;
516 +               atom = rb->nextAtom(ai)) {
517 +            potential +=  atom->getParticlePot();
518 +          }          
519 +        } else {
520 +          potential = integrableObject->getParticlePot();
521 +          cerr << "ppot = "  << potential << "\n";
522 +        }
523 +
524 +        potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
525 +        // The potential may not be a 1/2 factor
526 +        eatom = (kinetic + potential)/2.0;  // amu A^2/fs^2
527 +        heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3
528 +        heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3
529 +        heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3
530 +      }
531 +    }
532 +
533 +    std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl;
534 +
535 +    /* The J_v vector is reduced in fortan so everyone has the global
536 +     *  Jv. Jc is computed over the local atoms and must be reduced
537 +     *  among all processors.
538 +     */
539 + #ifdef IS_MPI
540 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
541 +                              MPI::SUM);
542 + #endif
543 +    
544 +    // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
545 +
546 +    Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
547 +      PhysicalConstants::energyConvert;
548 +    
549 +    std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl;
550 +    std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl;
551 +    
552 +    // Correct for the fact the flux is 1/V (Jc + Jv)
553 +    return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
554    }
555   } //end namespace OpenMD

Diff Legend

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