[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 1710 by gezelter, Fri May 18 21:44:02 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 {
+    RealType volume = this->getVolume();
+    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
-<
+    Mat3x3d tau = curSnapshot->statData.getTau();
->
+    Mat3x3d tau = curSnapshot->getTau();
+    pressureTensor =  (p_global + PhysicalConstants::energyConvert* tau)/volume;
+    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();
+    //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;
-+
+  }
+} //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 1710 by gezelter, Fri May 18 21:44:02 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 1710 by gezelter, Fri May 18 21:44:02 2012 UTC