[ViewVC] Diff of: OpenMD/branches/development/src/brains/Thermo.cpp

Comparing branches/development/src/brains/Thermo.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1733 by jmichalk, Tue Jun 5 17:48:40 2012 UTC

+ * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
+ * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
-<
+ * [4]  Vardeman & Gezelter, in progress (2009).
->
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
->
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
+ */
+#include <math.h>
+#include "primitives/Molecule.hpp"
+#include "utils/simError.h"
+#include "utils/PhysicalConstants.hpp"
-+
+#include "types/MultipoleAdapter.hpp"
+namespace OpenMD {
+    return temperature;
-+
+  RealType Thermo::getElectronicTemperature() {
-+
+    SimInfo::MoleculeIterator miter;
-+
+    std::vector<Atom*>::iterator iiter;
-+
+    Molecule* mol;
-+
+    Atom* atom;
-+
+    RealType cvel;
-+
+    RealType cmass;
-+
+    RealType kinetic = 0.0;
-+
+    RealType kinetic_global = 0.0;
-+
-+
+    for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) {
-+
+      for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL;
-+
+           atom = mol->nextFluctuatingCharge(iiter)) {
-+
+        cmass = atom->getChargeMass();
-+
+        cvel = atom->getFlucQVel();
-+
-+
+        kinetic += cmass * cvel * cvel;
-+
-+
+      }
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
-+
+    MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM,
-+
+                  MPI_COMM_WORLD);
-+
+    kinetic = kinetic_global;
-+
-+
+#endif //is_mpi
-+
-+
+    kinetic = kinetic * 0.5;
-+
+    return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb );
-+
+  }
-+
-+
-+
-+
+  RealType Thermo::getVolume() {
+    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
+    return curSnapshot->getVolume();
+    RealType volume = this->getVolume();
+    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-<
+    Mat3x3d tau = curSnapshot->statData.getTau();
->
+    Mat3x3d stressTensor = curSnapshot->getStressTensor();
-<
+    pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
->
+    pressureTensor =  (p_global +
->
+                       PhysicalConstants::energyConvert * stressTensor)/volume;
+    return pressureTensor;
+    stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1);
+    stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2);
-+
+    // grab the simulation box dipole moment if specified
-+
+    if (info_->getCalcBoxDipole()){
-+
+      Vector3d totalDipole = getBoxDipole();
-+
+      stat[Stats::BOX_DIPOLE_X] = totalDipole(0);
-+
+      stat[Stats::BOX_DIPOLE_Y] = totalDipole(1);
-+
+      stat[Stats::BOX_DIPOLE_Z] = totalDipole(2);
-+
+    }
+    Globals* simParams = info_->getSimParams();
-+
+    // grab the heat flux if desired
-+
+    if (simParams->havePrintHeatFlux()) {
-+
+      if (simParams->getPrintHeatFlux()){
-+
+        Vector3d heatFlux = getHeatFlux();
-+
+        stat[Stats::HEATFLUX_X] = heatFlux(0);
-+
+        stat[Stats::HEATFLUX_Y] = heatFlux(1);
-+
+        stat[Stats::HEATFLUX_Z] = heatFlux(2);
-+
+      }
-+
+    }
+    if (simParams->haveTaggedAtomPair() &&
+        simParams->havePrintTaggedPairDistance()) {
+    /**@todo need refactorying*/
+    //Conserved Quantity is set by integrator and time is set by setTime
-+
-+
+  }
-+
-+
-+
+  Vector3d Thermo::getBoxDipole() {
-+
+    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-+
+    SimInfo::MoleculeIterator miter;
-+
+    std::vector<Atom*>::iterator aiter;
-+
+    Molecule* mol;
-+
+    Atom* atom;
-+
+    RealType charge;
-+
+    RealType moment(0.0);
-+
+    Vector3d ri(0.0);
-+
+    Vector3d dipoleVector(0.0);
-+
+    Vector3d nPos(0.0);
-+
+    Vector3d pPos(0.0);
-+
+    RealType nChg(0.0);
-+
+    RealType pChg(0.0);
-+
+    int nCount = 0;
-+
+    int pCount = 0;
-+
-+
+    RealType chargeToC = 1.60217733e-19;
-+
+    RealType angstromToM = 1.0e-10;
-+
+    RealType debyeToCm = 3.33564095198e-30;
-+
+    for (mol = info_->beginMolecule(miter); mol != NULL;
-+
+         mol = info_->nextMolecule(miter)) {
-+
-+
+      for (atom = mol->beginAtom(aiter); atom != NULL;
-+
+           atom = mol->nextAtom(aiter)) {
-+
-+
+        if (atom->isCharge() ) {
-+
+          charge = 0.0;
-+
+          GenericData* data = atom->getAtomType()->getPropertyByName("Charge");
-+
+          if (data != NULL) {
-+
-+
+            charge = (dynamic_cast<DoubleGenericData*>(data))->getData();
-+
+            charge *= chargeToC;
-+
-+
+            ri = atom->getPos();
-+
+            currSnapshot->wrapVector(ri);
-+
+            ri *= angstromToM;
-+
-+
+            if (charge < 0.0) {
-+
+              nPos += ri;
-+
+              nChg -= charge;
-+
+              nCount++;
-+
+            } else if (charge > 0.0) {
-+
+              pPos += ri;
-+
+              pChg += charge;
-+
+              pCount++;
-+
+            }
-+
+          }
-+
+        }
-+
-+
+        MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
-+
+        if (ma.isDipole() ) {
-+
+          Vector3d u_i = atom->getElectroFrame().getColumn(2);
-+
+          moment = ma.getDipoleMoment();
-+
+          moment *= debyeToCm;
-+
+          dipoleVector += u_i * moment;
-+
+        }
-+
+      }
-+
+    }
-+
-+
-+
+#ifdef IS_MPI
-+
+    RealType pChg_global, nChg_global;
-+
+    int pCount_global, nCount_global;
-+
+    Vector3d pPos_global, nPos_global, dipVec_global;
-+
-+
+    MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM,
-+
+                  MPI_COMM_WORLD);
-+
+    pChg = pChg_global;
-+
+    MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM,
-+
+                  MPI_COMM_WORLD);
-+
+    nChg = nChg_global;
-+
+    MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM,
-+
+                  MPI_COMM_WORLD);
-+
+    pCount = pCount_global;
-+
+    MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM,
-+
+                  MPI_COMM_WORLD);
-+
+    nCount = nCount_global;
-+
+    MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3,
-+
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
-+
+    pPos = pPos_global;
-+
+    MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3,
-+
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
-+
+    nPos = nPos_global;
-+
+    MPI_Allreduce(dipoleVector.getArrayPointer(),
-+
+                  dipVec_global.getArrayPointer(), 3,
-+
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
-+
+    dipoleVector = dipVec_global;
-+
+#endif //is_mpi
-+
-+
+    // first load the accumulated dipole moment (if dipoles were present)
-+
+    Vector3d boxDipole = dipoleVector;
-+
+    // now include the dipole moment due to charges
-+
+    // use the lesser of the positive and negative charge totals
-+
+    RealType chg_value = nChg <= pChg ? nChg : pChg;
-+
-+
+    // find the average positions
-+
+    if (pCount > 0 && nCount > 0 ) {
-+
+      pPos /= pCount;
-+
+      nPos /= nCount;
-+
+    }
-+
-+
+    // dipole is from the negative to the positive (physics notation)
-+
+    boxDipole += (pPos - nPos) * chg_value;
-+
-+
+    return boxDipole;
-+
+  // Returns the Heat Flux Vector for the system
-+
+  Vector3d Thermo::getHeatFlux(){
-+
+    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-+
+    SimInfo::MoleculeIterator miter;
-+
+    std::vector<StuntDouble*>::iterator iiter;
-+
+    Molecule* mol;
-+
+    StuntDouble* integrableObject;
-+
+    RigidBody::AtomIterator ai;
-+
+    Atom* atom;
-+
+    Vector3d vel;
-+
+    Vector3d angMom;
-+
+    Mat3x3d I;
-+
+    int i;
-+
+    int j;
-+
+    int k;
-+
+    RealType mass;
-+
-+
+    Vector3d x_a;
-+
+    RealType kinetic;
-+
+    RealType potential;
-+
+    RealType eatom;
-+
+    RealType AvgE_a_ = 0;
-+
+    // Convective portion of the heat flux
-+
+    Vector3d heatFluxJc = V3Zero;
-+
-+
+    /* Calculate convective portion of the heat flux */
-+
+    for (mol = info_->beginMolecule(miter); mol != NULL;
-+
+         mol = info_->nextMolecule(miter)) {
-+
-+
+      for (integrableObject = mol->beginIntegrableObject(iiter);
-+
+           integrableObject != NULL;
-+
+           integrableObject = mol->nextIntegrableObject(iiter)) {
-+
-+
+        mass = integrableObject->getMass();
-+
+        vel = integrableObject->getVel();
-+
-+
+        kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]);
-+
-+
+        if (integrableObject->isDirectional()) {
-+
+          angMom = integrableObject->getJ();
-+
+          I = integrableObject->getI();
-+
-+
+          if (integrableObject->isLinear()) {
-+
+            i = integrableObject->linearAxis();
-+
+            j = (i + 1) % 3;
-+
+            k = (i + 2) % 3;
-+
+            kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
-+
+          } else {
-+
+            kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1)
-+
+              + angMom[2]*angMom[2]/I(2, 2);
-+
+          }
-+
+        }
-+
-+
+        potential = 0.0;
-+
-+
+        if (integrableObject->isRigidBody()) {
-+
+          RigidBody* rb = dynamic_cast<RigidBody*>(integrableObject);
-+
+          for (atom = rb->beginAtom(ai); atom != NULL;
-+
+               atom = rb->nextAtom(ai)) {
-+
+            potential +=  atom->getParticlePot();
-+
+          }
-+
+        } else {
-+
+          potential = integrableObject->getParticlePot();
-+
+          cerr << "ppot = "  << potential << "\n";
-+
+        }
-+
-+
+        potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2
-+
+        // The potential may not be a 1/2 factor
-+
+        eatom = (kinetic + potential)/2.0;  // amu A^2/fs^2
-+
+        heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3
-+
+        heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3
-+
+        heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3
-+
+      }
-+
+    }
-+
-+
+    std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl;
-+
-+
+    /* The J_v vector is reduced in fortan so everyone has the global
-+
+     *  Jv. Jc is computed over the local atoms and must be reduced
-+
+     *  among all processors.
-+
+     */
-+
+#ifdef IS_MPI
-+
+    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE,
-+
+                              MPI::SUM);
-+
+#endif
-+
-+
+    // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3
-+
-+
+    Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() *
-+
+      PhysicalConstants::energyConvert;
-+
-+
+    std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl;
-+
+    std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl;
-+
-+
+    // Correct for the fact the flux is 1/V (Jc + Jv)
-+
+    return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3
-+
+  }
+} //end namespace OpenMD

Diff Legend

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

Comparing branches/development/src/brains/Thermo.cpp (file contents): Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs. Revision 1733 by jmichalk, Tue Jun 5 17:48:40 2012 UTC

Diff Legend

Comparing branches/development/src/brains/Thermo.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1733 by jmichalk, Tue Jun 5 17:48:40 2012 UTC