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 1939 by gezelter, Thu Oct 31 15:32:17 2013 UTC vs.
Revision 1940 by gezelter, Fri Nov 1 19:31:41 2013 UTC

# Line 424 | Line 424 | namespace OpenMD {
424        OutputData angularVelocity;
425        angularVelocity.units = "angstroms^2/fs";
426        angularVelocity.title =  "AngularVelocity";  
427 <      angularVelocity.dataType = "Vector3d";
427 >      angularVelocity.dataType = "RealType";
428        angularVelocity.accumulator.reserve(nBins_);
429        for (int i = 0; i < nBins_; i++)
430 <        angularVelocity.accumulator.push_back( new VectorAccumulator() );
430 >        angularVelocity.accumulator.push_back( new Accumulator() );
431        data_[ANGULARVELOCITY] = angularVelocity;
432        outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
433  
# Line 1874 | 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);
# Line 1884 | Line 1887 | namespace OpenMD {
1887      vector<RealType> binPx(nBins_, 0.0);
1888      vector<RealType> binPy(nBins_, 0.0);
1889      vector<RealType> binPz(nBins_, 0.0);
1890 <    vector<RealType> binOmegax(nBins_, 0.0);
1888 <    vector<RealType> binOmegay(nBins_, 0.0);
1889 <    vector<RealType> binOmegaz(nBins_, 0.0);
1890 >    vector<RealType> binOmega(nBins_, 0.0);
1891      vector<RealType> binKE(nBins_, 0.0);
1892      vector<int> binDOF(nBins_, 0);
1893      vector<int> binCount(nBins_, 0);
# Line 1926 | Line 1927 | namespace OpenMD {
1927          binNo = int(rPos.length() / binWidth_);
1928        }
1929  
1930 +
1931        RealType mass = sd->getMass();
1932        Vector3d vel = sd->getVel();
1933        Vector3d rPos = sd->getPos() - coordinateOrigin_;
1934 <      Vector3d aVel = cross(rPos, vel);
1934 >      // Project the relative position onto a plane perpendicular to
1935 >      // the angularMomentumFluxVector:
1936 >      Vector3d rProj = rPos - dot(rPos, u) * u;
1937 >      // Project the velocity onto a plane perpendicular to the
1938 >      // angularMomentumFluxVector:
1939 >      Vector3d vProj = vel  - dot(vel, u) * u;
1940 >      // Compute angular velocity vector (should be nearly parallel to
1941 >      // angularMomentumFluxVector
1942 >      Vector3d aVel = cross(rProj, vProj);
1943        
1944        if (binNo >= 0 && binNo < nBins_)  {
1945          binCount[binNo]++;
# Line 1937 | Line 1947 | namespace OpenMD {
1947          binPx[binNo] += mass*vel.x();
1948          binPy[binNo] += mass*vel.y();
1949          binPz[binNo] += mass*vel.z();
1950 <        binOmegax[binNo] += aVel.x();
1941 <        binOmegay[binNo] += aVel.y();
1942 <        binOmegaz[binNo] += aVel.z();
1950 >        binOmega[binNo] += dot(aVel, u);
1951          binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
1952          binDOF[binNo] += 3;
1953          
# Line 1974 | Line 1982 | namespace OpenMD {
1982                                nBins_, MPI::REALTYPE, MPI::SUM);
1983      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
1984                                nBins_, MPI::REALTYPE, MPI::SUM);
1985 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
1978 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1979 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
1980 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1981 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
1985 >    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmega[0],
1986                                nBins_, MPI::REALTYPE, MPI::SUM);
1987      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
1988                                nBins_, MPI::REALTYPE, MPI::SUM);
# Line 1987 | Line 1991 | namespace OpenMD {
1991   #endif
1992  
1993      Vector3d vel;
1994 <    Vector3d aVel;
1994 >    RealType omega;
1995      RealType den;
1996      RealType temp;
1997      RealType z;
# Line 2007 | Line 2011 | namespace OpenMD {
2011        vel.x() = binPx[i] / binMass[i];
2012        vel.y() = binPy[i] / binMass[i];
2013        vel.z() = binPz[i] / binMass[i];
2014 <      aVel.x() = binOmegax[i] / binCount[i];
2011 <      aVel.y() = binOmegay[i] / binCount[i];
2012 <      aVel.z() = binOmegaz[i] / binCount[i];
2014 >      omega = binOmega[i] / binCount[i];
2015  
2016        if (binCount[i] > 0) {
2017          // only add values if there are things to add
# Line 2032 | Line 2034 | namespace OpenMD {
2034                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2035                break;
2036              case ANGULARVELOCITY:  
2037 <              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2037 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(omega);
2038                break;
2039              case DENSITY:
2040                dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);

Diff Legend

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