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 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/brains/Thermo.cpp (file contents), 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 238 | Line 277 | namespace OpenMD {
277      stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);      
278      stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);      
279  
280 +    // grab the simulation box dipole moment if specified
281 +    if (info_->getCalcBoxDipole()){
282 +      Vector3d totalDipole = getBoxDipole();
283 +      stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
284 +      stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
285 +      stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
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 301 | Line 356 | namespace OpenMD {
356        
357      /**@todo need refactorying*/
358      //Conserved Quantity is set by integrator and time is set by setTime
359 +    
360 +  }
361 +
362 +
363 +  Vector3d Thermo::getBoxDipole() {
364 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
365 +    SimInfo::MoleculeIterator miter;
366 +    std::vector<Atom*>::iterator aiter;
367 +    Molecule* mol;
368 +    Atom* atom;
369 +    RealType charge;
370 +    RealType moment(0.0);
371 +    Vector3d ri(0.0);
372 +    Vector3d dipoleVector(0.0);
373 +    Vector3d nPos(0.0);
374 +    Vector3d pPos(0.0);
375 +    RealType nChg(0.0);
376 +    RealType pChg(0.0);
377 +    int nCount = 0;
378 +    int pCount = 0;
379 +
380 +    RealType chargeToC = 1.60217733e-19;
381 +    RealType angstromToM = 1.0e-10;
382 +    RealType debyeToCm = 3.33564095198e-30;
383      
384 +    for (mol = info_->beginMolecule(miter); mol != NULL;
385 +         mol = info_->nextMolecule(miter)) {
386 +
387 +      for (atom = mol->beginAtom(aiter); atom != NULL;
388 +           atom = mol->nextAtom(aiter)) {
389 +        
390 +        if (atom->isCharge() ) {
391 +          charge = 0.0;
392 +          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
393 +          if (data != NULL) {
394 +
395 +            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
396 +            charge *= chargeToC;
397 +
398 +            ri = atom->getPos();
399 +            currSnapshot->wrapVector(ri);
400 +            ri *= angstromToM;
401 +
402 +            if (charge < 0.0) {
403 +              nPos += ri;
404 +              nChg -= charge;
405 +              nCount++;
406 +            } else if (charge > 0.0) {
407 +              pPos += ri;
408 +              pChg += charge;
409 +              pCount++;
410 +            }                      
411 +          }
412 +        }
413 +        
414 +        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
415 +        if (ma.isDipole() ) {
416 +          Vector3d u_i = atom->getElectroFrame().getColumn(2);
417 +          moment = ma.getDipoleMoment();
418 +          moment *= debyeToCm;
419 +          dipoleVector += u_i * moment;
420 +        }
421 +      }
422 +    }
423 +    
424 +                      
425 + #ifdef IS_MPI
426 +    RealType pChg_global, nChg_global;
427 +    int pCount_global, nCount_global;
428 +    Vector3d pPos_global, nPos_global, dipVec_global;
429 +
430 +    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
431 +                  MPI_COMM_WORLD);
432 +    pChg = pChg_global;
433 +    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
434 +                  MPI_COMM_WORLD);
435 +    nChg = nChg_global;
436 +    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
437 +                  MPI_COMM_WORLD);
438 +    pCount = pCount_global;
439 +    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
440 +                  MPI_COMM_WORLD);
441 +    nCount = nCount_global;
442 +    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
443 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
444 +    pPos = pPos_global;
445 +    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
446 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
447 +    nPos = nPos_global;
448 +    MPI_Allreduce(dipoleVector.getArrayPointer(),
449 +                  dipVec_global.getArrayPointer(), 3,
450 +                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
451 +    dipoleVector = dipVec_global;
452 + #endif //is_mpi
453 +
454 +    // first load the accumulated dipole moment (if dipoles were present)
455 +    Vector3d boxDipole = dipoleVector;
456 +    // now include the dipole moment due to charges
457 +    // use the lesser of the positive and negative charge totals
458 +    RealType chg_value = nChg <= pChg ? nChg : pChg;
459 +      
460 +    // find the average positions
461 +    if (pCount > 0 && nCount > 0 ) {
462 +      pPos /= pCount;
463 +      nPos /= nCount;
464 +    }
465 +
466 +    // dipole is from the negative to the positive (physics notation)
467 +    boxDipole += (pPos - nPos) * chg_value;
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

Comparing:
trunk/src/brains/Thermo.cpp (property svn:keywords), Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
branches/development/src/brains/Thermo.cpp (property svn:keywords), Revision 1723 by gezelter, Thu May 24 20:59:54 2012 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

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