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 1969 by gezelter, Wed Feb 26 14:14:50 2014 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 752 | Line 754 | namespace OpenMD {
754      }
755      
756   #ifdef IS_MPI    
757 <    int worldRank = MPI::COMM_WORLD.Get_rank();
758 <    
759 <    bool my_min_found = min_found;
760 <    bool my_max_found = max_found;
757 >    int worldRank;
758 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
759 >        
760 >    int my_min_found = min_found;
761 >    int my_max_found = max_found;
762  
763      // Even if we didn't find a minimum, did someone else?
764 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
764 >    MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
765 >                  MPI_COMM_WORLD);
766      // Even if we didn't find a maximum, did someone else?
767 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
767 >    MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
768 >                  MPI_COMM_WORLD);
769   #endif
770  
771      if (max_found && min_found) {
# Line 779 | Line 784 | namespace OpenMD {
784        min_vals.rank = worldRank;    
785        
786        // Who had the minimum?
787 <      MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
788 <                                1, MPI::REALTYPE_INT, MPI::MINLOC);
787 >      MPI_Allreduce(&min_vals, &min_vals,
788 >                    1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD);
789        min_val = min_vals.val;
790        
791        if (my_max_found) {
# Line 791 | Line 796 | namespace OpenMD {
796        max_vals.rank = worldRank;    
797        
798        // Who had the maximum?
799 <      MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
800 <                                1, MPI::REALTYPE_INT, MPI::MAXLOC);
799 >      MPI_Allreduce(&max_vals, &max_vals,
800 >                    1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
801        max_val = max_vals.val;
802   #endif
803        
# Line 852 | Line 857 | namespace OpenMD {
857            
858            Vector3d min_vel;
859            Vector3d max_vel = max_sd->getVel();
860 <          MPI::Status status;
860 >          MPI_Status* status;
861  
862            // point-to-point swap of the velocity vector
863 <          MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
864 <                                   min_vals.rank, 0,
865 <                                   min_vel.getArrayPointer(), 3, MPI::REALTYPE,
866 <                                   min_vals.rank, 0, status);
863 >          MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
864 >                       min_vals.rank, 0,
865 >                       min_vel.getArrayPointer(), 3, MPI_REALTYPE,
866 >                       min_vals.rank, 0, MPI_COMM_WORLD, status);
867            
868            switch(rnemdFluxType_) {
869            case rnemdKE :
# Line 869 | Line 874 | namespace OpenMD {
874                Vector3d max_angMom = max_sd->getJ();
875                
876                // point-to-point swap of the angular momentum vector
877 <              MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
878 <                                       MPI::REALTYPE, min_vals.rank, 1,
879 <                                       min_angMom.getArrayPointer(), 3,
880 <                                       MPI::REALTYPE, min_vals.rank, 1,
881 <                                       status);
877 >              MPI_Sendrecv(max_angMom.getArrayPointer(), 3,
878 >                           MPI_REALTYPE, min_vals.rank, 1,
879 >                           min_angMom.getArrayPointer(), 3,
880 >                           MPI_REALTYPE, min_vals.rank, 1,
881 >                           MPI_COMM_WORLD, status);
882                
883                max_sd->setJ(min_angMom);
884              }
# Line 898 | Line 903 | namespace OpenMD {
903            
904            Vector3d max_vel;
905            Vector3d min_vel = min_sd->getVel();
906 <          MPI::Status status;
906 >          MPI_Status* status;
907            
908            // point-to-point swap of the velocity vector
909 <          MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
910 <                                   max_vals.rank, 0,
911 <                                   max_vel.getArrayPointer(), 3, MPI::REALTYPE,
912 <                                   max_vals.rank, 0, status);
909 >          MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
910 >                       max_vals.rank, 0,
911 >                       max_vel.getArrayPointer(), 3, MPI_REALTYPE,
912 >                       max_vals.rank, 0, MPI_COMM_WORLD, status);
913            
914            switch(rnemdFluxType_) {
915            case rnemdKE :
# Line 915 | Line 920 | namespace OpenMD {
920                Vector3d max_angMom;
921                
922                // point-to-point swap of the angular momentum vector
923 <              MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
924 <                                       MPI::REALTYPE, max_vals.rank, 1,
925 <                                       max_angMom.getArrayPointer(), 3,
926 <                                       MPI::REALTYPE, max_vals.rank, 1,
927 <                                       status);
923 >              MPI_Sendrecv(min_angMom.getArrayPointer(), 3,
924 >                           MPI_REALTYPE, max_vals.rank, 1,
925 >                           max_angMom.getArrayPointer(), 3,
926 >                           MPI_REALTYPE, max_vals.rank, 1,
927 >                           MPI_COMM_WORLD, status);
928                
929                min_sd->setJ(max_angMom);
930              }
# Line 1088 | Line 1093 | namespace OpenMD {
1093      Kcw *= 0.5;
1094  
1095   #ifdef IS_MPI
1096 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
1097 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
1098 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
1099 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
1100 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
1101 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
1096 >    MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1097 >    MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1098 >    MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1099 >    MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1100 >    MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1101 >    MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1102  
1103 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
1104 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
1105 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
1106 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM);
1103 >    MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1104 >    MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1105 >    MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1106 >    MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1107  
1108 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
1109 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
1110 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
1111 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM);
1108 >    MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1109 >    MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1110 >    MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1111 >    MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1112   #endif
1113  
1114      //solve coldBin coeff's first
# Line 1580 | Line 1585 | namespace OpenMD {
1585      Kc *= 0.5;
1586      
1587   #ifdef IS_MPI
1588 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1589 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
1590 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM);
1591 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM);
1592 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM);
1593 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM);
1594 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM);
1595 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM);
1596 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9,
1597 <                              MPI::REALTYPE, MPI::SUM);
1598 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9,
1599 <                              MPI::REALTYPE, MPI::SUM);
1588 >    MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,
1589 >                  MPI_COMM_WORLD);
1590 >    MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,
1591 >                  MPI_COMM_WORLD);
1592 >    MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,
1593 >                  MPI_COMM_WORLD);
1594 >    MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,
1595 >                  MPI_COMM_WORLD);
1596 >    MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1597 >    MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1598 >    MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1599 >    MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1600 >    MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,
1601 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1602 >    MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,
1603 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1604   #endif
1605      
1606  
# Line 1729 | Line 1738 | namespace OpenMD {
1738      Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1739  
1740      if (hasSelectionA_) {
1741 <      int isd;
1742 <      StuntDouble* sd;
1743 <      vector<StuntDouble*> aSites;
1744 <      seleManA_.setSelectionSet(evaluatorA_.evaluate());
1745 <      for (sd = seleManA_.beginSelected(isd); sd != NULL;
1746 <           sd = seleManA_.nextSelected(isd)) {
1747 <        aSites.push_back(sd);
1748 <      }
1741 >
1742 >      if (evaluatorA_.hasSurfaceArea())
1743 >        areaA = evaluatorA_.getSurfaceArea();
1744 >      else {
1745 >        
1746 >        cerr << "selection A did not have surface area, recomputing\n";
1747 >        int isd;
1748 >        StuntDouble* sd;
1749 >        vector<StuntDouble*> aSites;
1750 >        seleManA_.setSelectionSet(evaluatorA_.evaluate());
1751 >        for (sd = seleManA_.beginSelected(isd); sd != NULL;
1752 >             sd = seleManA_.nextSelected(isd)) {
1753 >          aSites.push_back(sd);
1754 >        }
1755   #if defined(HAVE_QHULL)
1756 <      ConvexHull* surfaceMeshA = new ConvexHull();
1757 <      surfaceMeshA->computeHull(aSites);
1758 <      areaA = surfaceMeshA->getArea();
1759 <      delete surfaceMeshA;
1756 >        ConvexHull* surfaceMeshA = new ConvexHull();
1757 >        surfaceMeshA->computeHull(aSites);
1758 >        areaA = surfaceMeshA->getArea();
1759 >        delete surfaceMeshA;
1760   #else
1761 <      sprintf( painCave.errMsg,
1761 >        sprintf( painCave.errMsg,
1762                 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1763 <               "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1764 <      painCave.severity = OPENMD_ERROR;
1765 <      painCave.isFatal = 1;
1766 <      simError();
1763 >                 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1764 >        painCave.severity = OPENMD_ERROR;
1765 >        painCave.isFatal = 1;
1766 >        simError();
1767   #endif
1768 +      }
1769  
1770      } else {
1771        if (usePeriodicBoundaryConditions_) {
# Line 1765 | Line 1781 | namespace OpenMD {
1781      }
1782  
1783      if (hasSelectionB_) {
1784 <      int isd;
1785 <      StuntDouble* sd;
1786 <      vector<StuntDouble*> bSites;
1787 <      seleManB_.setSelectionSet(evaluatorB_.evaluate());
1772 <      for (sd = seleManB_.beginSelected(isd); sd != NULL;
1773 <           sd = seleManB_.nextSelected(isd)) {
1774 <        bSites.push_back(sd);
1775 <      }
1784 >      if (evaluatorB_.hasSurfaceArea())
1785 >        areaB = evaluatorB_.getSurfaceArea();
1786 >      else {
1787 >        cerr << "selection B did not have surface area, recomputing\n";
1788  
1789 +        int isd;
1790 +        StuntDouble* sd;
1791 +        vector<StuntDouble*> bSites;
1792 +        seleManB_.setSelectionSet(evaluatorB_.evaluate());
1793 +        for (sd = seleManB_.beginSelected(isd); sd != NULL;
1794 +             sd = seleManB_.nextSelected(isd)) {
1795 +          bSites.push_back(sd);
1796 +        }
1797 +        
1798   #if defined(HAVE_QHULL)
1799 <      ConvexHull* surfaceMeshB = new ConvexHull();    
1800 <      surfaceMeshB->computeHull(bSites);
1801 <      areaB = surfaceMeshB->getArea();
1802 <      delete surfaceMeshB;
1799 >        ConvexHull* surfaceMeshB = new ConvexHull();    
1800 >        surfaceMeshB->computeHull(bSites);
1801 >        areaB = surfaceMeshB->getArea();
1802 >        delete surfaceMeshB;
1803   #else
1804 <      sprintf( painCave.errMsg,
1805 <               "RNEMD::getDividingArea : Hull calculation is not possible\n"
1806 <               "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1807 <      painCave.severity = OPENMD_ERROR;
1808 <      painCave.isFatal = 1;
1809 <      simError();
1804 >        sprintf( painCave.errMsg,
1805 >                 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1806 >                 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1807 >        painCave.severity = OPENMD_ERROR;
1808 >        painCave.isFatal = 1;
1809 >        simError();
1810   #endif
1811 <
1812 <
1811 >      }
1812 >      
1813      } else {
1814        if (usePeriodicBoundaryConditions_) {
1815          // in periodic boundaries, the surface area is twice the x-y
# Line 1800 | Line 1821 | namespace OpenMD {
1821          areaB = 4.0 * M_PI * pow(sphereBRadius_, 2);
1822        }
1823      }
1824 <    
1824 >      
1825      dividingArea_ = min(areaA, areaB);
1826      hasDividingArea_ = true;
1827      return dividingArea_;
# Line 1860 | Line 1881 | namespace OpenMD {
1881      RealType area = getDividingArea();
1882      areaAccumulator_->add(area);
1883      Mat3x3d hmat = currentSnap_->getHmat();
1884 +    Vector3d u = angularMomentumFluxVector_;
1885 +    u.normalize();
1886 +
1887      seleMan_.setSelectionSet(evaluator_.evaluate());
1888  
1889      int selei(0);
1890      StuntDouble* sd;
1891      int binNo;
1892 +    RealType mass;
1893 +    Vector3d vel;
1894 +    Vector3d rPos;
1895 +    RealType KE;
1896 +    Vector3d L;
1897 +    Mat3x3d I;
1898 +    RealType r2;
1899  
1900      vector<RealType> binMass(nBins_, 0.0);
1901 <    vector<RealType> binPx(nBins_, 0.0);
1902 <    vector<RealType> binPy(nBins_, 0.0);
1903 <    vector<RealType> binPz(nBins_, 0.0);
1904 <    vector<RealType> binOmegax(nBins_, 0.0);
1874 <    vector<RealType> binOmegay(nBins_, 0.0);
1875 <    vector<RealType> binOmegaz(nBins_, 0.0);
1901 >    vector<Vector3d> binP(nBins_, V3Zero);
1902 >    vector<RealType> binOmega(nBins_, 0.0);
1903 >    vector<Vector3d> binL(nBins_, V3Zero);
1904 >    vector<Mat3x3d>  binI(nBins_);
1905      vector<RealType> binKE(nBins_, 0.0);
1906      vector<int> binDOF(nBins_, 0);
1907      vector<int> binCount(nBins_, 0);
# Line 1912 | Line 1941 | namespace OpenMD {
1941          binNo = int(rPos.length() / binWidth_);
1942        }
1943  
1944 <      RealType mass = sd->getMass();
1945 <      Vector3d vel = sd->getVel();
1946 <      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1947 <      Vector3d aVel = cross(rPos, vel);
1948 <      
1944 >      mass = sd->getMass();
1945 >      vel = sd->getVel();
1946 >      rPos = sd->getPos() - coordinateOrigin_;
1947 >      KE = 0.5 * mass * vel.lengthSquare();
1948 >      L = mass * cross(rPos, vel);
1949 >      I = outProduct(rPos, rPos) * mass;
1950 >      r2 = rPos.lengthSquare();
1951 >      I(0, 0) += mass * r2;
1952 >      I(1, 1) += mass * r2;
1953 >      I(2, 2) += mass * r2;
1954 >
1955 >      // Project the relative position onto a plane perpendicular to
1956 >      // the angularMomentumFluxVector:
1957 >      // Vector3d rProj = rPos - dot(rPos, u) * u;
1958 >      // Project the velocity onto a plane perpendicular to the
1959 >      // angularMomentumFluxVector:
1960 >      // Vector3d vProj = vel  - dot(vel, u) * u;
1961 >      // Compute angular velocity vector (should be nearly parallel to
1962 >      // angularMomentumFluxVector
1963 >      // Vector3d aVel = cross(rProj, vProj);
1964 >
1965        if (binNo >= 0 && binNo < nBins_)  {
1966          binCount[binNo]++;
1967          binMass[binNo] += mass;
1968 <        binPx[binNo] += mass*vel.x();
1969 <        binPy[binNo] += mass*vel.y();
1970 <        binPz[binNo] += mass*vel.z();
1971 <        binOmegax[binNo] += aVel.x();
1927 <        binOmegay[binNo] += aVel.y();
1928 <        binOmegaz[binNo] += aVel.z();
1929 <        binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
1968 >        binP[binNo] += mass*vel;
1969 >        binKE[binNo] += KE;
1970 >        binI[binNo] += I;
1971 >        binL[binNo] += L;
1972          binDOF[binNo] += 3;
1973          
1974          if (sd->isDirectional()) {
1975            Vector3d angMom = sd->getJ();
1976 <          Mat3x3d I = sd->getI();
1976 >          Mat3x3d Ia = sd->getI();
1977            if (sd->isLinear()) {
1978              int i = sd->linearAxis();
1979              int j = (i + 1) % 3;
1980              int k = (i + 2) % 3;
1981 <            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
1982 <                                   angMom[k] * angMom[k] / I(k, k));
1981 >            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
1982 >                                   angMom[k] * angMom[k] / Ia(k, k));
1983              binDOF[binNo] += 2;
1984            } else {
1985 <            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
1986 <                                   angMom[1] * angMom[1] / I(1, 1) +
1987 <                                   angMom[2] * angMom[2] / I(2, 2));
1985 >            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
1986 >                                   angMom[1] * angMom[1] / Ia(1, 1) +
1987 >                                   angMom[2] * angMom[2] / Ia(2, 2));
1988              binDOF[binNo] += 3;
1989            }
1990          }
# Line 1950 | Line 1992 | namespace OpenMD {
1992      }
1993      
1994   #ifdef IS_MPI
1995 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1996 <                              nBins_, MPI::INT, MPI::SUM);
1997 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0],
1998 <                              nBins_, MPI::REALTYPE, MPI::SUM);
1999 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0],
2000 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2001 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0],
2002 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2003 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
2004 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2005 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
2006 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2007 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
2008 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2009 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
2010 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2011 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
2012 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2013 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0],
2014 <                              nBins_, MPI::INT, MPI::SUM);
1995 >
1996 >    for (int i = 0; i < nBins_; i++) {
1997 >      
1998 >      MPI_Allreduce(MPI_IN_PLACE, &binCount[i],
1999 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2000 >      MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2001 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2002 >      MPI_Allreduce(MPI_IN_PLACE, &binP[i],
2003 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2004 >      MPI_Allreduce(MPI_IN_PLACE, &binL[i],
2005 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2006 >      MPI_Allreduce(MPI_IN_PLACE, &binI[i],
2007 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2008 >      MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2009 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2010 >      MPI_Allreduce(MPI_IN_PLACE, &binDOF[i],
2011 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2012 >      //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i],
2013 >      //                          1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2014 >    }
2015 >    
2016   #endif
2017  
2018 <    Vector3d vel;
1976 <    Vector3d aVel;
2018 >    Vector3d omega;
2019      RealType den;
2020      RealType temp;
2021      RealType z;
# Line 1990 | Line 2032 | namespace OpenMD {
2032          den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
2033            / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
2034        }
2035 <      vel.x() = binPx[i] / binMass[i];
2036 <      vel.y() = binPy[i] / binMass[i];
2037 <      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];
2035 >      vel = binP[i] / binMass[i];
2036 >
2037 >      omega = binI[i].inverse() * binL[i];
2038  
2039 +      // omega = binOmega[i] / binCount[i];
2040 +
2041        if (binCount[i] > 0) {
2042          // only add values if there are things to add
2043          temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
# Line 2018 | Line 2059 | namespace OpenMD {
2059                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2060                break;
2061              case ANGULARVELOCITY:  
2062 <              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2062 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega);
2063                break;
2064              case DENSITY:
2065                dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
# Line 2065 | Line 2106 | namespace OpenMD {
2106      
2107   #ifdef IS_MPI
2108      // If we're the root node, should we print out the results
2109 <    int worldRank = MPI::COMM_WORLD.Get_rank();
2109 >    int worldRank;
2110 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
2111 >
2112      if (worldRank == 0) {
2113   #endif
2114        rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc );

Diff Legend

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