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 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC vs.
Revision 2026 by gezelter, Wed Oct 22 12:23:59 2014 UTC

# Line 291 | Line 291 | namespace OpenMD {
291        kineticFlux_ = 0.0;
292      }
293      if (hasMomentumFluxVector) {
294 <      momentumFluxVector_ = rnemdParams->getMomentumFluxVector();
294 >      std::vector<RealType> mf = rnemdParams->getMomentumFluxVector();
295 >      if (mf.size() != 3) {
296 >          sprintf(painCave.errMsg,
297 >                  "RNEMD: Incorrect number of parameters specified for momentumFluxVector.\n"
298 >                  "\tthere should be 3 parameters, but %lu were specified.\n",
299 >                  mf.size());
300 >          painCave.isFatal = 1;
301 >          simError();      
302 >      }
303 >      momentumFluxVector_.x() = mf[0];
304 >      momentumFluxVector_.y() = mf[1];
305 >      momentumFluxVector_.z() = mf[2];
306      } else {
307        momentumFluxVector_ = V3Zero;
308        if (hasMomentumFlux) {
# Line 317 | Line 328 | namespace OpenMD {
328          }
329        }
330        if (hasAngularMomentumFluxVector) {
331 <        angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
331 >        std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
332 >        if (amf.size() != 3) {
333 >          sprintf(painCave.errMsg,
334 >                  "RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
335 >                  "\tthere should be 3 parameters, but %lu were specified.\n",
336 >                  amf.size());
337 >          painCave.isFatal = 1;
338 >          simError();      
339 >        }
340 >        angularMomentumFluxVector_.x() = amf[0];
341 >        angularMomentumFluxVector_.y() = amf[1];
342 >        angularMomentumFluxVector_.z() = amf[2];
343        } else {
344          angularMomentumFluxVector_ = V3Zero;
345          if (hasAngularMomentumFlux) {
# Line 348 | Line 370 | namespace OpenMD {
370        }
371  
372        if (hasCoordinateOrigin) {
373 <        coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
373 >        std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
374 >        if (co.size() != 3) {
375 >          sprintf(painCave.errMsg,
376 >                  "RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
377 >                  "\tthere should be 3 parameters, but %lu were specified.\n",
378 >                  co.size());
379 >          painCave.isFatal = 1;
380 >          simError();      
381 >        }
382 >        coordinateOrigin_.x() = co[0];
383 >        coordinateOrigin_.y() = co[1];
384 >        coordinateOrigin_.z() = co[2];
385        } else {
386          coordinateOrigin_ = V3Zero;
387        }
# Line 623 | Line 656 | namespace OpenMD {
656      StuntDouble* sd;
657  
658      RealType min_val;
659 <    bool min_found = false;  
659 >    int min_found = 0;  
660      StuntDouble* min_sd;
661  
662      RealType max_val;
663 <    bool max_found = false;
663 >    int max_found = 0;
664      StuntDouble* max_sd;
665  
666      for (sd = seleManA_.beginSelected(selei); sd != NULL;
# Line 682 | Line 715 | namespace OpenMD {
715        if (!max_found) {
716          max_val = value;
717          max_sd = sd;
718 <        max_found = true;
718 >        max_found = 1;
719        } else {
720          if (max_val < value) {
721            max_val = value;
# Line 744 | Line 777 | namespace OpenMD {
777        if (!min_found) {
778          min_val = value;
779          min_sd = sd;
780 <        min_found = true;
780 >        min_found = 1;
781        } else {
782          if (min_val > value) {
783            min_val = value;
# Line 754 | Line 787 | namespace OpenMD {
787      }
788      
789   #ifdef IS_MPI    
790 <    int worldRank = MPI::COMM_WORLD.Get_rank();
791 <    
792 <    bool my_min_found = min_found;
793 <    bool my_max_found = max_found;
790 >    int worldRank;
791 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
792 >        
793 >    int my_min_found = min_found;
794 >    int my_max_found = max_found;
795  
796      // Even if we didn't find a minimum, did someone else?
797 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
797 >    MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
798 >                  MPI_COMM_WORLD);
799      // Even if we didn't find a maximum, did someone else?
800 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
800 >    MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
801 >                  MPI_COMM_WORLD);
802   #endif
803  
804      if (max_found && min_found) {
# Line 781 | Line 817 | namespace OpenMD {
817        min_vals.rank = worldRank;    
818        
819        // Who had the minimum?
820 <      MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
821 <                                1, MPI::REALTYPE_INT, MPI::MINLOC);
820 >      MPI_Allreduce(&min_vals, &min_vals,
821 >                    1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD);
822        min_val = min_vals.val;
823        
824        if (my_max_found) {
# Line 793 | Line 829 | namespace OpenMD {
829        max_vals.rank = worldRank;    
830        
831        // Who had the maximum?
832 <      MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
833 <                                1, MPI::REALTYPE_INT, MPI::MAXLOC);
832 >      MPI_Allreduce(&max_vals, &max_vals,
833 >                    1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
834        max_val = max_vals.val;
835   #endif
836        
# Line 854 | Line 890 | namespace OpenMD {
890            
891            Vector3d min_vel;
892            Vector3d max_vel = max_sd->getVel();
893 <          MPI::Status status;
893 >          MPI_Status status;
894  
895            // point-to-point swap of the velocity vector
896 <          MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
897 <                                   min_vals.rank, 0,
898 <                                   min_vel.getArrayPointer(), 3, MPI::REALTYPE,
899 <                                   min_vals.rank, 0, status);
896 >          MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
897 >                       min_vals.rank, 0,
898 >                       min_vel.getArrayPointer(), 3, MPI_REALTYPE,
899 >                       min_vals.rank, 0, MPI_COMM_WORLD, &status);
900            
901            switch(rnemdFluxType_) {
902            case rnemdKE :
# Line 871 | Line 907 | namespace OpenMD {
907                Vector3d max_angMom = max_sd->getJ();
908                
909                // point-to-point swap of the angular momentum vector
910 <              MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
911 <                                       MPI::REALTYPE, min_vals.rank, 1,
912 <                                       min_angMom.getArrayPointer(), 3,
913 <                                       MPI::REALTYPE, min_vals.rank, 1,
914 <                                       status);
910 >              MPI_Sendrecv(max_angMom.getArrayPointer(), 3,
911 >                           MPI_REALTYPE, min_vals.rank, 1,
912 >                           min_angMom.getArrayPointer(), 3,
913 >                           MPI_REALTYPE, min_vals.rank, 1,
914 >                           MPI_COMM_WORLD, &status);
915                
916                max_sd->setJ(min_angMom);
917              }
# Line 900 | Line 936 | namespace OpenMD {
936            
937            Vector3d max_vel;
938            Vector3d min_vel = min_sd->getVel();
939 <          MPI::Status status;
939 >          MPI_Status status;
940            
941            // point-to-point swap of the velocity vector
942 <          MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
943 <                                   max_vals.rank, 0,
944 <                                   max_vel.getArrayPointer(), 3, MPI::REALTYPE,
945 <                                   max_vals.rank, 0, status);
942 >          MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
943 >                       max_vals.rank, 0,
944 >                       max_vel.getArrayPointer(), 3, MPI_REALTYPE,
945 >                       max_vals.rank, 0, MPI_COMM_WORLD, &status);
946            
947            switch(rnemdFluxType_) {
948            case rnemdKE :
# Line 917 | Line 953 | namespace OpenMD {
953                Vector3d max_angMom;
954                
955                // point-to-point swap of the angular momentum vector
956 <              MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
957 <                                       MPI::REALTYPE, max_vals.rank, 1,
958 <                                       max_angMom.getArrayPointer(), 3,
959 <                                       MPI::REALTYPE, max_vals.rank, 1,
960 <                                       status);
956 >              MPI_Sendrecv(min_angMom.getArrayPointer(), 3,
957 >                           MPI_REALTYPE, max_vals.rank, 1,
958 >                           max_angMom.getArrayPointer(), 3,
959 >                           MPI_REALTYPE, max_vals.rank, 1,
960 >                           MPI_COMM_WORLD, &status);
961                
962                min_sd->setJ(max_angMom);
963              }
# Line 1090 | Line 1126 | namespace OpenMD {
1126      Kcw *= 0.5;
1127  
1128   #ifdef IS_MPI
1129 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
1130 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
1131 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
1132 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
1133 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
1134 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
1129 >    MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1130 >    MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1131 >    MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1132 >    MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1133 >    MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1134 >    MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1135  
1136 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
1137 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
1138 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
1139 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM);
1136 >    MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1137 >    MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1138 >    MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1139 >    MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1140  
1141 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
1142 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
1143 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
1144 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM);
1141 >    MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1142 >    MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1143 >    MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1144 >    MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1145   #endif
1146  
1147      //solve coldBin coeff's first
# Line 1582 | Line 1618 | namespace OpenMD {
1618      Kc *= 0.5;
1619      
1620   #ifdef IS_MPI
1621 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1622 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
1623 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM);
1624 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM);
1625 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM);
1626 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM);
1627 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM);
1628 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM);
1629 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9,
1630 <                              MPI::REALTYPE, MPI::SUM);
1631 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9,
1632 <                              MPI::REALTYPE, MPI::SUM);
1621 >    MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,
1622 >                  MPI_COMM_WORLD);
1623 >    MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,
1624 >                  MPI_COMM_WORLD);
1625 >    MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,
1626 >                  MPI_COMM_WORLD);
1627 >    MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,
1628 >                  MPI_COMM_WORLD);
1629 >    MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1630 >    MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1631 >    MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1632 >    MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1633 >    MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,
1634 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1635 >    MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,
1636 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1637   #endif
1638      
1639  
# Line 1748 | Line 1788 | namespace OpenMD {
1788   #if defined(HAVE_QHULL)
1789          ConvexHull* surfaceMeshA = new ConvexHull();
1790          surfaceMeshA->computeHull(aSites);
1791 +        cerr << "flag1\n";
1792          areaA = surfaceMeshA->getArea();
1793 +        cerr << "Flag2 " << areaA << "\n";
1794          delete surfaceMeshA;
1795   #else
1796          sprintf( painCave.errMsg,
# Line 1874 | Line 1916 | namespace OpenMD {
1916      RealType area = getDividingArea();
1917      areaAccumulator_->add(area);
1918      Mat3x3d hmat = currentSnap_->getHmat();
1919 +    Vector3d u = angularMomentumFluxVector_;
1920 +    u.normalize();
1921 +
1922      seleMan_.setSelectionSet(evaluator_.evaluate());
1923  
1924      int selei(0);
1925      StuntDouble* sd;
1926      int binNo;
1927 +    RealType mass;
1928 +    Vector3d vel;
1929 +    Vector3d rPos;
1930 +    RealType KE;
1931 +    Vector3d L;
1932 +    Mat3x3d I;
1933 +    RealType r2;
1934  
1935      vector<RealType> binMass(nBins_, 0.0);
1936 <    vector<RealType> binPx(nBins_, 0.0);
1937 <    vector<RealType> binPy(nBins_, 0.0);
1938 <    vector<RealType> binPz(nBins_, 0.0);
1939 <    vector<RealType> binOmegax(nBins_, 0.0);
1888 <    vector<RealType> binOmegay(nBins_, 0.0);
1889 <    vector<RealType> binOmegaz(nBins_, 0.0);
1936 >    vector<Vector3d> binP(nBins_, V3Zero);
1937 >    vector<RealType> binOmega(nBins_, 0.0);
1938 >    vector<Vector3d> binL(nBins_, V3Zero);
1939 >    vector<Mat3x3d>  binI(nBins_);
1940      vector<RealType> binKE(nBins_, 0.0);
1941      vector<int> binDOF(nBins_, 0);
1942      vector<int> binCount(nBins_, 0);
# Line 1926 | Line 1976 | namespace OpenMD {
1976          binNo = int(rPos.length() / binWidth_);
1977        }
1978  
1979 <      RealType mass = sd->getMass();
1980 <      Vector3d vel = sd->getVel();
1981 <      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1982 <      Vector3d aVel = cross(rPos, vel);
1983 <      
1979 >      mass = sd->getMass();
1980 >      vel = sd->getVel();
1981 >      rPos = sd->getPos() - coordinateOrigin_;
1982 >      KE = 0.5 * mass * vel.lengthSquare();
1983 >      L = mass * cross(rPos, vel);
1984 >      I = outProduct(rPos, rPos) * mass;
1985 >      r2 = rPos.lengthSquare();
1986 >      I(0, 0) += mass * r2;
1987 >      I(1, 1) += mass * r2;
1988 >      I(2, 2) += mass * r2;
1989 >
1990 >      // Project the relative position onto a plane perpendicular to
1991 >      // the angularMomentumFluxVector:
1992 >      // Vector3d rProj = rPos - dot(rPos, u) * u;
1993 >      // Project the velocity onto a plane perpendicular to the
1994 >      // angularMomentumFluxVector:
1995 >      // Vector3d vProj = vel  - dot(vel, u) * u;
1996 >      // Compute angular velocity vector (should be nearly parallel to
1997 >      // angularMomentumFluxVector
1998 >      // Vector3d aVel = cross(rProj, vProj);
1999 >
2000        if (binNo >= 0 && binNo < nBins_)  {
2001          binCount[binNo]++;
2002          binMass[binNo] += mass;
2003 <        binPx[binNo] += mass*vel.x();
2004 <        binPy[binNo] += mass*vel.y();
2005 <        binPz[binNo] += mass*vel.z();
2006 <        binOmegax[binNo] += aVel.x();
1941 <        binOmegay[binNo] += aVel.y();
1942 <        binOmegaz[binNo] += aVel.z();
1943 <        binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
2003 >        binP[binNo] += mass*vel;
2004 >        binKE[binNo] += KE;
2005 >        binI[binNo] += I;
2006 >        binL[binNo] += L;
2007          binDOF[binNo] += 3;
2008          
2009          if (sd->isDirectional()) {
2010            Vector3d angMom = sd->getJ();
2011 <          Mat3x3d I = sd->getI();
2011 >          Mat3x3d Ia = sd->getI();
2012            if (sd->isLinear()) {
2013              int i = sd->linearAxis();
2014              int j = (i + 1) % 3;
2015              int k = (i + 2) % 3;
2016 <            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
2017 <                                   angMom[k] * angMom[k] / I(k, k));
2016 >            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
2017 >                                   angMom[k] * angMom[k] / Ia(k, k));
2018              binDOF[binNo] += 2;
2019            } else {
2020 <            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
2021 <                                   angMom[1] * angMom[1] / I(1, 1) +
2022 <                                   angMom[2] * angMom[2] / I(2, 2));
2020 >            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
2021 >                                   angMom[1] * angMom[1] / Ia(1, 1) +
2022 >                                   angMom[2] * angMom[2] / Ia(2, 2));
2023              binDOF[binNo] += 3;
2024            }
2025          }
# Line 1964 | Line 2027 | namespace OpenMD {
2027      }
2028      
2029   #ifdef IS_MPI
2030 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
2031 <                              nBins_, MPI::INT, MPI::SUM);
2032 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0],
2033 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2034 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0],
2035 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2036 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0],
2037 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2038 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
2039 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2040 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
2041 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2042 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
2043 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2044 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
2045 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2046 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
2047 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2048 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0],
2049 <                              nBins_, MPI::INT, MPI::SUM);
2030 >
2031 >    for (int i = 0; i < nBins_; i++) {
2032 >      
2033 >      MPI_Allreduce(MPI_IN_PLACE, &binCount[i],
2034 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2035 >      MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2036 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2037 >      MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(),
2038 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2039 >      MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(),
2040 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2041 >      MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(),
2042 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2043 >      MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2044 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2045 >      MPI_Allreduce(MPI_IN_PLACE, &binDOF[i],
2046 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2047 >      //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i],
2048 >      //                          1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2049 >    }
2050 >    
2051   #endif
2052  
2053 <    Vector3d vel;
1990 <    Vector3d aVel;
2053 >    Vector3d omega;
2054      RealType den;
2055      RealType temp;
2056      RealType z;
# Line 2004 | Line 2067 | namespace OpenMD {
2067          den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
2068            / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
2069        }
2070 <      vel.x() = binPx[i] / binMass[i];
2008 <      vel.y() = binPy[i] / binMass[i];
2009 <      vel.z() = binPz[i] / binMass[i];
2010 <      aVel.x() = binOmegax[i] / binCount[i];
2011 <      aVel.y() = binOmegay[i] / binCount[i];
2012 <      aVel.z() = binOmegaz[i] / binCount[i];
2070 >      vel = binP[i] / binMass[i];
2071  
2072 +      omega = binI[i].inverse() * binL[i];
2073 +
2074 +      // omega = binOmega[i] / binCount[i];
2075 +
2076        if (binCount[i] > 0) {
2077          // only add values if there are things to add
2078          temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
# Line 2032 | Line 2094 | namespace OpenMD {
2094                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2095                break;
2096              case ANGULARVELOCITY:  
2097 <              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2097 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega);
2098                break;
2099              case DENSITY:
2100                dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
# Line 2079 | Line 2141 | namespace OpenMD {
2141      
2142   #ifdef IS_MPI
2143      // If we're the root node, should we print out the results
2144 <    int worldRank = MPI::COMM_WORLD.Get_rank();
2144 >    int worldRank;
2145 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
2146 >
2147      if (worldRank == 0) {
2148   #endif
2149        rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc );
# Line 2210 | Line 2274 | namespace OpenMD {
2274        }        
2275  
2276        rnemdFile_ << "#######################################################\n";
2277 <      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
2277 >      rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
2278        rnemdFile_ << "#######################################################\n";
2279  
2280  
# Line 2219 | Line 2283 | namespace OpenMD {
2283          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2284            if (outputMask_[i]) {
2285              if (data_[i].dataType == "RealType")
2286 <              writeRealStdDev(i,j);
2286 >              writeRealErrorBars(i,j);
2287              else if (data_[i].dataType == "Vector3d")
2288 <              writeVectorStdDev(i,j);
2288 >              writeVectorErrorBars(i,j);
2289              else {
2290                sprintf( painCave.errMsg,
2291                         "RNEMD found an unknown data type for: %s ",
# Line 2292 | Line 2356 | namespace OpenMD {
2356      }
2357    }  
2358  
2359 <  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2359 >  void RNEMD::writeRealErrorBars(int index, unsigned int bin) {
2360      if (!doRNEMD_) return;
2361      assert(index >=0 && index < ENDINDEX);
2362      assert(int(bin) < nBins_);
# Line 2302 | Line 2366 | namespace OpenMD {
2366      count = data_[index].accumulator[bin]->count();
2367      if (count == 0) return;
2368      
2369 <    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2369 >    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2370      
2371      if (! isinf(s) && ! isnan(s)) {
2372        rnemdFile_ << "\t" << s;
# Line 2315 | Line 2379 | namespace OpenMD {
2379      }    
2380    }
2381    
2382 <  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2382 >  void RNEMD::writeVectorErrorBars(int index, unsigned int bin) {
2383      if (!doRNEMD_) return;
2384      assert(index >=0 && index < ENDINDEX);
2385      assert(int(bin) < nBins_);
# Line 2325 | Line 2389 | namespace OpenMD {
2389      count = data_[index].accumulator[bin]->count();
2390      if (count == 0) return;
2391  
2392 <    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2392 >    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2393      if (isinf(s[0]) || isnan(s[0]) ||
2394          isinf(s[1]) || isnan(s[1]) ||
2395          isinf(s[2]) || isnan(s[2]) ) {      

Diff Legend

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