ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/rnemd/RNEMD.cpp
(Generate patch)

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 1946 by gezelter, Tue Nov 12 02:18:35 2013 UTC

# Line 38 | Line 38
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41 + #ifdef IS_MPI
42 + #include <mpi.h>
43 + #endif
44  
45   #include <cmath>
46   #include <sstream>
# Line 54 | Line 57
57   #include "utils/Tuple.hpp"
58   #include "brains/Thermo.hpp"
59   #include "math/ConvexHull.hpp"
57 – #ifdef IS_MPI
58 – #include <mpi.h>
59 – #endif
60  
61   #ifdef _MSC_VER
62   #define isnan(x) _isnan((x))
# Line 547 | Line 547 | namespace OpenMD {
547      
548        if (hasSelectionB_) {
549          selectionB_ = rnemdParams->getSelectionB();
550 +
551        } else {
552          if (usePeriodicBoundaryConditions_) {    
553            Mat3x3d hmat = currentSnap_->getHmat();
# Line 573 | Line 574 | namespace OpenMD {
574              selectionB_ = selectionBstream.str();
575            } else {
576              selectionB_ = "select hull";
577 +            BisHull_ = true;
578              hasSelectionB_ = true;
579            }
580          }
# Line 1729 | Line 1731 | namespace OpenMD {
1731      Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1732  
1733      if (hasSelectionA_) {
1734 <      int isd;
1735 <      StuntDouble* sd;
1736 <      vector<StuntDouble*> aSites;
1737 <      seleManA_.setSelectionSet(evaluatorA_.evaluate());
1738 <      for (sd = seleManA_.beginSelected(isd); sd != NULL;
1739 <           sd = seleManA_.nextSelected(isd)) {
1740 <        aSites.push_back(sd);
1741 <      }
1742 < #if defined(HAVE_QHULL)
1743 <      ConvexHull* surfaceMeshA = new ConvexHull();
1744 <      surfaceMeshA->computeHull(aSites);
1745 <      areaA = surfaceMeshA->getArea();
1746 <      delete surfaceMeshA;
1747 < #else
1748 <      sprintf( painCave.errMsg,
1749 <               "RNEMD::getDividingArea : Hull calculation is not possible\n"
1750 <               "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1751 <      painCave.severity = OPENMD_ERROR;
1752 <      painCave.isFatal = 1;
1753 <      simError();
1734 >
1735 >      if (evaluatorA_.hasSurfaceArea())
1736 >        areaA = evaluatorA_.getSurfaceArea();
1737 >      else {
1738 >        
1739 >        cerr << "selection A did not have surface area, recomputing\n";
1740 >        int isd;
1741 >        StuntDouble* sd;
1742 >        vector<StuntDouble*> aSites;
1743 >        seleManA_.setSelectionSet(evaluatorA_.evaluate());
1744 >        for (sd = seleManA_.beginSelected(isd); sd != NULL;
1745 >             sd = seleManA_.nextSelected(isd)) {
1746 >          aSites.push_back(sd);
1747 >        }
1748 > #if defined(HAVE_QHULL)
1749 >        ConvexHull* surfaceMeshA = new ConvexHull();
1750 >        surfaceMeshA->computeHull(aSites);
1751 >        areaA = surfaceMeshA->getArea();
1752 >        delete surfaceMeshA;
1753 > #else
1754 >        sprintf( painCave.errMsg,
1755 >               "RNEMD::getDividingArea : Hull calculation is not possible\n"
1756 >                 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1757 >        painCave.severity = OPENMD_ERROR;
1758 >        painCave.isFatal = 1;
1759 >        simError();
1760   #endif
1761 +      }
1762  
1763      } else {
1764        if (usePeriodicBoundaryConditions_) {
# Line 1765 | Line 1774 | namespace OpenMD {
1774      }
1775  
1776      if (hasSelectionB_) {
1777 <      int isd;
1778 <      StuntDouble* sd;
1779 <      vector<StuntDouble*> bSites;
1780 <      seleManB_.setSelectionSet(evaluatorB_.evaluate());
1772 <      for (sd = seleManB_.beginSelected(isd); sd != NULL;
1773 <           sd = seleManB_.nextSelected(isd)) {
1774 <        bSites.push_back(sd);
1775 <      }
1777 >      if (evaluatorB_.hasSurfaceArea())
1778 >        areaB = evaluatorB_.getSurfaceArea();
1779 >      else {
1780 >        cerr << "selection B did not have surface area, recomputing\n";
1781  
1782 +        int isd;
1783 +        StuntDouble* sd;
1784 +        vector<StuntDouble*> bSites;
1785 +        seleManB_.setSelectionSet(evaluatorB_.evaluate());
1786 +        for (sd = seleManB_.beginSelected(isd); sd != NULL;
1787 +             sd = seleManB_.nextSelected(isd)) {
1788 +          bSites.push_back(sd);
1789 +        }
1790 +        
1791   #if defined(HAVE_QHULL)
1792 <      ConvexHull* surfaceMeshB = new ConvexHull();    
1793 <      surfaceMeshB->computeHull(bSites);
1794 <      areaB = surfaceMeshB->getArea();
1795 <      delete surfaceMeshB;
1792 >        ConvexHull* surfaceMeshB = new ConvexHull();    
1793 >        surfaceMeshB->computeHull(bSites);
1794 >        areaB = surfaceMeshB->getArea();
1795 >        delete surfaceMeshB;
1796   #else
1797 <      sprintf( painCave.errMsg,
1798 <               "RNEMD::getDividingArea : Hull calculation is not possible\n"
1799 <               "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1800 <      painCave.severity = OPENMD_ERROR;
1801 <      painCave.isFatal = 1;
1802 <      simError();
1797 >        sprintf( painCave.errMsg,
1798 >                 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1799 >                 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1800 >        painCave.severity = OPENMD_ERROR;
1801 >        painCave.isFatal = 1;
1802 >        simError();
1803   #endif
1804 <
1805 <
1804 >      }
1805 >      
1806      } else {
1807        if (usePeriodicBoundaryConditions_) {
1808          // in periodic boundaries, the surface area is twice the x-y
# Line 1800 | Line 1814 | namespace OpenMD {
1814          areaB = 4.0 * M_PI * pow(sphereBRadius_, 2);
1815        }
1816      }
1817 <    
1817 >      
1818      dividingArea_ = min(areaA, areaB);
1819      hasDividingArea_ = true;
1820      return dividingArea_;
# Line 1860 | Line 1874 | namespace OpenMD {
1874      RealType area = getDividingArea();
1875      areaAccumulator_->add(area);
1876      Mat3x3d hmat = currentSnap_->getHmat();
1877 +    Vector3d u = angularMomentumFluxVector_;
1878 +    u.normalize();
1879 +
1880      seleMan_.setSelectionSet(evaluator_.evaluate());
1881  
1882      int selei(0);
1883      StuntDouble* sd;
1884      int binNo;
1885 +    RealType mass;
1886 +    Vector3d vel;
1887 +    Vector3d rPos;
1888 +    RealType KE;
1889 +    Vector3d L;
1890 +    Mat3x3d I;
1891 +    RealType r2;
1892  
1893      vector<RealType> binMass(nBins_, 0.0);
1894 <    vector<RealType> binPx(nBins_, 0.0);
1895 <    vector<RealType> binPy(nBins_, 0.0);
1896 <    vector<RealType> binPz(nBins_, 0.0);
1897 <    vector<RealType> binOmegax(nBins_, 0.0);
1874 <    vector<RealType> binOmegay(nBins_, 0.0);
1875 <    vector<RealType> binOmegaz(nBins_, 0.0);
1894 >    vector<Vector3d> binP(nBins_, V3Zero);
1895 >    vector<RealType> binOmega(nBins_, 0.0);
1896 >    vector<Vector3d> binL(nBins_, V3Zero);
1897 >    vector<Mat3x3d>  binI(nBins_);
1898      vector<RealType> binKE(nBins_, 0.0);
1899      vector<int> binDOF(nBins_, 0);
1900      vector<int> binCount(nBins_, 0);
# Line 1912 | Line 1934 | namespace OpenMD {
1934          binNo = int(rPos.length() / binWidth_);
1935        }
1936  
1937 <      RealType mass = sd->getMass();
1938 <      Vector3d vel = sd->getVel();
1939 <      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1940 <      Vector3d aVel = cross(rPos, vel);
1941 <      
1937 >      mass = sd->getMass();
1938 >      vel = sd->getVel();
1939 >      rPos = sd->getPos() - coordinateOrigin_;
1940 >      KE = 0.5 * mass * vel.lengthSquare();
1941 >      L = mass * cross(rPos, vel);
1942 >      I = outProduct(rPos, rPos) * mass;
1943 >      r2 = rPos.lengthSquare();
1944 >      I(0, 0) += mass * r2;
1945 >      I(1, 1) += mass * r2;
1946 >      I(2, 2) += mass * r2;
1947 >
1948 >      // Project the relative position onto a plane perpendicular to
1949 >      // the angularMomentumFluxVector:
1950 >      // Vector3d rProj = rPos - dot(rPos, u) * u;
1951 >      // Project the velocity onto a plane perpendicular to the
1952 >      // angularMomentumFluxVector:
1953 >      // Vector3d vProj = vel  - dot(vel, u) * u;
1954 >      // Compute angular velocity vector (should be nearly parallel to
1955 >      // angularMomentumFluxVector
1956 >      // Vector3d aVel = cross(rProj, vProj);
1957 >
1958        if (binNo >= 0 && binNo < nBins_)  {
1959          binCount[binNo]++;
1960          binMass[binNo] += mass;
1961 <        binPx[binNo] += mass*vel.x();
1962 <        binPy[binNo] += mass*vel.y();
1963 <        binPz[binNo] += mass*vel.z();
1964 <        binOmegax[binNo] += aVel.x();
1927 <        binOmegay[binNo] += aVel.y();
1928 <        binOmegaz[binNo] += aVel.z();
1929 <        binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
1961 >        binP[binNo] += mass*vel;
1962 >        binKE[binNo] += KE;
1963 >        binI[binNo] += I;
1964 >        binL[binNo] += L;
1965          binDOF[binNo] += 3;
1966          
1967          if (sd->isDirectional()) {
1968            Vector3d angMom = sd->getJ();
1969 <          Mat3x3d I = sd->getI();
1969 >          Mat3x3d Ia = sd->getI();
1970            if (sd->isLinear()) {
1971              int i = sd->linearAxis();
1972              int j = (i + 1) % 3;
1973              int k = (i + 2) % 3;
1974 <            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
1975 <                                   angMom[k] * angMom[k] / I(k, k));
1974 >            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
1975 >                                   angMom[k] * angMom[k] / Ia(k, k));
1976              binDOF[binNo] += 2;
1977            } else {
1978 <            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
1979 <                                   angMom[1] * angMom[1] / I(1, 1) +
1980 <                                   angMom[2] * angMom[2] / I(2, 2));
1978 >            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
1979 >                                   angMom[1] * angMom[1] / Ia(1, 1) +
1980 >                                   angMom[2] * angMom[2] / Ia(2, 2));
1981              binDOF[binNo] += 3;
1982            }
1983          }
# Line 1950 | Line 1985 | namespace OpenMD {
1985      }
1986      
1987   #ifdef IS_MPI
1988 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1989 <                              nBins_, MPI::INT, MPI::SUM);
1990 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0],
1991 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1992 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0],
1993 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1994 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0],
1995 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1996 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
1997 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1998 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
1999 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2000 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
2001 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2002 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
2003 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2004 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
2005 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2006 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0],
2007 <                              nBins_, MPI::INT, MPI::SUM);
1988 >
1989 >    for (int i = 0; i < nBins_; i++) {
1990 >
1991 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[i],
1992 >                                1, MPI::INT, MPI::SUM);
1993 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[i],
1994 >                                1, MPI::REALTYPE, MPI::SUM);
1995 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binP[i],
1996 >                                3, MPI::REALTYPE, MPI::SUM);
1997 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binL[i],
1998 >                                3, MPI::REALTYPE, MPI::SUM);
1999 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binI[i],
2000 >                                9, MPI::REALTYPE, MPI::SUM);
2001 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[i],
2002 >                                1, MPI::REALTYPE, MPI::SUM);
2003 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[i],
2004 >                                1, MPI::INT, MPI::SUM);
2005 >      //MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmega[i],
2006 >      //                          1, MPI::REALTYPE, MPI::SUM);
2007 >    }
2008 >    
2009   #endif
2010  
2011 <    Vector3d vel;
1976 <    Vector3d aVel;
2011 >    Vector3d omega;
2012      RealType den;
2013      RealType temp;
2014      RealType z;
# Line 1990 | Line 2025 | namespace OpenMD {
2025          den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
2026            / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
2027        }
2028 <      vel.x() = binPx[i] / binMass[i];
1994 <      vel.y() = binPy[i] / binMass[i];
1995 <      vel.z() = binPz[i] / binMass[i];
1996 <      aVel.x() = binOmegax[i] / binCount[i];
1997 <      aVel.y() = binOmegay[i] / binCount[i];
1998 <      aVel.z() = binOmegaz[i] / binCount[i];
2028 >      vel = binP[i] / binMass[i];
2029  
2030 +      omega = binI[i].inverse() * binL[i];
2031 +
2032 +      // omega = binOmega[i] / binCount[i];
2033 +
2034        if (binCount[i] > 0) {
2035          // only add values if there are things to add
2036          temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
# Line 2018 | Line 2052 | namespace OpenMD {
2052                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2053                break;
2054              case ANGULARVELOCITY:  
2055 <              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2055 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega);
2056                break;
2057              case DENSITY:
2058                dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);

Diff Legend

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