| 71 | 
  | 
  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),  | 
| 72 | 
  | 
                                evaluatorA_(info), seleManA_(info),  | 
| 73 | 
  | 
                                commonA_(info), evaluatorB_(info),  | 
| 74 | 
< | 
                                seleManB_(info), commonB_(info), | 
| 74 | 
> | 
                                seleManB_(info), commonB_(info),  | 
| 75 | 
> | 
                                hasData_(false), hasDividingArea_(false), | 
| 76 | 
  | 
                                usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { | 
| 77 | 
  | 
 | 
| 78 | 
  | 
    trialCount_ = 0; | 
| 420 | 
  | 
        velocity.accumulator.push_back( new VectorAccumulator() ); | 
| 421 | 
  | 
      data_[VELOCITY] = velocity; | 
| 422 | 
  | 
      outputMap_["VELOCITY"] = VELOCITY; | 
| 423 | 
+ | 
 | 
| 424 | 
+ | 
      OutputData angularVelocity; | 
| 425 | 
+ | 
      angularVelocity.units = "angstroms^2/fs"; | 
| 426 | 
+ | 
      angularVelocity.title =  "AngularVelocity";   | 
| 427 | 
+ | 
      angularVelocity.dataType = "Vector3d"; | 
| 428 | 
+ | 
      angularVelocity.accumulator.reserve(nBins_); | 
| 429 | 
+ | 
      for (int i = 0; i < nBins_; i++)  | 
| 430 | 
+ | 
        angularVelocity.accumulator.push_back( new VectorAccumulator() ); | 
| 431 | 
+ | 
      data_[ANGULARVELOCITY] = angularVelocity; | 
| 432 | 
+ | 
      outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY; | 
| 433 | 
  | 
 | 
| 434 | 
  | 
      OutputData density; | 
| 435 | 
  | 
      density.units =  "g cm^-3"; | 
| 557 | 
  | 
            slabWidth_ = hmat(2,2) / 10.0; | 
| 558 | 
  | 
         | 
| 559 | 
  | 
          if (hasSlabBCenter)  | 
| 560 | 
< | 
            slabBCenter_ = rnemdParams->getSlabACenter(); | 
| 560 | 
> | 
            slabBCenter_ = rnemdParams->getSlabBCenter(); | 
| 561 | 
  | 
          else  | 
| 562 | 
  | 
            slabBCenter_ = hmat(2,2) / 2.0; | 
| 563 | 
  | 
         | 
| 578 | 
  | 
        } | 
| 579 | 
  | 
      } | 
| 580 | 
  | 
    } | 
| 581 | 
+ | 
 | 
| 582 | 
  | 
    // object evaluator: | 
| 583 | 
  | 
    evaluator_.loadScriptString(rnemdObjectSelection_); | 
| 584 | 
  | 
    seleMan_.setSelectionSet(evaluator_.evaluate()); | 
| 573 | 
– | 
     | 
| 585 | 
  | 
    evaluatorA_.loadScriptString(selectionA_); | 
| 586 | 
  | 
    evaluatorB_.loadScriptString(selectionB_); | 
| 576 | 
– | 
     | 
| 587 | 
  | 
    seleManA_.setSelectionSet(evaluatorA_.evaluate()); | 
| 588 | 
  | 
    seleManB_.setSelectionSet(evaluatorB_.evaluate()); | 
| 579 | 
– | 
     | 
| 589 | 
  | 
    commonA_ = seleManA_ & seleMan_; | 
| 590 | 
< | 
    commonB_ = seleManB_ & seleMan_;     | 
| 590 | 
> | 
    commonB_ = seleManB_ & seleMan_;   | 
| 591 | 
  | 
  } | 
| 592 | 
  | 
   | 
| 593 | 
  | 
     | 
| 604 | 
  | 
#ifdef IS_MPI | 
| 605 | 
  | 
    } | 
| 606 | 
  | 
#endif | 
| 607 | 
+ | 
 | 
| 608 | 
+ | 
    // delete all of the objects we created: | 
| 609 | 
+ | 
    delete areaAccumulator_;     | 
| 610 | 
+ | 
    data_.clear(); | 
| 611 | 
  | 
  } | 
| 612 | 
  | 
   | 
| 613 | 
  | 
  void RNEMD::doSwap(SelectionManager& smanA, SelectionManager& smanB) { | 
| 1151 | 
  | 
          //if w is in the right range, so should be x, y, z. | 
| 1152 | 
  | 
          vector<StuntDouble*>::iterator sdi; | 
| 1153 | 
  | 
          Vector3d vel; | 
| 1154 | 
< | 
          for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { | 
| 1154 | 
> | 
          for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) { | 
| 1155 | 
  | 
            if (rnemdFluxType_ == rnemdFullKE) { | 
| 1156 | 
  | 
              vel = (*sdi)->getVel() * c; | 
| 1157 | 
  | 
              (*sdi)->setVel(vel); | 
| 1162 | 
  | 
            } | 
| 1163 | 
  | 
          } | 
| 1164 | 
  | 
          w = sqrt(w); | 
| 1165 | 
< | 
          for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { | 
| 1165 | 
> | 
          for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) { | 
| 1166 | 
  | 
            if (rnemdFluxType_ == rnemdFullKE) { | 
| 1167 | 
  | 
              vel = (*sdi)->getVel(); | 
| 1168 | 
  | 
              vel.x() *= x; | 
| 1281 | 
  | 
      vector<RealType>::iterator ri; | 
| 1282 | 
  | 
      RealType r1, r2, alpha0; | 
| 1283 | 
  | 
      vector<pair<RealType,RealType> > rps; | 
| 1284 | 
< | 
      for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) { | 
| 1284 | 
> | 
      for (ri = realRoots.begin(); ri !=realRoots.end(); ++ri) { | 
| 1285 | 
  | 
        r2 = *ri; | 
| 1286 | 
  | 
        //check if FindRealRoots() give the right answer | 
| 1287 | 
  | 
        if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) { | 
| 1313 | 
  | 
        RealType diff; | 
| 1314 | 
  | 
        pair<RealType,RealType> bestPair = make_pair(1.0, 1.0); | 
| 1315 | 
  | 
        vector<pair<RealType,RealType> >::iterator rpi; | 
| 1316 | 
< | 
        for (rpi = rps.begin(); rpi != rps.end(); rpi++) { | 
| 1316 | 
> | 
        for (rpi = rps.begin(); rpi != rps.end(); ++rpi) { | 
| 1317 | 
  | 
          r1 = (*rpi).first; | 
| 1318 | 
  | 
          r2 = (*rpi).second; | 
| 1319 | 
  | 
          switch(rnemdFluxType_) { | 
| 1380 | 
  | 
        } | 
| 1381 | 
  | 
        vector<StuntDouble*>::iterator sdi; | 
| 1382 | 
  | 
        Vector3d vel; | 
| 1383 | 
< | 
        for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { | 
| 1383 | 
> | 
        for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) { | 
| 1384 | 
  | 
          vel = (*sdi)->getVel(); | 
| 1385 | 
  | 
          vel.x() *= x; | 
| 1386 | 
  | 
          vel.y() *= y; | 
| 1391 | 
  | 
        x = 1.0 + px * (1.0 - x); | 
| 1392 | 
  | 
        y = 1.0 + py * (1.0 - y); | 
| 1393 | 
  | 
        z = 1.0 + pz * (1.0 - z); | 
| 1394 | 
< | 
        for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { | 
| 1394 | 
> | 
        for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) { | 
| 1395 | 
  | 
          vel = (*sdi)->getVel(); | 
| 1396 | 
  | 
          vel.x() *= x; | 
| 1397 | 
  | 
          vel.y() *= y; | 
| 1448 | 
  | 
    RealType Mc = 0.0; | 
| 1449 | 
  | 
    Mat3x3d Ic(0.0); | 
| 1450 | 
  | 
    RealType Kc = 0.0; | 
| 1451 | 
+ | 
 | 
| 1452 | 
+ | 
    // Constraints can be on only the linear or angular momentum, but | 
| 1453 | 
+ | 
    // not both.  Usually, the user will specify which they want, but | 
| 1454 | 
+ | 
    // in case they don't, the use of periodic boundaries should make | 
| 1455 | 
+ | 
    // the choice for us. | 
| 1456 | 
+ | 
    bool doLinearPart = false; | 
| 1457 | 
+ | 
    bool doAngularPart = false; | 
| 1458 | 
+ | 
 | 
| 1459 | 
+ | 
    switch (rnemdFluxType_) { | 
| 1460 | 
+ | 
    case rnemdPx: | 
| 1461 | 
+ | 
    case rnemdPy: | 
| 1462 | 
+ | 
    case rnemdPz: | 
| 1463 | 
+ | 
    case rnemdPvector: | 
| 1464 | 
+ | 
    case rnemdKePx: | 
| 1465 | 
+ | 
    case rnemdKePy: | 
| 1466 | 
+ | 
    case rnemdKePvector: | 
| 1467 | 
+ | 
      doLinearPart = true; | 
| 1468 | 
+ | 
      break; | 
| 1469 | 
+ | 
    case rnemdLx: | 
| 1470 | 
+ | 
    case rnemdLy: | 
| 1471 | 
+ | 
    case rnemdLz: | 
| 1472 | 
+ | 
    case rnemdLvector: | 
| 1473 | 
+ | 
    case rnemdKeLx: | 
| 1474 | 
+ | 
    case rnemdKeLy: | 
| 1475 | 
+ | 
    case rnemdKeLz: | 
| 1476 | 
+ | 
    case rnemdKeLvector: | 
| 1477 | 
+ | 
      doAngularPart = true; | 
| 1478 | 
+ | 
      break; | 
| 1479 | 
+ | 
    case rnemdKE: | 
| 1480 | 
+ | 
    case rnemdRotKE: | 
| 1481 | 
+ | 
    case rnemdFullKE: | 
| 1482 | 
+ | 
    default: | 
| 1483 | 
+ | 
      if (usePeriodicBoundaryConditions_)  | 
| 1484 | 
+ | 
        doLinearPart = true; | 
| 1485 | 
+ | 
      else | 
| 1486 | 
+ | 
        doAngularPart = true; | 
| 1487 | 
+ | 
      break; | 
| 1488 | 
+ | 
    } | 
| 1489 | 
  | 
     | 
| 1490 | 
  | 
    for (sd = smanA.beginSelected(selei); sd != NULL;  | 
| 1491 | 
  | 
         sd = smanA.nextSelected(selei)) { | 
| 1594 | 
  | 
                              MPI::REALTYPE, MPI::SUM); | 
| 1595 | 
  | 
#endif | 
| 1596 | 
  | 
     | 
| 1597 | 
+ | 
 | 
| 1598 | 
+ | 
    Vector3d ac, acrec, bc, bcrec; | 
| 1599 | 
+ | 
    Vector3d ah, ahrec, bh, bhrec; | 
| 1600 | 
+ | 
 | 
| 1601 | 
  | 
    bool successfulExchange = false; | 
| 1602 | 
  | 
    if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty | 
| 1603 | 
  | 
      Vector3d vc = Pc / Mc; | 
| 1604 | 
< | 
      Vector3d ac = -momentumTarget_ / Mc + vc; | 
| 1605 | 
< | 
      Vector3d acrec = -momentumTarget_ / Mc; | 
| 1604 | 
> | 
      ac = -momentumTarget_ / Mc + vc; | 
| 1605 | 
> | 
      acrec = -momentumTarget_ / Mc; | 
| 1606 | 
  | 
       | 
| 1607 | 
  | 
      // We now need the inverse of the inertia tensor to calculate the  | 
| 1608 | 
  | 
      // angular velocity of the cold slab; | 
| 1609 | 
  | 
      Mat3x3d Ici = Ic.inverse(); | 
| 1610 | 
  | 
      Vector3d omegac = Ici * Lc; | 
| 1611 | 
< | 
      Vector3d bc  = -(Ici * angularMomentumTarget_) + omegac; | 
| 1612 | 
< | 
      Vector3d bcrec = bc - omegac; | 
| 1611 | 
> | 
      bc  = -(Ici * angularMomentumTarget_) + omegac; | 
| 1612 | 
> | 
      bcrec = bc - omegac; | 
| 1613 | 
  | 
       | 
| 1614 | 
< | 
      RealType cNumerator = Kc - kineticTarget_  | 
| 1615 | 
< | 
        - 0.5 * Mc * ac.lengthSquare() - 0.5 * ( dot(bc, Ic * bc)); | 
| 1614 | 
> | 
      RealType cNumerator = Kc - kineticTarget_; | 
| 1615 | 
> | 
      if (doLinearPart)  | 
| 1616 | 
> | 
        cNumerator -= 0.5 * Mc * ac.lengthSquare(); | 
| 1617 | 
> | 
       | 
| 1618 | 
> | 
      if (doAngularPart) | 
| 1619 | 
> | 
        cNumerator -= 0.5 * ( dot(bc, Ic * bc)); | 
| 1620 | 
> | 
 | 
| 1621 | 
  | 
      if (cNumerator > 0.0) { | 
| 1622 | 
  | 
         | 
| 1623 | 
< | 
        RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare()  | 
| 1624 | 
< | 
          - 0.5*(dot(omegac, Ic * omegac)); | 
| 1623 | 
> | 
        RealType cDenominator = Kc; | 
| 1624 | 
> | 
 | 
| 1625 | 
> | 
        if (doLinearPart) | 
| 1626 | 
> | 
          cDenominator -= 0.5 * Mc * vc.lengthSquare(); | 
| 1627 | 
> | 
 | 
| 1628 | 
> | 
        if (doAngularPart) | 
| 1629 | 
> | 
          cDenominator -= 0.5*(dot(omegac, Ic * omegac)); | 
| 1630 | 
  | 
         | 
| 1631 | 
  | 
        if (cDenominator > 0.0) { | 
| 1632 | 
  | 
          RealType c = sqrt(cNumerator / cDenominator); | 
| 1633 | 
  | 
          if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients | 
| 1634 | 
  | 
             | 
| 1635 | 
  | 
            Vector3d vh = Ph / Mh; | 
| 1636 | 
< | 
            Vector3d ah = momentumTarget_ / Mh + vh; | 
| 1637 | 
< | 
            Vector3d ahrec = momentumTarget_ / Mh; | 
| 1636 | 
> | 
            ah = momentumTarget_ / Mh + vh; | 
| 1637 | 
> | 
            ahrec = momentumTarget_ / Mh; | 
| 1638 | 
  | 
             | 
| 1639 | 
  | 
            // We now need the inverse of the inertia tensor to | 
| 1640 | 
  | 
            // calculate the angular velocity of the hot slab; | 
| 1641 | 
  | 
            Mat3x3d Ihi = Ih.inverse(); | 
| 1642 | 
  | 
            Vector3d omegah = Ihi * Lh; | 
| 1643 | 
< | 
            Vector3d bh  = (Ihi * angularMomentumTarget_) + omegah; | 
| 1644 | 
< | 
            Vector3d bhrec = bh - omegah; | 
| 1643 | 
> | 
            bh  = (Ihi * angularMomentumTarget_) + omegah; | 
| 1644 | 
> | 
            bhrec = bh - omegah; | 
| 1645 | 
> | 
             | 
| 1646 | 
> | 
            RealType hNumerator = Kh + kineticTarget_; | 
| 1647 | 
> | 
            if (doLinearPart)  | 
| 1648 | 
> | 
              hNumerator -= 0.5 * Mh * ah.lengthSquare(); | 
| 1649 | 
  | 
             | 
| 1650 | 
< | 
            RealType hNumerator = Kh + kineticTarget_ | 
| 1651 | 
< | 
              - 0.5 * Mh * ah.lengthSquare() - 0.5 * ( dot(bh, Ih * bh));; | 
| 1650 | 
> | 
            if (doAngularPart) | 
| 1651 | 
> | 
              hNumerator -= 0.5 * ( dot(bh, Ih * bh)); | 
| 1652 | 
> | 
               | 
| 1653 | 
  | 
            if (hNumerator > 0.0) { | 
| 1654 | 
  | 
               | 
| 1655 | 
< | 
              RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare()  | 
| 1656 | 
< | 
                - 0.5*(dot(omegah, Ih * omegah)); | 
| 1655 | 
> | 
              RealType hDenominator = Kh; | 
| 1656 | 
> | 
              if (doLinearPart)  | 
| 1657 | 
> | 
                hDenominator -= 0.5 * Mh * vh.lengthSquare(); | 
| 1658 | 
> | 
              if (doAngularPart) | 
| 1659 | 
> | 
                hDenominator -= 0.5*(dot(omegah, Ih * omegah)); | 
| 1660 | 
  | 
               | 
| 1661 | 
  | 
              if (hDenominator > 0.0) { | 
| 1662 | 
  | 
                RealType h = sqrt(hNumerator / hDenominator); | 
| 1666 | 
  | 
                  Vector3d vel; | 
| 1667 | 
  | 
                  Vector3d rPos; | 
| 1668 | 
  | 
                   | 
| 1669 | 
< | 
                  for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { | 
| 1669 | 
> | 
                  for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) { | 
| 1670 | 
  | 
                    //vel = (*sdi)->getVel(); | 
| 1671 | 
  | 
                    rPos = (*sdi)->getPos() - coordinateOrigin_; | 
| 1672 | 
< | 
                    vel = ((*sdi)->getVel() - vc - cross(omegac, rPos)) * c  | 
| 1673 | 
< | 
                      + ac + cross(bc, rPos); | 
| 1672 | 
> | 
                    if (doLinearPart) | 
| 1673 | 
> | 
                      vel = ((*sdi)->getVel() - vc) * c + ac; | 
| 1674 | 
> | 
                    if (doAngularPart) | 
| 1675 | 
> | 
                      vel = ((*sdi)->getVel() - cross(omegac, rPos)) * c + cross(bc, rPos); | 
| 1676 | 
> | 
 | 
| 1677 | 
  | 
                    (*sdi)->setVel(vel); | 
| 1678 | 
  | 
                    if (rnemdFluxType_ == rnemdFullKE) { | 
| 1679 | 
  | 
                      if ((*sdi)->isDirectional()) { | 
| 1682 | 
  | 
                      } | 
| 1683 | 
  | 
                    } | 
| 1684 | 
  | 
                  } | 
| 1685 | 
< | 
                  for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { | 
| 1685 | 
> | 
                  for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) { | 
| 1686 | 
  | 
                    //vel = (*sdi)->getVel(); | 
| 1687 | 
  | 
                    rPos = (*sdi)->getPos() - coordinateOrigin_; | 
| 1688 | 
< | 
                    vel = ((*sdi)->getVel() - vh - cross(omegah, rPos)) * h  | 
| 1689 | 
< | 
                      + ah + cross(bh, rPos);                    | 
| 1688 | 
> | 
                    if (doLinearPart) | 
| 1689 | 
> | 
                      vel = ((*sdi)->getVel() - vh) * h + ah;      | 
| 1690 | 
> | 
                    if (doAngularPart) | 
| 1691 | 
> | 
                      vel = ((*sdi)->getVel() - cross(omegah, rPos)) * h + cross(bh, rPos);      | 
| 1692 | 
> | 
 | 
| 1693 | 
  | 
                    (*sdi)->setVel(vel); | 
| 1694 | 
  | 
                    if (rnemdFluxType_ == rnemdFullKE) { | 
| 1695 | 
  | 
                      if ((*sdi)->isDirectional()) { | 
| 1732 | 
  | 
      int isd; | 
| 1733 | 
  | 
      StuntDouble* sd; | 
| 1734 | 
  | 
      vector<StuntDouble*> aSites; | 
| 1656 | 
– | 
      ConvexHull* surfaceMeshA = new ConvexHull(); | 
| 1735 | 
  | 
      seleManA_.setSelectionSet(evaluatorA_.evaluate()); | 
| 1736 | 
  | 
      for (sd = seleManA_.beginSelected(isd); sd != NULL;  | 
| 1737 | 
  | 
           sd = seleManA_.nextSelected(isd)) { | 
| 1738 | 
  | 
        aSites.push_back(sd); | 
| 1739 | 
  | 
      } | 
| 1740 | 
+ | 
#if defined(HAVE_QHULL) | 
| 1741 | 
+ | 
      ConvexHull* surfaceMeshA = new ConvexHull(); | 
| 1742 | 
  | 
      surfaceMeshA->computeHull(aSites); | 
| 1743 | 
  | 
      areaA = surfaceMeshA->getArea(); | 
| 1744 | 
+ | 
      delete surfaceMeshA; | 
| 1745 | 
+ | 
#else | 
| 1746 | 
+ | 
      sprintf( painCave.errMsg, | 
| 1747 | 
+ | 
               "RNEMD::getDividingArea : Hull calculation is not possible\n" | 
| 1748 | 
+ | 
               "\twithout libqhull. Please rebuild OpenMD with qhull enabled."); | 
| 1749 | 
+ | 
      painCave.severity = OPENMD_ERROR; | 
| 1750 | 
+ | 
      painCave.isFatal = 1; | 
| 1751 | 
+ | 
      simError(); | 
| 1752 | 
+ | 
#endif | 
| 1753 | 
+ | 
 | 
| 1754 | 
  | 
    } else { | 
| 1755 | 
  | 
      if (usePeriodicBoundaryConditions_) { | 
| 1756 | 
  | 
        // in periodic boundaries, the surface area is twice the x-y | 
| 1768 | 
  | 
      int isd; | 
| 1769 | 
  | 
      StuntDouble* sd; | 
| 1770 | 
  | 
      vector<StuntDouble*> bSites; | 
| 1681 | 
– | 
      ConvexHull* surfaceMeshB = new ConvexHull(); | 
| 1771 | 
  | 
      seleManB_.setSelectionSet(evaluatorB_.evaluate()); | 
| 1772 | 
  | 
      for (sd = seleManB_.beginSelected(isd); sd != NULL;  | 
| 1773 | 
  | 
           sd = seleManB_.nextSelected(isd)) { | 
| 1774 | 
  | 
        bSites.push_back(sd); | 
| 1775 | 
  | 
      } | 
| 1776 | 
+ | 
 | 
| 1777 | 
+ | 
#if defined(HAVE_QHULL) | 
| 1778 | 
+ | 
      ConvexHull* surfaceMeshB = new ConvexHull();     | 
| 1779 | 
  | 
      surfaceMeshB->computeHull(bSites); | 
| 1780 | 
  | 
      areaB = surfaceMeshB->getArea(); | 
| 1781 | 
+ | 
      delete surfaceMeshB; | 
| 1782 | 
+ | 
#else | 
| 1783 | 
+ | 
      sprintf( painCave.errMsg, | 
| 1784 | 
+ | 
               "RNEMD::getDividingArea : Hull calculation is not possible\n" | 
| 1785 | 
+ | 
               "\twithout libqhull. Please rebuild OpenMD with qhull enabled."); | 
| 1786 | 
+ | 
      painCave.severity = OPENMD_ERROR; | 
| 1787 | 
+ | 
      painCave.isFatal = 1; | 
| 1788 | 
+ | 
      simError(); | 
| 1789 | 
+ | 
#endif | 
| 1790 | 
+ | 
 | 
| 1791 | 
+ | 
 | 
| 1792 | 
  | 
    } else { | 
| 1793 | 
  | 
      if (usePeriodicBoundaryConditions_) { | 
| 1794 | 
  | 
        // in periodic boundaries, the surface area is twice the x-y | 
| 1854 | 
  | 
  void RNEMD::collectData() { | 
| 1855 | 
  | 
    if (!doRNEMD_) return; | 
| 1856 | 
  | 
    Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); | 
| 1857 | 
< | 
 | 
| 1857 | 
> | 
     | 
| 1858 | 
  | 
    // collectData can be called more frequently than the doRNEMD, so use the  | 
| 1859 | 
  | 
    // computed area from the last exchange time: | 
| 1860 | 
< | 
 | 
| 1861 | 
< | 
    areaAccumulator_->add(getDividingArea()); | 
| 1860 | 
> | 
    RealType area = getDividingArea(); | 
| 1861 | 
> | 
    areaAccumulator_->add(area); | 
| 1862 | 
  | 
    Mat3x3d hmat = currentSnap_->getHmat(); | 
| 1863 | 
  | 
    seleMan_.setSelectionSet(evaluator_.evaluate()); | 
| 1864 | 
  | 
 | 
| 1917 | 
  | 
      Vector3d rPos = sd->getPos() - coordinateOrigin_; | 
| 1918 | 
  | 
      Vector3d aVel = cross(rPos, vel); | 
| 1919 | 
  | 
       | 
| 1920 | 
< | 
      if (binNo < nBins_)  { | 
| 1920 | 
> | 
      if (binNo >= 0 && binNo < nBins_)  { | 
| 1921 | 
  | 
        binCount[binNo]++; | 
| 1922 | 
  | 
        binMass[binNo] += mass; | 
| 1923 | 
  | 
        binPx[binNo] += mass*vel.x(); | 
| 1993 | 
  | 
      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]; | 
| 1997 | 
< | 
      aVel.y() = binOmegay[i]; | 
| 1998 | 
< | 
      aVel.z() = binOmegaz[i]; | 
| 1996 | 
> | 
      aVel.x() = binOmegax[i] / binCount[i]; | 
| 1997 | 
> | 
      aVel.y() = binOmegay[i] / binCount[i]; | 
| 1998 | 
> | 
      aVel.z() = binOmegaz[i] / binCount[i]; | 
| 1999 | 
  | 
 | 
| 2000 | 
  | 
      if (binCount[i] > 0) { | 
| 2001 | 
  | 
        // only add values if there are things to add | 
| 2017 | 
  | 
            case VELOCITY: | 
| 2018 | 
  | 
              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel); | 
| 2019 | 
  | 
              break; | 
| 2020 | 
< | 
            case ANGULARVELOCITY: | 
| 2020 | 
> | 
            case ANGULARVELOCITY:   | 
| 2021 | 
  | 
              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel); | 
| 2022 | 
  | 
              break; | 
| 2023 | 
  | 
            case DENSITY: | 
| 2028 | 
  | 
        } | 
| 2029 | 
  | 
      } | 
| 2030 | 
  | 
    } | 
| 2031 | 
+ | 
    hasData_ = true; | 
| 2032 | 
  | 
  } | 
| 2033 | 
  | 
 | 
| 2034 | 
  | 
  void RNEMD::getStarted() { | 
| 2061 | 
  | 
   | 
| 2062 | 
  | 
  void RNEMD::writeOutputFile() { | 
| 2063 | 
  | 
    if (!doRNEMD_) return; | 
| 2064 | 
+ | 
    if (!hasData_) return; | 
| 2065 | 
  | 
     | 
| 2066 | 
  | 
#ifdef IS_MPI | 
| 2067 | 
  | 
    // If we're the root node, should we print out the results | 
| 2083 | 
  | 
      RealType time = currentSnap_->getTime(); | 
| 2084 | 
  | 
      RealType avgArea; | 
| 2085 | 
  | 
      areaAccumulator_->getAverage(avgArea); | 
| 1981 | 
– | 
      RealType Jz = kineticExchange_ / (time * avgArea)  | 
| 1982 | 
– | 
        / PhysicalConstants::energyConvert; | 
| 1983 | 
– | 
      Vector3d JzP = momentumExchange_ / (time * avgArea);       | 
| 1984 | 
– | 
      Vector3d JzL = angularMomentumExchange_ / (time * avgArea);       | 
| 2086 | 
  | 
 | 
| 2087 | 
+ | 
      RealType Jz(0.0); | 
| 2088 | 
+ | 
      Vector3d JzP(V3Zero); | 
| 2089 | 
+ | 
      Vector3d JzL(V3Zero); | 
| 2090 | 
+ | 
      if (time >= info_->getSimParams()->getDt()) { | 
| 2091 | 
+ | 
        Jz = kineticExchange_ / (time * avgArea) | 
| 2092 | 
+ | 
          / PhysicalConstants::energyConvert; | 
| 2093 | 
+ | 
        JzP = momentumExchange_ / (time * avgArea); | 
| 2094 | 
+ | 
        JzL = angularMomentumExchange_ / (time * avgArea); | 
| 2095 | 
+ | 
      } | 
| 2096 | 
+ | 
 | 
| 2097 | 
  | 
      rnemdFile_ << "#######################################################\n"; | 
| 2098 | 
  | 
      rnemdFile_ << "# RNEMD {\n"; | 
| 2099 | 
  | 
 | 
| 2180 | 
  | 
          if (outputMask_[i]) { | 
| 2181 | 
  | 
            if (data_[i].dataType == "RealType") | 
| 2182 | 
  | 
              writeReal(i,j); | 
| 2183 | 
< | 
            else if (data_[i].dataType == "Vector3d") | 
| 2183 | 
> | 
            else if (data_[i].dataType == "Vector3d")  | 
| 2184 | 
  | 
              writeVector(i,j); | 
| 2185 | 
  | 
            else { | 
| 2186 | 
  | 
              sprintf( painCave.errMsg, | 
| 2237 | 
  | 
    RealType s; | 
| 2238 | 
  | 
    int count; | 
| 2239 | 
  | 
     | 
| 2240 | 
< | 
    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count(); | 
| 2240 | 
> | 
    count = data_[index].accumulator[bin]->count(); | 
| 2241 | 
  | 
    if (count == 0) return; | 
| 2242 | 
  | 
     | 
| 2243 | 
  | 
    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s); | 
| 2246 | 
  | 
      rnemdFile_ << "\t" << s; | 
| 2247 | 
  | 
    } else{ | 
| 2248 | 
  | 
      sprintf( painCave.errMsg, | 
| 2249 | 
< | 
               "RNEMD detected a numerical error writing: %s for bin %d", | 
| 2249 | 
> | 
               "RNEMD detected a numerical error writing: %s for bin %u", | 
| 2250 | 
  | 
               data_[index].title.c_str(), bin); | 
| 2251 | 
  | 
      painCave.isFatal = 1; | 
| 2252 | 
  | 
      simError(); | 
| 2260 | 
  | 
    Vector3d s; | 
| 2261 | 
  | 
    int count; | 
| 2262 | 
  | 
     | 
| 2263 | 
< | 
    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count(); | 
| 2263 | 
> | 
    count = data_[index].accumulator[bin]->count(); | 
| 2264 | 
> | 
 | 
| 2265 | 
  | 
    if (count == 0) return; | 
| 2266 | 
  | 
 | 
| 2267 | 
  | 
    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s); | 
| 2269 | 
  | 
        isinf(s[1]) || isnan(s[1]) ||  | 
| 2270 | 
  | 
        isinf(s[2]) || isnan(s[2]) ) {       | 
| 2271 | 
  | 
      sprintf( painCave.errMsg, | 
| 2272 | 
< | 
               "RNEMD detected a numerical error writing: %s for bin %d", | 
| 2272 | 
> | 
               "RNEMD detected a numerical error writing: %s for bin %u", | 
| 2273 | 
  | 
               data_[index].title.c_str(), bin); | 
| 2274 | 
  | 
      painCave.isFatal = 1; | 
| 2275 | 
  | 
      simError(); | 
| 2285 | 
  | 
    RealType s; | 
| 2286 | 
  | 
    int count; | 
| 2287 | 
  | 
     | 
| 2288 | 
< | 
    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count(); | 
| 2288 | 
> | 
    count = data_[index].accumulator[bin]->count(); | 
| 2289 | 
  | 
    if (count == 0) return; | 
| 2290 | 
  | 
     | 
| 2291 | 
  | 
    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s); | 
| 2294 | 
  | 
      rnemdFile_ << "\t" << s; | 
| 2295 | 
  | 
    } else{ | 
| 2296 | 
  | 
      sprintf( painCave.errMsg, | 
| 2297 | 
< | 
               "RNEMD detected a numerical error writing: %s std. dev. for bin %d", | 
| 2297 | 
> | 
               "RNEMD detected a numerical error writing: %s std. dev. for bin %u", | 
| 2298 | 
  | 
               data_[index].title.c_str(), bin); | 
| 2299 | 
  | 
      painCave.isFatal = 1; | 
| 2300 | 
  | 
      simError(); | 
| 2308 | 
  | 
    Vector3d s; | 
| 2309 | 
  | 
    int count; | 
| 2310 | 
  | 
     | 
| 2311 | 
< | 
    count = dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->count(); | 
| 2311 | 
> | 
    count = data_[index].accumulator[bin]->count(); | 
| 2312 | 
  | 
    if (count == 0) return; | 
| 2313 | 
  | 
 | 
| 2314 | 
  | 
    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s); | 
| 2316 | 
  | 
        isinf(s[1]) || isnan(s[1]) ||  | 
| 2317 | 
  | 
        isinf(s[2]) || isnan(s[2]) ) {       | 
| 2318 | 
  | 
      sprintf( painCave.errMsg, | 
| 2319 | 
< | 
               "RNEMD detected a numerical error writing: %s std. dev. for bin %d", | 
| 2319 | 
> | 
               "RNEMD detected a numerical error writing: %s std. dev. for bin %u", | 
| 2320 | 
  | 
               data_[index].title.c_str(), bin); | 
| 2321 | 
  | 
      painCave.isFatal = 1; | 
| 2322 | 
  | 
      simError(); |